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 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.99716
Run 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
iasva.sv <- iasva.res$sv
plot(iasva.sv[, 1], iasva.sv[, 2], xlab = "SV1", ylab = "SV2")
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])
cor(geo_lib_size, iasva.sv[, 1])
## [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).
marker.counts <- find_markers(summ_exp, as.matrix(iasva.sv[,1]))
## # 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 Cell2
pheatmap(log(marker.counts + 1), show_colnames = FALSE,
clustering_method = "ward.D2", cutree_cols = 2,
annotation_col = anno.col)
Run tSNE to detect the hidden heterogeneity.
For comparison purposes, we applied tSNE on read counts of all genes
to identify the hidden heterogeneity. We used the Rtsne R package with
default settings.
set.seed(100)
tsne.res <- Rtsne(t(lcounts), dims = 2)
plot(tsne.res$Y, main = "tSNE", 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")
As shown above, tSNE fails to detect the outlier cells that are
identified by IA-SVA when all genes are used. Same color coding is used
as above.
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.
iasva.res <- fast_iasva(summ_exp, mod[, -1], num.sv = 5)
## fast IA-SVA running...
##
## SV 1 Detected!
##
## SV 2 Detected!
##
## SV 3 Detected!
##
## SV 4 Detected!
##
## SV 5 Detected!
##
## # of obtained surrogate variables: 5
Tuning 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.
study_res <- study_R2(summ_exp, iasva.sv)
## # 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.