SS2 (for AG Hefendehl): Integration
Desirée Brösamle
Dec 2021
Last updated: 2022-02-21
Checks: 7 0
Knit directory: S:/KJP_Biolabor/Projects/scSeq_Hefendehl/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20220131) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 73aabec. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rproj.user/
Ignored: data/ReloadAllData_Hefendehl_Stroke_Dec'21.RData
Ignored: data/Sample_Tables/
Ignored: data/counts.csv
Ignored: data/genecounts.csv
Ignored: data/microglia_protein.rds
Ignored: data/samples.integrated.RData
Ignored: data/tx2genes.csv
Ignored: output/Descriptives.Rmd
Ignored: output/Descriptives.docx
Untracked files:
Untracked: app/
Untracked: geneviewer/
Untracked: genplotter/
Untracked: workflow_helper.R
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
There are no past versions. Publish this analysis with wflow_publish() to start tracking its development.
Preprocessing
QC metrics from AG Hefendehl (Stroke)
# Visualize QC metrics as a violin plot
Idents(plates_hefendehl) <- "Plate"
VlnPlot(plates_hefendehl, features = c("nFeature_RNA", "nCount_RNA", "percent.mito","percent.ribo"),
group.by = "Plate", pt.size = 1, ncol = 4, )
plot1 | plot2
QC metrics from AG Neher (own additional WT control)
# plates_neher <- subset(plates_neher, subset = Genotype == "wt")
# Visualize QC metrics as a violin plot
Idents(plates_neher) <- "Plate"
VlnPlot(plates_neher, features = c("nFeature_RNA", "nCount_RNA", "percent.mito","percent.ribo"),
group.by = "Plate",pt.size =1, ncol = 4)
plot1a | plot2a
Combine both seurat objects and filter cells that have unique feature counts (gene number) over 7000 or less than 200 and cells with >5% mitochondrial counts.
QC metrics after filtering
plot3 | plot4
Perform integrated analysis
PCA, tSNE and UMAP from the first 30 dimensions and with a resolution of 0.8
# Visualization
Idents(samples.integrated) <- "seurat_clusters"
DimPlot(samples.integrated, reduction = "umap", split.by = "Treatment")
p1 + p1a + p1b
DimPlot(samples.integrated, group.by = "seurat_clusters", pt.size =1, label = T) + NoLegend()
DimPlot(samples.integrated, reduction = "umap", split.by = "Plate", pt.size = 1)
DimPlot(samples.integrated, reduction = "umap", split.by = "Genotype_Treatment", pt.size =1)
Cell cycle scoring
First, a score is assigned to each cell (Tirosh et al. 2016), based on its expression of G2/M and S phase markers. These markers should be anticorrelated in their expression levels and cells expression neither are likely not cycling and in G1 phase. Note: For downstream cell cylce regression the quantitative scores for G2/M and S phase are used, not the dicrete classification.
Cluster distribution
DimPlot(object = samples.integrated, pt.size = 1,reduction = "umap", group.by="Age",label = F) +
ggtitle("Cluster distribution according to Age") | DimPlot(object = samples.integrated, pt.size = 1,reduction = "umap", group.by="Genotype_Treatment",label = F) +
ggtitle("Cluster distribution according to Genotype + Treatment")
DimPlot(object = samples.integrated, pt.size = 1,reduction = "umap", group.by="Plate",label = F) +
ggtitle("Cluster distribution according to Plate") | DimPlot(object = samples.integrated, pt.size = 1,reduction = "umap", group.by = "seurat_clusters", label = T) +
ggtitle("Clustering")
Idents(samples.integrated) <- "seurat_clusters"
# use same colors for clusters as in plots
require(scales)
identities <- levels(samples.integrated$seurat_clusters) # Create vector with levels of object@ident
cluster_colors <- hue_pal()(length(identities)) # Create vector of default ggplot2 colors
# number of cells in each cluster
cluster_nCell <- as.data.frame.matrix(table(samples.integrated$seurat_clusters,
samples.integrated$Genotype_Age))
cluster_nCell["Total" ,] = colSums(cluster_nCell)
cluster_nCell APPPS1+_10mo APPPS1+_9mo MX04+_10mo MX04+_9mo wt_10mo wt_17mo wt_9mo
0 42 57 20 18 15 20 34
1 23 38 3 7 9 68 16
2 22 29 2 3 2 31 18
3 14 25 12 5 6 29 9
4 4 6 1 0 1 11 3
5 3 4 5 1 5 3 3
6 6 10 0 0 3 0 3
Total 114 169 43 34 41 162 86# % of cells in each cluster , grouped by genotype_Age
cluster_percent_Cell <- data.frame(round((prop.table(x = table(samples.integrated$seurat_clusters,
samples.integrated$Genotype_Age), margin = 2)*100),2))
colnames(cluster_percent_Cell) <- c("Genotype_Age", "Celltype", "Frequency")
cluster_percent_Cell Genotype_Age Celltype Frequency
1 0 APPPS1+_10mo 36.84
2 1 APPPS1+_10mo 20.18
3 2 APPPS1+_10mo 19.30
4 3 APPPS1+_10mo 12.28
5 4 APPPS1+_10mo 3.51
6 5 APPPS1+_10mo 2.63
7 6 APPPS1+_10mo 5.26
8 0 APPPS1+_9mo 33.73
9 1 APPPS1+_9mo 22.49
10 2 APPPS1+_9mo 17.16
11 3 APPPS1+_9mo 14.79
12 4 APPPS1+_9mo 3.55
13 5 APPPS1+_9mo 2.37
14 6 APPPS1+_9mo 5.92
15 0 MX04+_10mo 46.51
16 1 MX04+_10mo 6.98
17 2 MX04+_10mo 4.65
18 3 MX04+_10mo 27.91
19 4 MX04+_10mo 2.33
20 5 MX04+_10mo 11.63
21 6 MX04+_10mo 0.00
22 0 MX04+_9mo 52.94
23 1 MX04+_9mo 20.59
24 2 MX04+_9mo 8.82
25 3 MX04+_9mo 14.71
26 4 MX04+_9mo 0.00
27 5 MX04+_9mo 2.94
28 6 MX04+_9mo 0.00
29 0 wt_10mo 36.59
30 1 wt_10mo 21.95
31 2 wt_10mo 4.88
32 3 wt_10mo 14.63
33 4 wt_10mo 2.44
34 5 wt_10mo 12.20
35 6 wt_10mo 7.32
36 0 wt_17mo 12.35
37 1 wt_17mo 41.98
38 2 wt_17mo 19.14
39 3 wt_17mo 17.90
40 4 wt_17mo 6.79
41 5 wt_17mo 1.85
42 6 wt_17mo 0.00
43 0 wt_9mo 39.53
44 1 wt_9mo 18.60
45 2 wt_9mo 20.93
46 3 wt_9mo 10.47
47 4 wt_9mo 3.49
48 5 wt_9mo 3.49
49 6 wt_9mo 3.49# Grouped barchart of cell proportions
ggplot(cluster_percent_Cell, aes(fill=Genotype_Age, y=Frequency, x=Celltype)) +
geom_bar(position="dodge", stat="identity")+
ggtitle("Cell distribution according to Genotype_Age [%]") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1)) +
xlab("")
# % of cells in each cluster , grouped by genotype_treatment
cluster_percent_Cell <- data.frame(round((prop.table(x = table(samples.integrated$seurat_clusters,
samples.integrated$Genotype_Treatment), margin = 2)*100),2))
colnames(cluster_percent_Cell) <- c("Genotype_Treatment", "Celltype", "Frequency")
cluster_percent_Cell Genotype_Treatment Celltype Frequency
1 0 APPPS1+_Ctrl 30.32
2 1 APPPS1+_Ctrl 23.87
3 2 APPPS1+_Ctrl 21.94
4 3 APPPS1+_Ctrl 14.84
5 4 APPPS1+_Ctrl 5.81
6 5 APPPS1+_Ctrl 1.94
7 6 APPPS1+_Ctrl 1.29
8 0 APPPS1+_Stroke 40.62
9 1 APPPS1+_Stroke 18.75
10 2 APPPS1+_Stroke 13.28
11 3 APPPS1+_Stroke 12.50
12 4 APPPS1+_Stroke 0.78
13 5 APPPS1+_Stroke 3.12
14 6 APPPS1+_Stroke 10.94
15 0 MX04+_Ctrl 50.68
16 1 MX04+_Ctrl 10.96
17 2 MX04+_Ctrl 6.85
18 3 MX04+_Ctrl 21.92
19 4 MX04+_Ctrl 1.37
20 5 MX04+_Ctrl 8.22
21 6 MX04+_Ctrl 0.00
22 0 MX04+_Stroke 25.00
23 1 MX04+_Stroke 50.00
24 2 MX04+_Stroke 0.00
25 3 MX04+_Stroke 25.00
26 4 MX04+_Stroke 0.00
27 5 MX04+_Stroke 0.00
28 6 MX04+_Stroke 0.00
29 0 wt_Ctrl 12.35
30 1 wt_Ctrl 41.98
31 2 wt_Ctrl 19.14
32 3 wt_Ctrl 17.90
33 4 wt_Ctrl 6.79
34 5 wt_Ctrl 1.85
35 6 wt_Ctrl 0.00
36 0 wt_Stroke 38.58
37 1 wt_Stroke 19.69
38 2 wt_Stroke 15.75
39 3 wt_Stroke 11.81
40 4 wt_Stroke 3.15
41 5 wt_Stroke 6.30
42 6 wt_Stroke 4.72# Grouped barchart of cell proportions
ggplot(cluster_percent_Cell, aes(fill=Genotype_Treatment, y=Frequency, x=Celltype)) +
geom_bar(position="dodge", stat="identity")+
ggtitle("Cell distribution according to Genotype_Treatment [%]") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1)) +
xlab("")
Cell type identification
Finding differentially expressed features (cluster biomarkers)
Idents(samples.integrated) <- samples.integrated$seurat_clusters
DoHeatmap(samples.integrated, features = top20$gene) + NoLegend() + ggtitle("Top20 cluster marker genes")
SingleR
SingleR is an automatic annotation method for (scRNAseq) data (Aran et al. 2019). Given a reference dataset of samples (single-cell or bulk) with known labels, it labels new cells from a test dataset based on similarity to the reference set. Here we use the built-in references “Immgen” (830 microarray samples of sorted hematopoetic and immune cell populations) and “Mouse RNA-Seq” (358 non-specific mouse RNA-seq samples).
Immgen reference
Section Skipped as pred.immgen object is not available
MouseRNA-Seq reference Section Skipped as pred.mouseRNA object is not available
# number of cells in each cluster
cluster_nCell <- data.frame(table(samples.integrated$MouseRNASeq_sc_labels,samples.integrated$Genotype_Treatment))
colnames(cluster_nCell) <- c("Genotype_Treatment", "MouseRNASeq_sc_labels", "Number")
#cluster_nCell["Total" ,] = colSums(cluster_nCell)
# cluster_nCell
# % of cells in each cluster
cluster_percent_Cell <- data.frame(round((prop.table(x = table(samples.integrated$MouseRNASeq_sc_labels,
samples.integrated$Genotype_Treatment), margin = 2)*100),2))
colnames(cluster_percent_Cell) <- c("Genotype_Treatment", "MouseRNASeq_sc_labels", "Frequency")
# cluster_percent_Cell
# Grouped barchart of absolute cell numbers
ggplot(cluster_nCell, aes(fill=Genotype_Treatment, y=Number, x=MouseRNASeq_sc_labels)) +
geom_bar(position="dodge", stat="identity") +
ggtitle("Cell distribution according to MouseRNASeq reference (absolut values)") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=0.5)) +
xlab("")
# Grouped barchart of cell proportions
ggplot(cluster_percent_Cell, aes(fill=Genotype_Treatment, y=Frequency, x=MouseRNASeq_sc_labels)) +
geom_bar(position="dodge", stat="identity")+
ggtitle("Cell distribution according to MouseRNASeq reference [%]") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=0.5)) +
xlab("")
# number of cells in each cluster
cluster_nCell <- data.frame(table(samples.integrated$Immgen_sc_labels,samples.integrated$Treatment))
colnames(cluster_nCell) <- c("Treatment", "Immgen_sc_labels", "Number")
#cluster_nCell["Total" ,] = colSums(cluster_nCell)
# cluster_nCell
# % of cells in each cluster
cluster_percent_Cell <- data.frame(round((prop.table(x = table(samples.integrated$Immgen_sc_labels,
samples.integrated$Treatment), margin = 2)*100),2))
colnames(cluster_percent_Cell) <- c("Treatment", "Immgen_sc_labels", "Frequency")
# cluster_percent_Cell
# Grouped barchart of absolute cell numbers
ggplot(cluster_nCell, aes(fill=Treatment, y=Number, x=Immgen_sc_labels)) +
geom_bar(position="dodge", stat="identity") +
ggtitle("Cell distribution according to Immgen reference (absolut values)") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=0.5)) +
xlab("")
# Grouped barchart of cell proportions
ggplot(cluster_percent_Cell, aes(fill=Treatment, y=Frequency, x=Immgen_sc_labels)) +
geom_bar(position="dodge", stat="identity")+
ggtitle("Cell distribution according to Immgen reference [%]") +
theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=0.5)) +
xlab("")
Renaming Clusters
DimPlot(samples.integrated, reduction = "umap", label=T) + NoLegend()
Zoom into Microglia clusters
A second Seurat cluster is generated, where all non-microglial immune cells are removed according to the mouseRNASeq_sc_labels.
p3 + p3a
DAMs (Cell. 2017 Jun 15;169(7):1276-1290.e17. doi: 10.1016/j.cell.2017.05.018)
DAM plotting
p5 | p5aWarning: Use of `expr_df[[axes[1]]]` is discouraged. Use `.data[[axes[1]]]`
instead.Warning: Use of `expr_df[[axes[2]]]` is discouraged. Use `.data[[axes[2]]]`
instead.
Renaming clusters
Idents(microglia) <- "Celltype"
DimPlot(microglia, reduction = "umap", label = TRUE, pt.size = 1)+ NoLegend()
Signature enrichment
DEGs and GSEA between microglia of different genotypes
WT vs App
WT vs MX04+
MX04+ vs App
DEGs and GSEA of microglia of different genotype and treatment
WT_Ctrl vs WT_Stroke
APP_Ctrl vs APP_Stroke
WT_Ctrl vs APP_Ctrl
WT_Stroke vs APP_Stroke
DEGs and GSEA of each cluster
Microglia_0
Only works if ReloadAllData_Harmony_Hefendehl_Stroke.RData is loaded. Otherwise, an error message will occur.
Microglia_1
Microglia_2
Microglia_4
DAM
Microglia_5
Only works if ReloadAllData_Harmony_Hefendehl_Stroke.RData is loaded. Otherwise, an error message will occur.
Genotype + Treatment comparison in each cluster
Microglia_0
Microglia_1
Microglia_2
DAM?
Microglia_4
Microglia_5
sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
[5] LC_TIME=German_Germany.1252
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] scales_1.1.1 pheatmap_1.0.12
[3] EnhancedVolcano_1.12.0 ggrepel_0.9.1
[5] SingleR_1.8.1 SummarizedExperiment_1.24.0
[7] Biobase_2.54.0 GenomicRanges_1.46.1
[9] GenomeInfoDb_1.30.1 IRanges_2.28.0
[11] S4Vectors_0.32.3 BiocGenerics_0.40.0
[13] MatrixGenerics_1.6.0 matrixStats_0.61.0
[15] ggplot2_3.3.5 SeuratObject_4.0.4
[17] Seurat_4.1.0 workflowr_1.7.0
loaded via a namespace (and not attached):
[1] scattermore_0.7 tidyr_1.1.4
[3] bit64_4.0.5 knitr_1.37
[5] irlba_2.3.5 DelayedArray_0.20.0
[7] data.table_1.14.2 rpart_4.1.16
[9] KEGGREST_1.34.0 RCurl_1.98-1.5
[11] generics_0.1.2 ScaledMatrix_1.2.0
[13] callr_3.7.0 cowplot_1.1.1
[15] RSQLite_2.2.9 shadowtext_0.1.1
[17] RANN_2.6.1 future_1.23.0
[19] bit_4.0.4 enrichplot_1.14.1
[21] spatstat.data_2.1-2 httpuv_1.6.5
[23] assertthat_0.2.1 viridis_0.6.2
[25] xfun_0.29 jquerylib_0.1.4
[27] evaluate_0.14 promises_1.2.0.1
[29] fansi_1.0.2 igraph_1.2.11
[31] DBI_1.1.2 htmlwidgets_1.5.4
[33] spatstat.geom_2.3-1 purrr_0.3.4
[35] ellipsis_0.3.2 dplyr_1.0.7
[37] deldir_1.0-6 sparseMatrixStats_1.6.0
[39] vctrs_0.3.8 Cairo_1.5-14
[41] ROCR_1.0-11 abind_1.4-5
[43] cachem_1.0.6 withr_2.4.3
[45] ggforce_0.3.3 sctransform_0.3.3
[47] treeio_1.18.1 goftest_1.2-3
[49] cluster_2.1.2 DOSE_3.20.1
[51] ape_5.6-1 lazyeval_0.2.2
[53] crayon_1.4.2 pkgconfig_2.0.3
[55] labeling_0.4.2 tweenr_1.0.2
[57] nlme_3.1-155 vipor_0.4.5
[59] rlang_1.0.0 globals_0.14.0
[61] lifecycle_1.0.1 miniUI_0.1.1.1
[63] downloader_0.4 extrafontdb_1.0
[65] rsvd_1.0.5 ggrastr_1.0.1
[67] rprojroot_2.0.2 polyclip_1.10-0
[69] lmtest_0.9-39 Matrix_1.4-0
[71] aplot_0.1.2 zoo_1.8-9
[73] beeswarm_0.4.0 whisker_0.4
[75] ggridges_0.5.3 processx_3.5.2
[77] png_0.1-7 viridisLite_0.4.0
[79] bitops_1.0-7 getPass_0.2-2
[81] KernSmooth_2.23-20 Biostrings_2.62.0
[83] blob_1.2.2 DelayedMatrixStats_1.16.0
[85] stringr_1.4.0 qvalue_2.26.0
[87] parallelly_1.30.0 gridGraphics_0.5-1
[89] beachmat_2.10.0 memoise_2.0.1
[91] magrittr_2.0.2 plyr_1.8.6
[93] ica_1.0-2 zlibbioc_1.40.0
[95] compiler_4.1.2 scatterpie_0.1.7
[97] RColorBrewer_1.1-2 ash_1.0-15
[99] fitdistrplus_1.1-6 cli_3.1.1
[101] XVector_0.34.0 listenv_0.8.0
[103] patchwork_1.1.1 pbapply_1.5-0
[105] ps_1.6.0 MASS_7.3-55
[107] mgcv_1.8-38 tidyselect_1.1.1
[109] stringi_1.7.6 highr_0.9
[111] proj4_1.0-10.1 yaml_2.2.2
[113] GOSemSim_2.20.0 BiocSingular_1.10.0
[115] grid_4.1.2 sass_0.4.0
[117] fastmatch_1.1-3 tools_4.1.2
[119] future.apply_1.8.1 parallel_4.1.2
[121] rstudioapi_0.13 git2r_0.29.0
[123] gridExtra_2.3 farver_2.1.0
[125] Rtsne_0.15 ggraph_2.0.5
[127] digest_0.6.29 shiny_1.7.1
[129] Rcpp_1.0.8 ggalt_0.4.0
[131] later_1.3.0 RcppAnnoy_0.0.19
[133] httr_1.4.2 AnnotationDbi_1.56.2
[135] colorspace_2.0-2 fs_1.5.2
[137] tensor_1.5 reticulate_1.24
[139] splines_4.1.2 uwot_0.1.11
[141] yulab.utils_0.0.4 tidytree_0.3.7
[143] spatstat.utils_2.3-0 graphlayouts_0.8.0
[145] ggplotify_0.1.0 plotly_4.10.0
[147] xtable_1.8-4 jsonlite_1.7.3
[149] ggtree_3.2.1 tidygraph_1.2.0
[151] ggfun_0.0.5 R6_2.5.1
[153] pillar_1.7.0 htmltools_0.5.2
[155] mime_0.12 glue_1.6.1
[157] fastmap_1.1.0 clusterProfiler_4.2.2
[159] BiocParallel_1.28.3 BiocNeighbors_1.12.0
[161] codetools_0.2-18 maps_3.4.0
[163] fgsea_1.20.0 utf8_1.2.2
[165] lattice_0.20-45 bslib_0.3.1
[167] spatstat.sparse_2.1-0 tibble_3.1.6
[169] ggbeeswarm_0.6.0 leiden_0.3.9
[171] GO.db_3.14.0 Rttf2pt1_1.3.9
[173] survival_3.2-13 rmarkdown_2.11
[175] munsell_0.5.0 DO.db_2.9
[177] GenomeInfoDbData_1.2.7 reshape2_1.4.4
[179] gtable_0.3.0 spatstat.core_2.3-2
[181] extrafont_0.17