Differential gene expression
In this tutorial we will cover about Differetial gene expression, which comprises an extensive range of topics and methods. In single cell, differential expresison can have multiple functionalities such as of identifying marker genes for cell populations, as well as differentially regulated genes across conditions (healthy vs control). We will also exercise on how to account the batch information in your test.
We can first load the data from the clustering session. Moreover, we can already decide which clustering resolution to use. First let’s define using the louvain_2 clustering to identifying differentially expressed genes.
suppressPackageStartupMessages({
library(scater)
library(scran)
# library(venn)
library(cowplot)
library(ggplot2)
library(rafalib)
library(pheatmap)
library(igraph)
library(dplyr)
})
sce <- readRDS("data/3pbmc_qc_dr_int_cl.rds")Cell marker genes
Let us first compute a ranking for the highly differential genes in each cluster. There are many different tests and parameters to be chosen that can be used to refine your results. When looking for marker genes, we want genes that are positivelly expressed in a cell type and possibly not expressed in the others.
# Compute differentiall expression
markers_genes <- findMarkers(x = sce, clusters = sce$kmeans_5, lfc = 0.5, pval.type = "all", direction = "up")
# List of dataFrames with the results for each cluster
markers_genes## DataFrameList of length 5
## names(5): 1 2 3 4 5## DataFrame with 16157 rows and 6 columns
## p.value FDR logFC.2 logFC.3 logFC.4 logFC.5
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## CDKN1C 7.35183149155511e-06 0.118783541409056 0.979144387210788 0.998949056933944 0.969011598924564 0.991495303801578
## TCF7L2 0.00012879246392497 1 0.815057815932671 1.03779310507903 1.02556879447223 1.00178886845701
## SMIM25 0.00124542854863882 1 0.839447625697549 1.5731720846299 1.567609623075 1.5731720846299
## CSF1R 0.00975722085979975 1 0.718391145224288 1.28866050964937 1.30643924791123 1.3020534882496
## HES4 0.0159258780596399 1 0.66438917326519 0.726808230474354 0.730094609542395 0.729703856857455
## ... ... ... ... ... ... ...
## COL6A2 1 1 -0.000772237217109488 -0.00358712990462427 -0.199559233102298 -0.470972301125049
## C21orf58 1 1 -0.00365849065946177 -0.0587826143860415 -0.0291564527523501 0.00234371542516463
## BX004987.1 1 1 0.00232376715326658 0.00232376715326658 -0.000699738532482444 0.00232376715326658
## AC011043.1 1 1 0.000793415321419471 -0.00240030584993707 8.78541764944237e-05 -0.0202713946352216
## AL592183.1 1 1 -0.00166339863031148 -0.0122041610684985 -0.00490734614857844 -0.00423474851640378We can now select the top 25 up regulated genes for plotting.
# Colect the top 25 genes for each cluster and put the into a single table
top25 <- lapply(names(markers_genes), function(x) {
temp <- markers_genes[[x]][1:25, 1:2]
temp$gene <- rownames(markers_genes[[x]])[1:25]
temp$cluster <- x
return(temp)
})
top25 <- as_tibble(do.call(rbind, top25))
top25We can now select the top 25 up regulated genes for plotting.
mypar(1, 5, mar = c(4, 6, 3, 1))
for (i in unique(top25$cluster)) {
barplot(sort(setNames(-log10(top25$p.value), top25$gene)[top25$cluster == i], F), horiz = T, las = 1, main = paste0(i, " vs. rest"), border = "white", yaxs = "i", xlab = "-log10FC")
abline(v = c(0, -log10(0.05)), lty = c(1, 2))
}
We can visualize them as a heatmap. Here we are selecting the top 5.
top5 <- as_tibble(top25) %>% group_by(cluster) %>% top_n(-5, p.value)
scater::plotHeatmap(sce[, order(sce$kmeans_5)], features = unique(top5$gene), center = T, zlim = c(-3, 3), colour_columns_by = "kmeans_5", show_colnames = F, cluster_cols = F, fontsize_row = 6,
color = colorRampPalette(c("purple", "black", "yellow"))(90))
We can also plot a violin plot for each gene.
Differential expression across conditions
The second way of computing differential expression is to answer which genes are differentially expressed within a cluster. For example, in our case we have libraries comming from 2 different library preparation methods (batches) and we would like to know which genes are influenced the most in a particular cell type. The same concenpt applies if you have instead two or more biological groups (control vs treated, time#0 vs time#1 vs time#2, etc).
For this end, we will first subset our data for the desired cell cluster, then change the cell identities to the variable of comparison (which now in our case is the “Chemistry”).
# Filter cells from that cluster
cell_selection <- sce[, sce$kmeans_5 == 4]
cell_selection$Chemistry <- ifelse(cell_selection$sample_id == "v2.1k", "v2", "v3")
# Compute differentiall expression
DGE_cell_selection <- findMarkers(x = cell_selection, clusters = cell_selection$Chemistry, lfc = 0.5, pval.type = "all", direction = "down")
top5_cell_selection <- lapply(names(DGE_cell_selection), function(x) {
temp <- DGE_cell_selection[[x]][1:5, 1:2]
temp$gene <- rownames(DGE_cell_selection[[x]])[1:5]
temp$cluster <- x
return(temp)
})
top5_cell_selection <- as_tibble(do.call(rbind, top5_cell_selection))
top5_cell_selectionWe can now plot the expression across the “Chemistry”.
We can clearly see some patterns across them. Those are the genes that impact the most on your batches (see the dimensionality reduction and integration exercises for more details). We can plot those genes using the integrated and non-integrated UMAP for ilustration.
plotlist <- list()
for (i in c("JUND", "RPS17", "GNAS")) {
plotlist[[i]] <- plotReducedDim(sce, use_dimred = "UMAP_on_PCA", colour_by = i, by_exprs_values = "logcounts", add_ticks = F) + scale_fill_gradientn(colours = colorRampPalette(c("grey90",
"orange3", "firebrick", "firebrick", "red", "red"))(10)) + ggtitle(label = i) + theme(plot.title = element_text(size = 20))
}
plot_grid(ncol = 3, plotlist = plotlist)
Finally, lets save the integrated data for further analysis.
Session Info
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] venn_1.7 Seurat_3.0.2 dplyr_0.8.0.1 igraph_1.2.4.1 pheatmap_1.0.12 rafalib_1.0.0
## [7] cowplot_0.9.4 scran_1.10.2 scater_1.10.1 ggplot2_3.1.1 SingleCellExperiment_1.4.1 SummarizedExperiment_1.12.0
## [13] DelayedArray_0.8.0 BiocParallel_1.16.6 matrixStats_0.54.0 Biobase_2.42.0 GenomicRanges_1.34.0 GenomeInfoDb_1.18.2
## [19] IRanges_2.16.0 S4Vectors_0.20.1 BiocGenerics_0.28.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 ggbeeswarm_0.6.0 colorspace_1.4-1 ggridges_0.5.1 dynamicTreeCut_1.63-1 XVector_0.22.0 BiocNeighbors_1.0.0
## [8] rstudioapi_0.10 listenv_0.7.0 npsurv_0.4-0 ggrepel_0.8.0 codetools_0.2-16 splines_3.5.1 R.methodsS3_1.7.1
## [15] lsei_1.2-0 knitr_1.26 jsonlite_1.6 ica_1.0-2 cluster_2.0.9 png_0.1-7 R.oo_1.22.0
## [22] sctransform_0.2.0 HDF5Array_1.10.1 httr_1.4.0 compiler_3.5.1 assertthat_0.2.1 Matrix_1.2-17 lazyeval_0.2.2
## [29] limma_3.38.3 formatR_1.6 htmltools_0.4.0 tools_3.5.1 rsvd_1.0.0 gtable_0.3.0 glue_1.3.1
## [36] GenomeInfoDbData_1.2.0 RANN_2.6.1 reshape2_1.4.3 Rcpp_1.0.3 gdata_2.18.0 ape_5.3 nlme_3.1-139
## [43] DelayedMatrixStats_1.4.0 gbRd_0.4-11 lmtest_0.9-37 xfun_0.11 stringr_1.4.0 globals_0.12.4 irlba_2.3.3
## [50] gtools_3.8.1 statmod_1.4.30 future_1.12.0 edgeR_3.24.3 zlibbioc_1.28.0 MASS_7.3-51.4 zoo_1.8-5
## [57] scales_1.0.0 rhdf5_2.26.2 RColorBrewer_1.1-2 yaml_2.2.0 reticulate_1.12 pbapply_1.4-0 gridExtra_2.3
## [64] stringi_1.4.3 caTools_1.17.1.2 bibtex_0.4.2 Rdpack_0.11-0 SDMTools_1.1-221.1 rlang_0.4.1 pkgconfig_2.0.2
## [71] bitops_1.0-6 evaluate_0.14 lattice_0.20-38 ROCR_1.0-7 purrr_0.3.2 Rhdf5lib_1.4.3 htmlwidgets_1.3
## [78] labeling_0.3 tidyselect_0.2.5 plyr_1.8.4 magrittr_1.5 R6_2.4.1 gplots_3.0.1.1 pillar_1.3.1
## [85] withr_2.1.2 fitdistrplus_1.0-14 survival_2.44-1.1 RCurl_1.95-4.12 tsne_0.1-3 tibble_2.1.1 future.apply_1.2.0
## [92] crayon_1.3.4 KernSmooth_2.23-15 plotly_4.9.0 rmarkdown_1.17 viridis_0.5.1 locfit_1.5-9.1 grid_3.5.1
## [99] data.table_1.12.2 metap_1.1 digest_0.6.22 tidyr_0.8.3 R.utils_2.8.0 munsell_0.5.0 beeswarm_0.2.3
## [106] viridisLite_0.3.0 vipor_0.4.5