Detecting hidden heterogeneity in single cell RNA-Seq data
The iasva package can be used to detect hidden heterogenity within bulk or single cell sequencing data. To illustrate how to use the iasva package for heterogenity detection, we use real-world single cell RNA sequencing (scRNA-Seq) data obtained from human pancreatic islet samples (Lawlor et. al., 2016).
Load packages
Load the islet single cell RNA-Seq data
To illustrate how IA-SVA can be used to detect hidden heterogeneity within a homogenous cell population (i.e., alpha cells), we use read counts of alpha cells from healthy (non-diabetic) subjects (n = 101).
counts_file <- system.file("extdata", "iasva_counts_test.Rds",
package = "iasva")
# matrix of read counts where genes are rows, samples are columns
counts <- readRDS(counts_file)
# matrix of sample annotations/metadata
anns_file <- system.file("extdata", "iasva_anns_test.Rds",
package = "iasva")
anns <- readRDS(anns_file)Calculate geometric library size, i.e., library size of log-transfromed read counts.
It is well known that the geometric library size (i.e., library size of log-transfromed read counts) or proportion of expressed genes in each cell explains a very large portion of variability of scRNA-Seq data (Hicks et. al. 2015 BioRxiv, McDavid et. al. 2016 Nature Biotechnology). Frequently, the first principal component of log-transformed scRNA-Seq read counts is highly correlated with the geometric library size (r ~ 0.9). Here, we calculate the geometric library size vector, which will be used as a known factor in the IA-SVA algorithm.
geo_lib_size <- colSums(log(counts + 1))
barplot(geo_lib_size, xlab = "Cell", ylab = "Geometric Lib Size", las = 2)
lcounts <- log(counts + 1)
# PC1 and Geometric library size correlation
pc1 <- irlba(lcounts - rowMeans(lcounts), 1)$v[, 1]
cor(geo_lib_size, pc1)## [1] -0.99716Run IA-SVA
Here, we run IA-SVA using patient_id and geo_lib_size as known factors and identify five hidden factors. SVs are plotted in a pairwise fashion to uncover which SVs can seperate cell types.
set.seed(100)
patient_id <- anns$Patient_ID
mod <- model.matrix(~patient_id + geo_lib_size)
# create a summarizedexperiment class
summ_exp <- SummarizedExperiment(assays = counts)
iasva.res<- iasva(summ_exp, mod[, -1],verbose = FALSE,
permute = FALSE, num.sv = 5)## IA-SVA running...##
## SV 1 Detected!##
## SV 2 Detected!##
## SV 3 Detected!##
## SV 4 Detected!##
## SV 5 Detected!##
## # of significant surrogate variables: 5
cell_type <- as.factor(iasva.sv[, 1] > -0.1)
levels(cell_type) <- c("Cell1", "Cell2")
table(cell_type)## cell_type
## Cell1 Cell2
## 6 95# We identified 6 outlier cells based on SV1 that are marked in red
pairs(iasva.sv, main = "IA-SVA", pch = 21, col = color.vec[cell_type],
bg = color.vec[cell_type], oma = c(4,4,6,12))
legend("right", levels(cell_type), fill = color.vec, bty = "n")
plot(iasva.sv[, 1:2], main = "IA-SVA", pch = 21, xlab = "SV1", ylab = "SV2",
col = color.vec[cell_type], bg = color.vec[cell_type])
## [1] -0.1469422
As shown in the above figure, SV1 clearly separates alpha cells into two groups: 6 outlier cells (marked in red) and the rest of the alpha cells (marked in blue). ## Find marker genes for the detected heterogeneity (SV1). Here, using the find_markers() function we find marker genes that are significantly associated with SV1 (multiple testing adjusted p-value < 0.05, default significance cutoff, and R-squared value > 0.3, default R-squared cutoff).
## # of markers (): 33## total # of unique markers: 33## [1] 33## [1] "PMEPA1" "FAM198B" "FLT1" "ENG" "SOX4" "ITGA5"
## [7] "PXDN" "PRDM1" "ERG" "CLIC4" "A2M" "PPAP2B"
## [13] "THBS1" "CLIC2" "S100A16" "STC1" "ACVRL1" "COL4A1"
## [19] "MSN" "TNFAIP2" "MMP2" "SERPINE1" "SPARC" "SPARCL1"
## [25] "ESAM" "KDR" "CD9" "CXCR4" "PODXL" "PLVAP"
## [31] "CALD1" "MMP1" "ADAMTS4"anno.col <- data.frame(cell_type = cell_type)
rownames(anno.col) <- colnames(marker.counts)
head(anno.col)## cell_type
## 4th-C63_S30 Cell2
## 4th-C66_S36 Cell2
## 4th-C18_S31 Cell2
## 4th-C57_S18 Cell1
## 4th-C56_S17 Cell2
## 4th-C68_S41 Cell2pheatmap(log(marker.counts + 1), show_colnames = FALSE,
clustering_method = "ward.D2", cutree_cols = 2,
annotation_col = anno.col)
Run tSNE post IA-SVA analyses, i.e., run tSNE on marker genes associated with SV1 as detected by IA-SVA.
Here, we apply tSNE on the marker genes for SV1 obtained from IA-SVA
set.seed(100)
tsne.res <- Rtsne(unique(t(log(marker.counts + 1))), dims = 2)
plot(tsne.res$Y, main = "tSNE post IA-SVA", xlab = "tSNE Dim1",
ylab = "tSNE Dim2", pch = 21, col = color.vec[cell_type],
bg = color.vec[cell_type], oma = c(4, 4, 6, 12))
legend("bottomright", levels(cell_type), fill = color.vec, bty = "n")
tSNE using SV1 marker genes better seperate these ourlier cells. This analyses suggest that gene selection using IA-SVA combined with tSNE analyses can be a powerful way to detect rare cells introducing variability in the single cell gene expression data.
Using a faster implementation of IA-SVA (fast_iasva)
Here, we run a faster implementation of IA-SVA using the same known factors (patient_id and geo_lib_size) as demonstrated above. This function is useful when working with particularly large datasets.
## fast IA-SVA running...##
## SV 1 Detected!##
## SV 2 Detected!##
## SV 3 Detected!##
## SV 4 Detected!##
## SV 5 Detected!##
## # of obtained surrogate variables: 5Tuning parameters for IA-SVA
The R-squared thresholds used to identify marker genes (find_markers) can greatly influence the 1) number of marker genes identified and 2) the quality of clustering results. With the study_R2() function, users can visualize how different R-squared thresholds influence both factors.
## # of markers (): 274## total # of unique markers: 274## # of markers (): 177## total # of unique markers: 177## # of markers (): 123## total # of unique markers: 123## # of markers (): 84## total # of unique markers: 84## # of markers (): 54## total # of unique markers: 54## # of markers (): 39## total # of unique markers: 39## # of markers (): 30## total # of unique markers: 30## # of markers (): 23## total # of unique markers: 23## # of markers (): 7## total # of unique markers: 7## # of markers (): 3## total # of unique markers: 3## # of markers (): 0## total # of unique markers: 0
This function produces a plot of the number of genes selected vs. the
cluster quality (average silhouette score) for different R-squared
values.
Session Info
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.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=C
## [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 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SummarizedExperiment_1.35.5 Biobase_2.67.0
## [3] GenomicRanges_1.57.2 GenomeInfoDb_1.41.2
## [5] IRanges_2.39.2 S4Vectors_0.43.2
## [7] BiocGenerics_0.53.0 MatrixGenerics_1.17.1
## [9] matrixStats_1.4.1 RColorBrewer_1.1-3
## [11] DescTools_0.99.57 corrplot_0.95
## [13] pheatmap_1.0.12 Rtsne_0.17
## [15] sva_3.53.0 BiocParallel_1.41.0
## [17] genefilter_1.87.0 mgcv_1.9-1
## [19] nlme_3.1-166 iasva_1.25.0
## [21] irlba_2.3.5.1 Matrix_1.7-1
## [23] rmarkdown_2.28
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 gld_2.6.6 readxl_1.4.3
## [4] rlang_1.1.4 e1071_1.7-16 compiler_4.4.1
## [7] RSQLite_2.3.7 png_0.1-8 vctrs_0.6.5
## [10] crayon_1.5.3 fastmap_1.2.0 XVector_0.45.0
## [13] UCSC.utils_1.1.0 bit_4.5.0 xfun_0.48
## [16] zlibbioc_1.51.2 cachem_1.1.0 jsonlite_1.8.9
## [19] blob_1.2.4 highr_0.11 DelayedArray_0.33.1
## [22] parallel_4.4.1 cluster_2.1.6 R6_2.5.1
## [25] bslib_0.8.0 limma_3.61.12 boot_1.3-31
## [28] jquerylib_0.1.4 cellranger_1.1.0 Rcpp_1.0.13
## [31] knitr_1.48 splines_4.4.1 rstudioapi_0.17.1
## [34] abind_1.4-8 yaml_2.3.10 codetools_0.2-20
## [37] lattice_0.22-6 withr_3.0.2 KEGGREST_1.45.1
## [40] evaluate_1.0.1 survival_3.7-0 proxy_0.4-27
## [43] Biostrings_2.75.0 munsell_0.5.1 scales_1.3.0
## [46] rootSolve_1.8.2.4 xtable_1.8-4 class_7.3-22
## [49] glue_1.8.0 lmom_3.2 maketools_1.3.1
## [52] tools_4.4.1 sys_3.4.3 data.table_1.16.2
## [55] annotate_1.85.0 locfit_1.5-9.10 Exact_3.3
## [58] buildtools_1.0.0 mvtnorm_1.3-1 XML_3.99-0.17
## [61] grid_4.4.1 AnnotationDbi_1.69.0 edgeR_4.3.21
## [64] colorspace_2.1-1 GenomeInfoDbData_1.2.13 cli_3.6.3
## [67] expm_1.0-0 S4Arrays_1.5.11 gtable_0.3.6
## [70] sass_0.4.9 digest_0.6.37 SparseArray_1.5.45
## [73] farver_2.1.2 memoise_2.0.1 htmltools_0.5.8.1
## [76] lifecycle_1.0.4 httr_1.4.7 statmod_1.5.0
## [79] bit64_4.5.2 MASS_7.3-61- Load packages
- Load the islet single cell RNA-Seq data
- Calculate geometric library size, i.e., library size of log-transfromed read counts.
- Run IA-SVA
- Run tSNE post IA-SVA analyses, i.e., run tSNE on marker genes associated with SV1 as detected by IA-SVA.
- Using a faster implementation of IA-SVA (fast_iasva)
- Tuning parameters for IA-SVA
- Session Info