In this tutorial, we’ll demonstrate how to infer gene regulatory network using single-cell multimodal data which measures gene expression and chromatin accessibility profiles from the same single cells. We will use a publicly available 10x Genomic Multiome data set for human PBMCs.

We first download the required data.

Run the following commands to download the data:

if(!dir.exists('./10x_pbmc')){
    dir.create('./10x_pbmc')
}

download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/10x_multiome/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz', 
              destfile = './10x_pbmc/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz', 
              method = 'wget', extra = '--no-check-certificate')

download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/10x_multiome/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz.tbi', 
               destfile = './10x_pbmc/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz.tbi', 
               method = 'wget', extra = '--no-check-certificate')
              
download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/10x_multiome/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5', 
              destfile = './10x_pbmc/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5', 
              method = 'wget', extra = '--no-check-certificate')

Next, we load all necessary packages:

Data pre-processing

Let’s load the data into memory and extract scRNA and scATAC-seq data:

inputdata.10x <- Read10X_h5("./10x_pbmc/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
# extract RNA and ATAC data
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# filter peaks by chromosome
atac_counts <- atac_counts[grep("chr", rownames(atac_counts)), ]

Next create a Seurat object and filter low-quality cells in scRNA-seq data

obj.rna <- CreateSeuratObject(counts = rna_counts)
obj.rna[["percent.mt"]] <- PercentageFeatureSet(obj.rna, pattern = "^MT-")

We can visualize the data quality

# Visualize QC metrics as a violin plot
VlnPlot(obj.rna, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0)

obj.rna <- subset(obj.rna, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mt < 20)

obj.rna
## An object of class Seurat 
## 36601 features across 10504 samples within 1 assay 
## Active assay: RNA (36601 features, 0 variable features)
##  2 layers present: counts, data

Create Seurat object for scATAC-seq data:

# create seurat object
chrom_assay <- CreateChromatinAssay(
    counts = atac_counts,
    sep = c(":", "-"),
    min.cells = 1,
    genome = 'hg38',
    fragments = './10x_pbmc/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz'
)
## Computing hash
obj.atac <- CreateSeuratObject(
    counts = chrom_assay,
    assay = "ATAC")

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86, 
                                   verbose = FALSE)
## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)

## Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
##   suppressWarnings() to suppress this warning.)
# change to UCSC style since the data was mapped to hg38
seqlevelsStyle(annotations) <- 'UCSC'

# add the gene information to the object
Annotation(obj.atac) <- annotations

Make sure that the two modalities have the same cells:

cell.sel <- intersect(colnames(obj.rna), colnames(obj.atac))

obj.rna <- subset(obj.rna, cells = cell.sel)
obj.atac <- subset(obj.atac, cells = cell.sel)

We next process the scRNA-seq and scATAC-seq data using standard Seurat and Signac analysis pipeline:

# normalization followed by dimensionality reduction
obj.rna <- obj.rna %>%
    NormalizeData(verbose=FALSE) %>%
    FindVariableFeatures(nfeatures=3000, verbose=F) %>%
    ScaleData() %>%
    RunPCA(verbose = FALSE) %>%
    RunUMAP(dims = 1:30, verbose = FALSE)
## Warning: The following arguments are not used: nselect
## Centering and scaling data matrix
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
obj.atac <- obj.atac %>%
    RunTFIDF() %>%
    FindTopFeatures() %>%
    RunSVD() %>%
    RunUMAP(reduction = 'lsi', dims = 2:30, verbose = FALSE)
## Performing TF-IDF normalization
## Running SVD
## Scaling cell embeddings
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'

Visualize the scRNA-seq and scATAC-seq separately.

p1 <- DimPlot(obj.rna) + NoLegend()
p2 <- DimPlot(obj.atac) + NoLegend()

p1 + p2

Cell type annotation

We next will annotate the cell types in our multimodal data via label transfer approach implemented by Seurat. The reference dataset is found here

Download data

download.file(url = 'https://costalab.ukaachen.de/open_data/scMEGA/10x_multiome/seurat.rds', 
              destfile = './10x_pbmc/seurat.rds', 
              method = 'wget', extra = '--no-check-certificate')

Load data

reference <- readRDS("./10x_pbmc/seurat.rds")

reference <- UpdateSeuratObject(object = reference)
## Validating object structure
## Updating object slots
## Ensuring keys are in the proper structure
## Warning: Assay RNA changing from Assay to Assay
## Warning: Assay ATAC changing from Assay to Assay
## Ensuring keys are in the proper structure
## Ensuring feature names don't have underscores or pipes
## Updating slots in RNA
## Updating slots in ATAC
## Validating object structure for Assay 'RNA'
## Validating object structure for Assay 'ATAC'
## Object representation is consistent with the most current Seurat version
reference[["RNA"]] <- as(reference[["RNA"]], "Assay5")
## Warning: Assay RNA changing from Assay to Assay5
# we only keep cells with annotated cell type
reference <- subset(reference, subset = celltype != "NA")
## Warning: Removing 1877 cells missing data for vars requested
# run sctransform
reference <- reference %>%
    NormalizeData(verbose=FALSE) %>%
    FindVariableFeatures(nfeatures=3000, verbose=F) %>%
    ScaleData() %>%
    RunPCA(verbose = FALSE) %>%
    RunUMAP(dims = 1:30, verbose = FALSE)
## Warning: The following arguments are not used: layer
## Centering and scaling data matrix
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
p1 <- DimPlot(reference, label = TRUE, repel = TRUE, 
        reduction = "umap", group.by = "broad_celltype") + NoLegend()

p2 <- DimPlot(reference, label = TRUE, repel = TRUE, 
        reduction = "umap", group.by = "celltype") + NoLegend()

p1 + p2

Transfer cell type labels from reference to query:

transfer_anchors <- FindTransferAnchors(
  reference = reference,
  query = obj.rna,
  reference.reduction = "pca",
  dims = 1:30
)
## Projecting cell embeddings
## Finding neighborhoods
## Finding anchors
##  Found 26315 anchors
predictions <- TransferData(
  anchorset = transfer_anchors, 
  refdata = reference$celltype,
  weight.reduction = obj.rna[['pca']],
  dims = 1:30
)
## Finding integration vectors
## Finding integration vector weights
## Predicting cell labels
obj.rna <- AddMetaData(
  object = obj.rna,
  metadata = predictions
)

obj.atac <- AddMetaData(
  object = obj.atac,
  metadata = predictions
)

Visualize the predicted cell types:

p1 <- DimPlot(obj.rna, label = TRUE, repel = TRUE, 
        reduction = "umap", group.by = "predicted.id") + NoLegend()

p2 <- DimPlot(obj.atac, label = TRUE, repel = TRUE, 
        reduction = "umap", group.by = "predicted.id") + NoLegend()

p1 + p2

Integration of multimodal single-cell data using MOJITOO

We next need to project the cells into a low-dimensional space. Here we can use MOJITOO to do this job.

Let’s first create a single Seurat object including both modalities and predicted cell types that we generated in the above step:

meta.data <- obj.rna@meta.data %>%
    as.data.frame()

# create a Seurat object containing the RNA adata
pbmc <- CreateSeuratObject(
  counts = obj.rna@assays$RNA@counts,
  assay = "RNA",
    meta.data = meta.data
)

# create ATAC assay and add it to the object
pbmc[["ATAC"]] <- CreateChromatinAssay(
  counts = obj.atac@assays$ATAC@counts,
  sep = c(":", "-"),
    min.cells = 1,
    genome = 'hg38',
    fragments = './10x_pbmc/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz'
)
## Computing hash

MOJITOO takes as input the individual low-dimensional space for each modality, so here we do this again:

## RNA pre-processing and PCA dimension reduction
DefaultAssay(pbmc) <- "RNA"

pbmc <- pbmc %>%
    NormalizeData(verbose=F) %>%
    FindVariableFeatures(nfeatures=3000, verbose=F) %>%
    ScaleData(verbose=F) %>%
    RunPCA(npcs=50, reduction.name="RNA_PCA", verbose=F)
## Warning: The following arguments are not used: nselect
## ATAC pre-processing and LSI dimension reduction
DefaultAssay(pbmc) <- "ATAC"

pbmc <- pbmc %>%
    RunTFIDF(verbose=F) %>%
    FindTopFeatures(min.cutoff = 'q0', verbose=F) %>%
    RunSVD(verbose=F)

Run MOJITOO:

pbmc <- mojitoo(
     object = pbmc,
     reduction.list = list("RNA_PCA", "lsi"),
     dims.list = list(1:50, 2:50), ## exclude 1st dimension of LSI
     reduction.name = 'MOJITOO',
     assay = "RNA"
)
## processing RNA_PCA
## adding lsi
## 1 round cc 45
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from MOJITOO to MOJITOO_

We can generate another UMAP representation based on the MOJITOO results:

DefaultAssay(pbmc) <- "RNA"
embedd <- Embeddings(pbmc[["MOJITOO"]])
pbmc <- RunUMAP(pbmc, 
                reduction="MOJITOO", 
                reduction.name="MOJITOO_UMAP", 
                dims=1:ncol(embedd), verbose=F)
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
## Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
## Also defined by 'spam'
DimPlot(pbmc, group.by = "predicted.id", 
        shuffle = TRUE, label = TRUE, reduction = "MOJITOO_UMAP") + NoLegend()

Trajectory analysis

We next infer a trajectory from naive CD4 T cells to memory CD4 T cells to characterize CD4+ T cell activation.

pbmc <- AddTrajectory(object = pbmc, 
                      trajectory = c("naive CD4 T cells", 
                                     "memory CD4 T cells"),
                      group.by = "predicted.id", 
                          reduction = "MOJITOO_UMAP",
                          dims = 1:2, 
                          use.all = FALSE)

# we only plot the cells that are in this trajectory
pbmc.t.cells <- pbmc[, !is.na(pbmc$Trajectory)]

The results can be visualized as:

p1 <- DimPlot(object = pbmc.t.cells, 
              group.by = "predicted.id", 
              reduction = "MOJITOO_UMAP",
             label = TRUE) + NoLegend()

p2 <- TrajectoryPlot(object = pbmc.t.cells, 
                    reduction = "MOJITOO_UMAP",
                    continuousSet = "blueYellow",
                    size = 1,
                   addArrow = FALSE)
p1 + p2

TF and gene selection

We next select candidate TFs and genes for building a meaningful gene regulatory network.

Select TFs

To identify potential regulator (i.e., TFs), we first estimate an acitivty score for each TF in each cell. This is done by first performing motif matching and then computing deviation scores using chromVAR.

# Get a list of motif position frequency matrices from the JASPAR database
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)

# add motif information
pbmc.t.cells <- AddMotifs(
  object = pbmc.t.cells,
  genome = BSgenome.Hsapiens.UCSC.hg38,
  pfm = pfm,
    assay = "ATAC"
)
## Building motif matrix
## Finding motif positions
## Creating Motif object
# run chromVAR
pbmc.t.cells <- RunChromVAR(
  object = pbmc.t.cells,
  genome = BSgenome.Hsapiens.UCSC.hg38,
    assay = "ATAC"
)
## Computing GC bias per region
## Selecting background regions
## Computing deviations from background
## Constructing chromVAR assay
sel.tfs <- SelectTFs(object = pbmc.t.cells, 
                 return.heatmap = TRUE,
                cor.cutoff = 0.4)
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
## Creating Trajectory Group Matrix..
## Smoothing...
## Find 559 shared features!
df.cor <- sel.tfs$tfs
ht <- sel.tfs$heatmap

draw(ht)

Select genes

sel.genes <- SelectGenes(object = pbmc.t.cells,
                  labelTop1 = 0,
                  labelTop2 = 0)
## Creating Trajectory Group Matrix..
## Smoothing...
## Creating Trajectory Group Matrix..
## Smoothing...
## Linking cis-regulatory elements to genes...
df.p2g <- sel.genes$p2g
ht <- sel.genes$heatmap

draw(ht)

Gene regulatory network inference and visualization

We here will try to predict a gene regulatory network based on the correlation of the TF binding activity as estimated from snATAC-seq and gene expression as measured by snRNA-seq along the trajectory.

tf.gene.cor <- GetTFGeneCorrelation(object = pbmc.t.cells, 
                                    tf.use = df.cor$tfs, 
                                    gene.use = unique(df.p2g$gene),
                                    tf.assay = "chromvar", 
                                    gene.assay = "RNA",
                                    trajectory.name = "Trajectory")
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
## Creating Trajectory Group Matrix..
## Smoothing...

We can then visualize this correlation matrix by heatmap. Also, we can cluster the genes and TFs to identify different regulatory modules for the predefined sub-populations.

ht <- GRNHeatmap(tf.gene.cor, 
                 tf.timepoint = df.cor$time_point)
## `use_raster` is automatically set to TRUE for a matrix with more than
## 2000 rows. You can control `use_raster` argument by explicitly setting
## TRUE/FALSE to it.
## 
## Set `ht_opt$message = FALSE` to turn off this message.
## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.
ht

To associate genes to TFs, we will use the peak-to-gene links and TF binding sites information. Specifically, if a gene is regulated by a peak and this peak is bound by a TF, then we consider this gene to be a target of this TF.

motif.matching <- pbmc.t.cells@assays$ATAC@motifs@data
colnames(motif.matching) <- pbmc.t.cells@assays$ATAC@motifs@motif.names
motif.matching <-
    motif.matching[unique(df.p2g$peak), unique(tf.gene.cor$tf)]


df.grn <- GetGRN(motif.matching = motif.matching, 
                 df.cor = tf.gene.cor, 
                 df.p2g = df.p2g)
## Filtering network by peak-to-gene links...
## Filtering network by TF binding site prediction...

Finally, we can visualize our network as the last step of this analysis

# define colors for nodes representing TFs (i.e., regulators)
df.cor <- df.cor[order(df.cor$time_point), ]
tfs.timepoint <- df.cor$time_point
names(tfs.timepoint) <- df.cor$tfs

# plot the graph, here we can highlight some genes
df.grn2 <- df.grn %>%
    subset(correlation > 0.5) %>%
    select(c(tf, gene, correlation)) %>%
    rename(weights = correlation)

p <- GRNPlot(df.grn2, 
             tfs.timepoint = tfs.timepoint,
             show.tf.labels = TRUE,
             seed = 42, 
             plot.importance = FALSE,
            min.importance = 2,
            remove.isolated = FALSE)

print(p)
## Warning: Using alpha for a discrete variable is not advised.

GRN visualization

Once we generated the gene regulatory network, we can visualize individual TFs in terms of binding activity, expression, and target expression along the pseudotime trajectory.

Here we select two TFs for visualization.

pbmc.t.cells <- AddTargetAssay(object = pbmc.t.cells, df.grn = df.grn2)

p1 <- PseudotimePlot(object = pbmc.t.cells, tf.use = "SOX4")
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
## Creating Trajectory Group Matrix..
## Smoothing...
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
p2 <- PseudotimePlot(object = pbmc.t.cells, tf.use = "TBX21")
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
## Creating Trajectory Group Matrix..
## Smoothing...
## Creating Trajectory Group Matrix..
## Some values are below 0, this could be the Motif activity matrix in which scaleTo should be set = NULL.
## Continuing without depth normalization!
## Smoothing...
p1 + p2
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

The x-axis in above plots present pseudotime point along the trajectory, and the y-axis represent TF binding acitivty, TF expression, and TF target expression after z-score transformation.

# Check session information
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /data/pinello/SHARED_SOFTWARE/anaconda_latest/envs/zl_scmega/lib/libopenblasp-r0.3.24.so;  LAPACK version 3.11.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] circlize_0.4.15                   ComplexHeatmap_2.16.0            
##  [3] MOJITOO_1.0                       ggraph_2.1.0                     
##  [5] igraph_1.5.1                      TFBSTools_1.38.0                 
##  [7] JASPAR2020_0.99.10                dplyr_1.1.3                      
##  [9] EnsDb.Hsapiens.v86_2.99.0         ensembldb_2.24.0                 
## [11] AnnotationFilter_1.24.0           GenomicFeatures_1.52.2           
## [13] AnnotationDbi_1.62.2              BSgenome.Hsapiens.UCSC.hg38_1.4.5
## [15] BSgenome_1.68.0                   rtracklayer_1.60.1               
## [17] Biostrings_2.68.1                 XVector_0.40.0                   
## [19] Nebulosa_1.10.0                   patchwork_1.1.3                  
## [21] scMEGA_1.0.2                      Signac_1.11.9000                 
## [23] Seurat_4.9.9.9060                 SeuratObject_4.9.9.9091          
## [25] sp_2.0-0                          rhdf5_2.44.0                     
## [27] SummarizedExperiment_1.30.2       Biobase_2.60.0                   
## [29] MatrixGenerics_1.12.3             Rcpp_1.0.11                      
## [31] Matrix_1.6-1                      GenomicRanges_1.52.0             
## [33] GenomeInfoDb_1.36.3               IRanges_2.34.1                   
## [35] S4Vectors_0.38.1                  BiocGenerics_0.46.0              
## [37] matrixStats_1.0.0                 data.table_1.14.8                
## [39] stringr_1.5.0                     plyr_1.8.8                       
## [41] magrittr_2.0.3                    ggplot2_3.4.3                    
## [43] gtable_0.3.4                      gtools_3.9.4                     
## [45] gridExtra_2.3                     ArchR_1.0.2                      
## 
## loaded via a namespace (and not attached):
##   [1] dichromat_2.0-0.1           R.methodsS3_1.8.2          
##   [3] progress_1.2.2              nnet_7.3-19                
##   [5] poweRlaw_0.70.6             DT_0.29                    
##   [7] goftest_1.2-3               vctrs_0.6.3                
##   [9] spatstat.random_3.1-6       corpcor_1.6.10             
##  [11] shape_1.4.6                 digest_0.6.33              
##  [13] png_0.1-8                   proxy_0.4-27               
##  [15] pcaPP_2.0-3                 ggrepel_0.9.3              
##  [17] deldir_1.0-9                parallelly_1.36.0          
##  [19] hdrcde_3.4                  MASS_7.3-60                
##  [21] pkgdown_2.0.7               reshape2_1.4.4             
##  [23] foreach_1.5.2               httpuv_1.6.11              
##  [25] withr_2.5.0                 xfun_0.40                  
##  [27] ellipsis_0.3.2              survival_3.5-7             
##  [29] memoise_2.0.1               hexbin_1.28.3              
##  [31] systemfonts_1.0.4           GlobalOptions_0.1.2        
##  [33] ragg_1.2.5                  zoo_1.8-12                 
##  [35] pbapply_1.7-2               R.oo_1.25.0                
##  [37] DEoptimR_1.1-2              Formula_1.2-5              
##  [39] prettyunits_1.1.1           KEGGREST_1.40.0            
##  [41] promises_1.2.1              scatterplot3d_0.3-44       
##  [43] httr_1.4.7                  restfulr_0.0.15            
##  [45] globals_0.16.2              fitdistrplus_1.1-11        
##  [47] rhdf5filters_1.12.1         rstudioapi_0.15.0          
##  [49] miniUI_0.1.1.1              generics_0.1.3             
##  [51] base64enc_0.1-3             curl_5.0.2                 
##  [53] zlibbioc_1.46.0             polyclip_1.10-4            
##  [55] chromVAR_1.22.1             GenomeInfoDbData_1.2.10    
##  [57] RcppEigen_0.3.3.9.3         xtable_1.8-4               
##  [59] desc_1.4.2                  pracma_2.4.2               
##  [61] doParallel_1.0.17           evaluate_0.21              
##  [63] S4Arrays_1.0.6              BiocFileCache_2.8.0        
##  [65] hms_1.1.3                   irlba_2.3.5.1              
##  [67] colorspace_2.1-0            filelock_1.0.2             
##  [69] hdf5r_1.3.8                 ROCR_1.0-11                
##  [71] reticulate_1.32.0           spatstat.data_3.0-1        
##  [73] lmtest_0.9-40               readr_2.1.4                
##  [75] later_1.3.1                 viridis_0.6.4              
##  [77] lattice_0.21-8              spatstat.geom_3.2-5        
##  [79] future.apply_1.11.0         robustbase_0.99-0          
##  [81] scattermore_1.2             XML_3.99-0.14              
##  [83] cowplot_1.1.1               RcppAnnoy_0.0.21           
##  [85] xts_0.13.1                  Hmisc_5.1-1                
##  [87] class_7.3-22                pillar_1.9.0               
##  [89] nlme_3.1-163                iterators_1.0.14           
##  [91] caTools_1.18.2              compiler_4.3.1             
##  [93] RSpectra_0.16-1             stringi_1.7.12             
##  [95] tensor_1.5                  GenomicAlignments_1.36.0   
##  [97] fda_6.1.4                   crayon_1.5.2               
##  [99] abind_1.4-5                 BiocIO_1.10.0              
## [101] graphlayouts_1.0.1          bit_4.0.5                  
## [103] pcaMethods_1.92.0           fastmatch_1.1-4            
## [105] codetools_0.2-19            textshaping_0.3.6          
## [107] TTR_0.24.3                  bslib_0.5.1                
## [109] textshape_1.7.3             e1071_1.7-13               
## [111] destiny_3.14.0              fds_1.8                    
## [113] biovizBase_1.48.0           GetoptLong_1.0.5           
## [115] ggplot.multistats_1.0.0     plotly_4.10.2              
## [117] mime_0.12                   splines_4.3.1              
## [119] ramify_0.3.3                fastDummies_1.7.3          
## [121] dbplyr_2.3.3                interp_1.1-4               
## [123] knitr_1.44                  blob_1.2.4                 
## [125] utf8_1.2.3                  here_1.0.1                 
## [127] clue_0.3-64                 seqLogo_1.66.0             
## [129] fs_1.6.3                    checkmate_2.2.0            
## [131] listenv_0.9.0               Gviz_1.44.1                
## [133] tibble_3.2.1                tzdb_0.4.0                 
## [135] tweenr_2.0.2                pkgconfig_2.0.3            
## [137] tools_4.3.1                 cachem_1.0.8               
## [139] RSQLite_2.3.1               viridisLite_0.4.2          
## [141] smoother_1.1                DBI_1.1.3                  
## [143] fastmap_1.1.1               rmarkdown_2.25             
## [145] scales_1.2.1                ica_1.0-3                  
## [147] Rsamtools_2.16.0            sass_0.4.7                 
## [149] VariantAnnotation_1.46.0    dotCall64_1.0-2            
## [151] carData_3.0-5               rpart_4.1.19               
## [153] RANN_2.6.1                  farver_2.1.1               
## [155] mgcv_1.9-0                  tidygraph_1.2.3            
## [157] yaml_2.3.7                  deSolve_1.38               
## [159] latticeExtra_0.6-30         foreign_0.8-85             
## [161] ggthemes_4.2.4              cli_3.6.1                  
## [163] purrr_1.0.2                 motifmatchr_1.22.0         
## [165] leiden_0.4.3                lifecycle_1.0.3            
## [167] uwot_0.1.16                 rainbow_3.7                
## [169] mvtnorm_1.2-3               backports_1.4.1            
## [171] BiocParallel_1.34.2         annotate_1.78.0            
## [173] rjson_0.2.21                ggridges_0.5.4             
## [175] progressr_0.14.0            parallel_4.3.1             
## [177] jsonlite_1.8.7              RcppHNSW_0.4.1             
## [179] bitops_1.0-7                assertthat_0.2.1           
## [181] bit64_4.0.5                 Rtsne_0.16                 
## [183] nabor_0.5.0                 spatstat.utils_3.0-3       
## [185] ranger_0.15.1               CNEr_1.36.0                
## [187] jquerylib_0.1.4             R.utils_2.12.2             
## [189] lazyeval_0.2.2              shiny_1.7.5                
## [191] htmltools_0.5.6             GO.db_3.17.0               
## [193] sctransform_0.3.5           rappdirs_0.3.3             
## [195] glue_1.6.2                  factoextra_1.0.7           
## [197] TFMPvalue_0.0.9             spam_2.9-1                 
## [199] VIM_6.2.2                   RCurl_1.98-1.12            
## [201] rprojroot_2.0.3             mclust_6.0.0               
## [203] ks_1.14.1                   jpeg_0.1-10                
## [205] boot_1.3-28.1               R6_2.5.1                   
## [207] tidyr_1.3.0                 SingleCellExperiment_1.22.0
## [209] labeling_0.4.3              RcppRoll_0.3.0             
## [211] vcd_1.4-11                  cluster_2.1.4              
## [213] Rhdf5lib_1.22.1             DirichletMultinomial_1.42.0
## [215] DelayedArray_0.26.7         tidyselect_1.2.0           
## [217] ProtGenerics_1.32.0         htmlTable_2.4.1            
## [219] ggforce_0.4.1               xml2_1.3.5                 
## [221] car_3.1-2                   future_1.33.0              
## [223] munsell_0.5.0               KernSmooth_2.23-22         
## [225] laeken_0.5.2                htmlwidgets_1.6.2          
## [227] RColorBrewer_1.1-3          biomaRt_2.56.1             
## [229] rlang_1.1.1                 spatstat.sparse_3.0-2      
## [231] spatstat.explore_3.2-3      fansi_1.0.4