pbmc_10x_multiome.Rmd
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:
suppressMessages(library(ArchR))
suppressMessages(library(Seurat))
suppressMessages(library(Signac))
suppressMessages(library(scMEGA))
suppressMessages(library(Nebulosa))
suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
suppressMessages(library(GenomeInfoDb))
suppressMessages(library(EnsDb.Hsapiens.v86))
suppressMessages(library(ggplot2))
suppressMessages(library(dplyr))
suppressMessages(library(JASPAR2020))
suppressMessages(library(TFBSTools))
suppressMessages(library(igraph))
suppressMessages(library(ggraph))
suppressMessages(library(MOJITOO))
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)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.
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)
## 1 layer present: counts
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)
## 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
obj.atac <- obj.atac %>%
RunTFIDF() %>%
FindTopFeatures() %>%
RunSVD() %>%
RunUMAP(reduction = 'lsi', dims = 2:30, verbose = FALSE)
## Performing TF-IDF normalization
## Running SVD
## Scaling cell embeddings
Visualize the scRNA-seq and scATAC-seq separately.
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)
## Centering and scaling data matrix
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:
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 = GetAssayData(obj.rna, layer = "counts"),
assay = "RNA",
meta.data = meta.data
)
# create ATAC assay and add it to the object
pbmc[["ATAC"]] <- CreateChromatinAssay(
counts = GetAssayData(obj.atac, layer = "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)
## 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)
DimPlot(pbmc, group.by = "predicted.id",
shuffle = TRUE, label = TRUE, reduction = "MOJITOO_UMAP") + NoLegend()
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
We next select candidate TFs and genes for building a meaningful gene regulatory network.
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
## Warning: Layer counts isn't present in the assay object; returning NULL
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)
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)
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.
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)
## Warning: Layer counts isn't present in the assay object; returning NULL
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.3 (2024-02-29)
## 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_envs/zl_scmega/lib/libopenblasp-r0.3.27.so; LAPACK version 3.12.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: Etc/UTC
## 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.16 ComplexHeatmap_2.18.0
## [3] MOJITOO_1.0 ggraph_2.2.1
## [5] igraph_2.0.3 TFBSTools_1.40.0
## [7] JASPAR2020_0.99.10 dplyr_1.1.4
## [9] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.26.0
## [11] AnnotationFilter_1.26.0 GenomicFeatures_1.54.1
## [13] AnnotationDbi_1.64.1 BSgenome.Hsapiens.UCSC.hg38_1.4.5
## [15] BSgenome_1.70.2 rtracklayer_1.62.0
## [17] BiocIO_1.12.0 Biostrings_2.70.3
## [19] XVector_0.42.0 Nebulosa_1.12.1
## [21] patchwork_1.2.0 scMEGA_1.0.2
## [23] Signac_1.13.0 Seurat_5.1.0
## [25] SeuratObject_5.0.2 sp_2.1-4
## [27] rhdf5_2.46.1 SummarizedExperiment_1.32.0
## [29] Biobase_2.62.0 MatrixGenerics_1.14.0
## [31] Rcpp_1.0.12 Matrix_1.6-5
## [33] GenomicRanges_1.54.1 GenomeInfoDb_1.38.8
## [35] IRanges_2.36.0 S4Vectors_0.40.2
## [37] BiocGenerics_0.48.1 matrixStats_1.3.0
## [39] data.table_1.15.4 stringr_1.5.1
## [41] plyr_1.8.9 magrittr_2.0.3
## [43] ggplot2_3.5.1 gtable_0.3.5
## [45] gtools_3.9.5 gridExtra_2.3
## [47] ArchR_1.0.2
##
## loaded via a namespace (and not attached):
## [1] ica_1.0-3 plotly_4.10.4
## [3] Formula_1.2-5 zlibbioc_1.48.2
## [5] tidyselect_1.2.1 bit_4.0.5
## [7] doParallel_1.0.17 clue_0.3-65
## [9] lattice_0.22-6 rjson_0.2.21
## [11] factoextra_1.0.7 blob_1.2.4
## [13] S4Arrays_1.2.1 parallel_4.3.3
## [15] dichromat_2.0-0.1 hdrcde_3.4
## [17] seqLogo_1.68.0 png_0.1-8
## [19] cli_3.6.3 ProtGenerics_1.34.0
## [21] goftest_1.2-3 VIM_6.2.2
## [23] textshape_1.7.5 pkgdown_2.0.9
## [25] textshaping_0.3.7 purrr_1.0.2
## [27] chromVAR_1.24.0 uwot_0.2.2
## [29] curl_5.2.1 deSolve_1.40
## [31] mime_0.12 evaluate_0.24.0
## [33] leiden_0.4.3.1 stringi_1.8.4
## [35] backports_1.5.0 desc_1.4.3
## [37] XML_3.99-0.17 httpuv_1.6.15
## [39] rappdirs_0.3.3 splines_4.3.3
## [41] RcppRoll_0.3.0 mclust_6.1.1
## [43] jpeg_0.1-10 rainbow_3.8
## [45] pcaPP_2.0-4 DT_0.33
## [47] sctransform_0.4.1 ggbeeswarm_0.7.2
## [49] DBI_1.2.3 jquerylib_0.1.4
## [51] smoother_1.3 withr_3.0.0
## [53] corpcor_1.6.10 class_7.3-22
## [55] systemfonts_1.0.5 rprojroot_2.0.4
## [57] lmtest_0.9-40 tidygraph_1.3.1
## [59] htmlwidgets_1.6.4 fs_1.6.4
## [61] biomaRt_2.58.0 SingleCellExperiment_1.24.0
## [63] ggrepel_0.9.5 labeling_0.4.3
## [65] SparseArray_1.2.4 ranger_0.16.0
## [67] DEoptimR_1.1-3 annotate_1.80.0
## [69] reticulate_1.38.0 VariantAnnotation_1.48.1
## [71] hexbin_1.28.3 zoo_1.8-12
## [73] knitr_1.47 ggplot.multistats_1.0.0
## [75] TFMPvalue_0.0.9 foreach_1.5.2
## [77] fansi_1.0.6 fda_6.1.8
## [79] caTools_1.18.2 R.oo_1.26.0
## [81] poweRlaw_0.80.0 RSpectra_0.16-1
## [83] irlba_2.3.5.1 ggrastr_1.0.2
## [85] fastDummies_1.7.3 lazyeval_0.2.2
## [87] yaml_2.3.9 survival_3.7-0
## [89] scattermore_1.2 crayon_1.5.3
## [91] RcppAnnoy_0.0.22 RColorBrewer_1.1-3
## [93] tidyr_1.3.1 progressr_0.14.0
## [95] tweenr_2.0.3 later_1.3.2
## [97] ggridges_0.5.6 fds_1.8
## [99] codetools_0.2-20 base64enc_0.1-3
## [101] GlobalOptions_0.1.2 KEGGREST_1.42.0
## [103] Rtsne_0.17 shape_1.4.6.1
## [105] Rsamtools_2.18.0 filelock_1.0.3
## [107] foreign_0.8-87 pkgconfig_2.0.3
## [109] xml2_1.3.6 GenomicAlignments_1.38.2
## [111] scatterplot3d_0.3-44 spatstat.sparse_3.1-0
## [113] viridisLite_0.4.2 biovizBase_1.50.0
## [115] xtable_1.8-4 interp_1.1-6
## [117] car_3.1-2 highr_0.11
## [119] httr_1.4.7 tools_4.3.3
## [121] globals_0.16.3 beeswarm_0.4.0
## [123] htmlTable_2.4.2 checkmate_2.3.1
## [125] nlme_3.1-165 dbplyr_2.5.0
## [127] hdf5r_1.3.10 assertthat_0.2.1
## [129] digest_0.6.36 farver_2.1.2
## [131] tzdb_0.4.0 reshape2_1.4.4
## [133] ks_1.14.2 viridis_0.6.5
## [135] DirichletMultinomial_1.44.0 rpart_4.1.23
## [137] glue_1.7.0 cachem_1.1.0
## [139] BiocFileCache_2.10.1 polyclip_1.10-6
## [141] ramify_0.3.3 Hmisc_5.1-3
## [143] generics_0.1.3 mvtnorm_1.2-5
## [145] parallelly_1.37.1 here_1.0.1
## [147] RcppHNSW_0.6.0 ragg_1.3.2
## [149] carData_3.0-5 pbapply_1.7-2
## [151] nabor_0.5.0 spam_2.10-0
## [153] utf8_1.2.4 graphlayouts_1.1.1
## [155] RcppEigen_0.3.4.0.0 shiny_1.8.1.1
## [157] GenomeInfoDbData_1.2.11 R.utils_2.12.3
## [159] rhdf5filters_1.14.1 RCurl_1.98-1.14
## [161] memoise_2.0.1 rmarkdown_2.27
## [163] scales_1.3.0 R.methodsS3_1.8.2
## [165] future_1.33.2 RANN_2.6.1
## [167] Cairo_1.6-2 spatstat.data_3.1-2
## [169] rstudioapi_0.16.0 cluster_2.1.6
## [171] spatstat.utils_3.0-5 hms_1.1.3
## [173] fitdistrplus_1.1-11 munsell_0.5.1
## [175] cowplot_1.1.3 colorspace_2.1-0
## [177] rlang_1.1.4 xts_0.14.0
## [179] dotCall64_1.1-1 ggforce_0.4.2
## [181] mgcv_1.9-1 laeken_0.5.3
## [183] xfun_0.45 e1071_1.7-14
## [185] CNEr_1.38.0 iterators_1.0.14
## [187] abind_1.4-5 tibble_3.2.1
## [189] Rhdf5lib_1.24.2 motifmatchr_1.24.0
## [191] readr_2.1.5 bitops_1.0-7
## [193] promises_1.3.0 RSQLite_2.3.7
## [195] DelayedArray_0.28.0 proxy_0.4-27
## [197] GO.db_3.18.0 compiler_4.3.3
## [199] prettyunits_1.2.0 boot_1.3-30
## [201] listenv_0.9.1 tensor_1.5
## [203] MASS_7.3-60 progress_1.2.3
## [205] BiocParallel_1.36.0 spatstat.random_3.2-3
## [207] R6_2.5.1 fastmap_1.2.0
## [209] fastmatch_1.1-4 vipor_0.4.7
## [211] TTR_0.24.4 ROCR_1.0-11
## [213] vcd_1.4-12 nnet_7.3-19
## [215] KernSmooth_2.23-24 latticeExtra_0.6-30
## [217] miniUI_0.1.1.1 deldir_2.0-4
## [219] htmltools_0.5.8.1 ggthemes_5.1.0
## [221] bit64_4.0.5 spatstat.explore_3.2-7
## [223] lifecycle_1.0.4 destiny_3.16.0
## [225] Gviz_1.46.1 restfulr_0.0.15
## [227] sass_0.4.9 vctrs_0.6.5
## [229] robustbase_0.99-2 spatstat.geom_3.2-9
## [231] future.apply_1.11.2 pracma_2.4.4
## [233] bslib_0.7.0 pillar_1.9.0
## [235] pcaMethods_1.94.0 jsonlite_1.8.8
## [237] GetoptLong_1.0.5