Finding co-accessible networks with Cicero

In this vignette we will demonstrate how to find cis-co-accessible networks with Cicero using single-cell ATAC-seq data. Please see the Cicero website for information about Cicero.

To facilitate conversion between the Seurat (used by Signac) and CellDataSet (used by Cicero) formats, we will use a conversion function in the SeuratWrappers package available on GitHub.

Data loading

We will use a single-cell ATAC-seq dataset containing human CD34+ hematopoietic stem and progenitor cells published by Satpathy and Granja et al. (2019, Nature Biotechnology). The processed datasets are available on NCBI GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129785

This is the same dataset we used in the trajectory vignette, and we’ll start by loading the dataset that was created in that vignette. See the trajectory vignette for the code used to create the object from raw data.

First we will load their dataset and perform some standard preprocessing using Signac.

library(Signac)
library(Seurat)
library(SeuratWrappers)
library(ggplot2)
library(patchwork)

# load the object created in the Monocle 3 vignette
bone <- readRDS("cd34.rds")

Create the Cicero object

We can find cis-co-accessible networks (CCANs) using Cicero.

The Cicero developers have developed a separate branch of the package that works with a Monocle 3 CellDataSet object. We will first make sure this branch is installed, then convert our Seurat object for the whole bone marrow dataset to CellDataSet format.

# Install Cicero
if (!requireNamespace("remotes", quietly = TRUE))
    install.packages("remotes")
remotes::install_github("cole-trapnell-lab/cicero-release", ref = "monocle3")

## monocle3  (e658515b2... -> 98402ed0c...) [GitHub]
## openssl   (2.2.0        -> 2.2.1       ) [CRAN]
## tinytex   (0.51         -> 0.52        ) [CRAN]
## bslib     (0.7.0        -> 0.8.0       ) [CRAN]
## yaml      (2.3.8        -> 2.3.10      ) [CRAN]
## xfun      (0.44         -> 0.47        ) [CRAN]
## evaluate  (0.23         -> 0.24.0      ) [CRAN]
## checkmate (2.3.1        -> 2.3.2       ) [CRAN]
## XML       (3.99-0.16.1  -> 3.99-0.17   ) [CRAN]
## knitr     (1.45         -> 1.48        ) [CRAN]
## rmarkdown (2.27         -> 2.28        ) [CRAN]
## htmlTable (2.4.2        -> 2.4.3       ) [CRAN]
## DBI       (1.2.2        -> 1.2.3       ) [CRAN]
## Hmisc     (5.1-2        -> 5.1-3       ) [CRAN]
## 
## The downloaded binary packages are in
##  /var/folders/n5/1l8rdp0951gcwznhr_sxwj0m0000gp/T//Rtmp2ZNmZ2/downloaded_packages
## wk         (0.9.1    -> 0.9.2   ) [CRAN]
## s2         (1.1.6    -> 1.1.7   ) [CRAN]
## dqrng      (0.4.0    -> 0.4.1   ) [CRAN]
## RSpectra   (0.16-1   -> 0.16-2  ) [CRAN]
## biglm      (0.9-2.1  -> 0.9-3   ) [CRAN]
## spData     (2.3.0    -> 2.3.1   ) [CRAN]
## nloptr     (2.0.3    -> 2.1.1   ) [CRAN]
## minqa      (1.2.7    -> 1.2.8   ) [CRAN]
## uwot       (0.1.16   -> 0.2.2   ) [CRAN]
## spdep      (1.3-3    -> 1.3-5   ) [CRAN]
## slam       (0.1-50   -> 0.1-52  ) [CRAN]
## shiny      (1.8.1.1  -> 1.9.1   ) [CRAN]
## lme4       (1.1-35.3 -> 1.1-35.5) [CRAN]
## leidenbase (0.1.27   -> 0.1.30  ) [CRAN]
## 
## The downloaded binary packages are in
##  /var/folders/n5/1l8rdp0951gcwznhr_sxwj0m0000gp/T//Rtmp2ZNmZ2/downloaded_packages
## ── R CMD build ─────────────────────────────────────────────────────────────────
## * checking for file ‘/private/var/folders/n5/1l8rdp0951gcwznhr_sxwj0m0000gp/T/Rtmp2ZNmZ2/remotes480c3f5e994e/cole-trapnell-lab-monocle3-98402ed/DESCRIPTION’ ... OK
## * preparing ‘monocle3’:
## * checking DESCRIPTION meta-information ... OK
## * cleaning src
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## Omitted ‘LazyData’ from DESCRIPTION
## * building ‘monocle3_1.3.7.tar.gz’
## 
## ── R CMD build ─────────────────────────────────────────────────────────────────
## * checking for file ‘/private/var/folders/n5/1l8rdp0951gcwznhr_sxwj0m0000gp/T/Rtmp2ZNmZ2/remotes480c22118ef8/cole-trapnell-lab-cicero-release-2b0ee07/DESCRIPTION’ ... OK
## * preparing ‘cicero’:
## * checking DESCRIPTION meta-information ... OK
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## * looking to see if a ‘data/datalist’ file should be added
## * building ‘cicero_1.3.9.tar.gz’

library(cicero)

# convert to CellDataSet format and make the cicero object
bone.cds <- as.cell_data_set(x = bone)
bone.cicero <- make_cicero_cds(bone.cds, reduced_coordinates = reducedDims(bone.cds)$UMAP)

Find Cicero connections

We’ll demonstrate running Cicero here using just one chromosome to save some time, but the same workflow can be applied to find CCANs for the whole genome.

Here we demonstrate the most basic workflow for running Cicero. This workflow can be broken down into several steps, each with parameters that can be changed from their defaults to fine-tune the Cicero algorithm depending on your data. We highly recommend that you explore the Cicero website, paper, and documentation for more information.

# get the chromosome sizes from the Seurat object
genome <- seqlengths(bone)

# use chromosome 1 to save some time
# omit this step to run on the whole genome
genome <- genome[1]

# convert chromosome sizes to a dataframe
genome.df <- data.frame("chr" = names(genome), "length" = genome)

# run cicero
conns <- run_cicero(bone.cicero, genomic_coords = genome.df, sample_num = 100)

## [1] "Starting Cicero"
## [1] "Calculating distance_parameter value"
## [1] "Running models"
## [1] "Assembling connections"
## [1] "Successful cicero models:  755"
## [1] "Other models: "
## 
##   Too many elements in range Zero or one element in range 
##                          157                           86 
## [1] "Models with errors:  0"
## [1] "Done"

head(conns)

##                      Peak1                  Peak2 coaccess
## 1 chr1-100003337-100003837 chr1-99791719-99792219        0
## 2 chr1-100003337-100003837 chr1-99828699-99829199        0
## 3 chr1-100003337-100003837 chr1-99835542-99836042        0
## 4 chr1-100003337-100003837 chr1-99836217-99836717        0
## 5 chr1-100003337-100003837 chr1-99839576-99840076        0
## 6 chr1-100003337-100003837 chr1-99840640-99841140        0

Find cis-co-accessible networks (CCANs)

Now that we’ve found pairwise co-accessibility scores for each peak, we can now group these pairwise connections into larger co-accessible networks using the generate_ccans() function from Cicero.

ccans <- generate_ccans(conns)

## [1] "Coaccessibility cutoff used: 0.19"

head(ccans)

##                                              Peak CCAN
## chr1-10009702-10010202     chr1-10009702-10010202    1
## chr1-100151188-100151688 chr1-100151188-100151688    2
## chr1-100164787-100165287 chr1-100164787-100165287    2
## chr1-100165566-100166066 chr1-100165566-100166066    2
## chr1-100202505-100203005 chr1-100202505-100203005    3
## chr1-100215491-100215991 chr1-100215491-100215991    3

Add links to a Seurat object

We can add the co-accessible links found by Cicero to the ChromatinAssay object in Seurat. Using the ConnectionsToLinks() function in Signac we can convert the outputs of Cicero to the format needed to store in the links slot in the ChromatinAssay, and add this to the object using the Links<- assignment function.

links <- ConnectionsToLinks(conns = conns, ccans = ccans)
Links(bone) <- links

We can now visualize these links along with DNA accessibility information by running CoveragePlot() for a region:

CoveragePlot(bone, region = "chr1-40189344-40252549")

Acknowledgements

Thanks to the developers of Cicero, especially Cole Trapnell, Hannah Pliner, and members of the Trapnell lab. If you use Cicero please cite the Cicero paper.

Session Info

sessionInfo()