## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----getPackage, eval=FALSE---------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
#  BiocManager::install("cicero")

## ----eval = FALSE-------------------------------------------------------------
# BiocManager::install(cole-trapnell-lab/cicero)

## ----Load, message=FALSE------------------------------------------------------
 library(cicero)

## -----------------------------------------------------------------------------
data(cicero_data)

## ----eval=TRUE----------------------------------------------------------------
input_cds <- make_atac_cds(cicero_data, binarize = TRUE)

## ----eval=TRUE----------------------------------------------------------------
set.seed(2017)
input_cds <- detectGenes(input_cds)
input_cds <- estimateSizeFactors(input_cds)
input_cds <- reduceDimension(input_cds, max_components = 2, num_dim=6,
                        reduction_method = 'tSNE', norm_method = "none")

## ----eval=FALSE---------------------------------------------------------------
# tsne_coords <- t(reducedDimA(input_cds))
# row.names(tsne_coords) <- row.names(pData(input_cds))
# cicero_cds <- make_cicero_cds(input_cds, reduced_coordinates = tsne_coords)

## ----eval=FALSE---------------------------------------------------------------
# data("human.hg19.genome")
# sample_genome <- subset(human.hg19.genome, V1 == "chr18")
# sample_genome$V2[1] <- 10000000
# conns <- run_cicero(cicero_cds, sample_genome, sample_num = 2) # Takes a few minutes to run
# head(conns)

## ----fig.width = 7, fig.height = 4, fig.align='center', eval=FALSE------------
# data(gene_annotation_sample)
# plot_connections(conns, "chr18", 8575097, 8839855,
#                  gene_model = gene_annotation_sample,
#                  coaccess_cutoff = .25,
#                  connection_width = .5,
#                  collapseTranscripts = "longest" )

## ----eval=FALSE---------------------------------------------------------------
# chia_conns <-  data.frame(Peak1 = c("chr18_10000_10200", "chr18_10000_10200",
#                                     "chr18_49500_49600"),
#                           Peak2 = c("chr18_10600_10700", "chr18_111700_111800",
#                                     "chr18_10600_10700"))
# 
# conns$in_chia <- compare_connections(conns, chia_conns)

## ----eval=FALSE---------------------------------------------------------------
# conns$in_chia_100 <- compare_connections(conns, chia_conns, maxgap=100)
# 
# head(conns)

## ----eval=FALSE---------------------------------------------------------------
# # Add a column of 1s called "coaccess"
# chia_conns <-  data.frame(Peak1 = c("chr18_10000_10200", "chr18_10000_10200",
#                                     "chr18_49500_49600"),
#                           Peak2 = c("chr18_10600_10700", "chr18_111700_111800",
#                                     "chr18_10600_10700"),
#                           coaccess = c(1, 1, 1))
# 
# plot_connections(conns, "chr18", 10000, 112367,
#                  gene_model = gene_annotation_sample,
#                  coaccess_cutoff = 0,
#                  connection_width = .5,
#                  comparison_track = chia_conns,
#                  include_axis_track = FALSE,
#                  collapseTranscripts = "longest")

## ----eval=FALSE---------------------------------------------------------------
# CCAN_assigns <- generate_ccans(conns)
# 
# head(CCAN_assigns)

## ----eval=FALSE---------------------------------------------------------------
# 
# #### Add a column for the pData table indicating the gene if a peak is a promoter ####
# # Create a gene annotation set that only marks the transcription start sites of
# # the genes. We use this as a proxy for promoters.
# # To do this we need the first exon of each transcript
# pos <- subset(gene_annotation_sample, strand == "+")
# pos <- pos[order(pos$start),]
# pos <- pos[!duplicated(pos$transcript),] # remove all but the first exons per transcript
# pos$end <- pos$start + 1 # make a 1 base pair marker of the TSS
# 
# neg <- subset(gene_annotation_sample, strand == "-")
# neg <- neg[order(neg$start, decreasing = TRUE),]
# neg <- neg[!duplicated(neg$transcript),] # remove all but the first exons per transcript
# neg$start <- neg$end - 1
# 
# gene_annotation_sub <- rbind(pos, neg)
# 
# # Make a subset of the TSS annotation columns containing just the coordinates
# # and the gene name
# gene_annotation_sub <- gene_annotation_sub[,c(1:3, 8)]
# 
# # Rename the gene symbol column to "gene"
# names(gene_annotation_sub)[4] <- "gene"
# 
# input_cds <- annotate_cds_by_site(input_cds, gene_annotation_sub)
# 
# head(fData(input_cds))
# 
# 
# #### Generate gene activity scores ####
# # generate unnormalized gene activity matrix
# unnorm_ga <- build_gene_activity_matrix(input_cds, conns)
# 
# # remove any rows/columns with all zeroes
# unnorm_ga <- unnorm_ga[!Matrix::rowSums(unnorm_ga) == 0, !Matrix::colSums(unnorm_ga) == 0]
# 
# # make a list of num_genes_expressed
# num_genes <- pData(input_cds)$num_genes_expressed
# names(num_genes) <- row.names(pData(input_cds))
# 
# # normalize
# cicero_gene_activities <- normalize_gene_activities(unnorm_ga, num_genes)
# 
# # if you had two datasets to normalize, you would pass both:
# # num_genes should then include all cells from both sets
# unnorm_ga2 <- unnorm_ga
# cicero_gene_activities <- normalize_gene_activities(list(unnorm_ga, unnorm_ga2), num_genes)
# 

## ----eval=TRUE----------------------------------------------------------------
data("cicero_data")
input_cds <- make_atac_cds(cicero_data)

# Add some cell meta-data
data("cell_data")
pData(input_cds) <- cbind(pData(input_cds), cell_data[row.names(pData(input_cds)),])
pData(input_cds)$cell <- NULL

agg_cds <- aggregate_nearby_peaks(input_cds, distance = 10000)
agg_cds <- detectGenes(agg_cds)
agg_cds <- estimateSizeFactors(agg_cds)
agg_cds <- estimateDispersions(agg_cds)

## ----eval=TRUE----------------------------------------------------------------
# This takes a few minutes to run
diff_timepoint <- differentialGeneTest(agg_cds,
                      fullModelFormulaStr="~timepoint + num_genes_expressed")

# We chose a very high q-value cutoff because there are
# so few sites in the sample dataset, in general a q-value
# cutoff in the range of 0.01 to 0.1 would be appropriate
ordering_sites <- row.names(subset(diff_timepoint, qval < .5))
length(ordering_sites)

## ----fig.show='hold', eval=FALSE----------------------------------------------
# plot_pc_variance_explained(agg_cds, return_all = FALSE) #Choose 2 PCs
# agg_cds <- reduceDimension(agg_cds,
#                               max_components = 2,
#                               norm_method = 'log',
#                               num_dim = 2,
#                               reduction_method = 'tSNE',
#                               verbose = TRUE)
# 
# agg_cds <- clusterCells(agg_cds, verbose = FALSE)
# 
# plot_cell_clusters(agg_cds, color_by = 'as.factor(Cluster)') + theme(text = element_text(size=8))
# clustering_DA_sites <- differentialGeneTest(agg_cds[1:10,], #Subset for example only
#                                              fullModelFormulaStr = '~Cluster')
# 
# # Not run because using Option 1 to continue
# # ordering_sites <-
# #  row.names(clustering_DA_sites)[order(clustering_DA_sites$qval)][1:1000]
# 

## ----eval=TRUE----------------------------------------------------------------
agg_cds <- setOrderingFilter(agg_cds, ordering_sites)

## ----fig.align='center', fig.height=4, fig.width=4, eval=TRUE-----------------
agg_cds <- reduceDimension(agg_cds, max_components = 2,
          residualModelFormulaStr="~num_genes_expressed",
          reduction_method = 'DDRTree')

# skipped due to monocle bug
# agg_cds <- suppressWarnings(orderCells(agg_cds))

plot_cell_trajectory(agg_cds, color_by = "timepoint")

## ----fig.align='center', fig.height=4, fig.width=4, eval=TRUE-----------------
# skipped due to monocle bug
#plot_cell_trajectory(agg_cds, color_by = "State")

## ----fig.align='center', fig.height=4, fig.width=4, eval=TRUE-----------------
# skipped due to monocle bug
#agg_cds <- suppressWarnings(orderCells(agg_cds, root_state = 4))
#plot_cell_trajectory(agg_cds, color_by = "Pseudotime")

## ----fig.align='center', fig.height=4, fig.width=4, eval=TRUE-----------------
# skipped due to monocle bug
#pData(input_cds)$Pseudotime <- pData(agg_cds)[colnames(input_cds),]$Pseudotime
#pData(input_cds)$State <- pData(agg_cds)[colnames(input_cds),]$State

## ----fig.width = 3, fig.height = 4, fig.align='center', eval=TRUE-------------
# skipped due to monocle bug
#input_cds_lin <- input_cds[,row.names(subset(pData(input_cds), State  != 5))]

#plot_accessibility_in_pseudotime(input_cds_lin[c("chr18_38156577_38158261", 
#                                                 "chr18_48373358_48374180", 
#                                                 "chr18_60457956_60459080")])

## ----eval=TRUE----------------------------------------------------------------
# skipped due to monocle bug
#pData(input_cds_lin)$cell_subtype <- cut(pData(input_cds_lin)$Pseudotime, 10)
#binned_input_lin <- aggregate_by_cell_bin(input_cds_lin, "cell_subtype")

## ----eval=TRUE----------------------------------------------------------------
# skipped due to monocle bug
#diff_test_res <- differentialGeneTest(binned_input_lin[1:10,], #Subset for example only
#    fullModelFormulaStr="~sm.ns(Pseudotime, df=3) + sm.ns(num_genes_expressed, df=3)",
#    reducedModelFormulaStr="~sm.ns(num_genes_expressed, df=3)", cores=1)

#head(diff_test_res)

## -----------------------------------------------------------------------------
citation("cicero")

## -----------------------------------------------------------------------------
sessionInfo()

