Last updated: 2026-04-23
Checks: 2 0
Knit directory: scSeq_Hefendehl/
This reproducible R Markdown analysis was created with workflowr (version 1.7.2). 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! 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 23714d5. 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: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/archive/
Ignored: analysis/output/
Ignored: data
Ignored: output/Res_202602/
Ignored: output/archive/
Unstaged changes:
Modified: wflowhelper.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.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/index.Rmd) and HTML
(docs/index.html) files. If you’ve configured a remote Git
repository (see ?wflow_git_remote), click on the hyperlinks
in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 3475e86 | achiocch | 2026-04-13 | Build site. |
| Rmd | 375b529 | achiocch | 2026-03-20 | index update |
| html | 375b529 | achiocch | 2026-03-20 | index update |
| html | 75e9a44 | achiocch | 2026-03-19 | Build site. |
| Rmd | f7ac121 | achiocch | 2026-03-19 | Updated Analysis cluster DEX and GO |
| html | f7ac121 | achiocch | 2026-03-19 | Updated Analysis cluster DEX and GO |
| html | 0d6e2cb | achiocch | 2026-03-18 | Build site. |
| Rmd | 05a87a1 | achiocch | 2026-03-18 | workflowr::wflow_publish(c("./analysis/sc", "./analysis/index")) |
| html | 8f13584 | achiocch | 2024-02-12 | Build site. |
| html | 23643ee | achiocch | 2024-02-09 | Build site. |
| Rmd | dcc077d | achiocch | 2024-02-09 | workflowr::wflow_publish(c("./analysis/03", "./analysis/index")) |
| Rmd | 2290b23 | Andreas Chiocchetti | 2023-11-03 | backup 20231103 |
| html | 1f4af15 | Andreas Chiocchetti | 2023-08-11 | Build site. |
| html | 87fbacb | Andreas Chiocchetti | 2023-08-07 | Build site. |
| html | 6d3e973 | Andreas Chiocchetti | 2023-07-14 | Build site. |
| html | 982b9f6 | Andreas Chiocchetti | 2023-02-20 | Build site. |
| Rmd | 438f7ad | Andreas Chiocchetti | 2023-02-20 | workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*", |
| html | 0e20f22 | Andreas Chiocchetti | 2023-02-16 | Build site. |
| Rmd | eaa373d | Andreas Chiocchetti | 2023-02-16 | workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*", |
| html | 7764484 | achiocch | 2022-06-27 | Build site. |
| html | 14e1ee0 | achiocch | 2022-06-27 | Build site. |
| html | 3451e15 | achiocch | 2022-06-22 | Build site. |
| Rmd | cfbfcf6 | achiocch | 2022-06-22 | wflow_publish(c("analysis/", "docs/", "code/*")) |
| html | cfbfcf6 | achiocch | 2022-06-22 | wflow_publish(c("analysis/", "docs/", "code/*")) |
| html | faf20e9 | achiocch | 2022-05-16 | Build site. |
| html | 51b35c8 | achiocch | 2022-04-26 | Build site. |
| html | 41d6cd0 | achiocch | 2022-04-26 | Build site. |
| html | 30f02eb | achiocch | 2022-04-22 | Build site. |
| html | 9b778b1 | achiocch | 2022-04-19 | Build site. |
| html | 17114be | achiocch | 2022-04-08 | Build site. |
| html | cf395f4 | achiocch | 2022-03-30 | Build site. |
| html | c998828 | achiocch | 2022-03-30 | Build site. |
| html | 87439a3 | achiocch | 2022-02-21 | Build site. |
| html | d6f2105 | achiocch | 2022-02-21 | Build site. |
| html | 79dd0b1 | achiocch | 2022-02-21 | Build site. |
| html | ded601d | achiocch | 2022-02-03 | Build site. |
| html | 50fe211 | achiocch | 2022-02-03 | Build site. |
| Rmd | 64cd37c | achiocch | 2022-02-03 | wflow_publish(c("analysis/", "docs/", "code/*")) |
| html | 5fbdba4 | achiocch | 2022-02-02 | Build site. |
| html | 9044798 | Andreas Geburtig-Chiocchetti | 2022-02-02 | workflowr automated |
| Rmd | 57c9e70 | achiocch | 2022-02-02 | adds initial html files |
| html | 57c9e70 | achiocch | 2022-02-02 | adds initial html files |
| Rmd | ddb6373 | achiocch | 2022-02-02 | first commit |
| html | f010e79 | achiocch | 2022-01-31 | Build site. |
| html | 4c4711a | achiocch | 2022-01-31 | wflow_publish(c("analysis/", "docs/", "code/*")) |
| Rmd | 730bee8 | achiocch | 2022-01-31 | Start workflowr project. |
‘# These authors contributed equally to this work’
‘* Corresponding author:’ hefendehl@bio.uni-frankfurt.de
Single-cell RNA sequencing analysis of microglial responses in a
mouse model of Alzheimer’s disease (APPPS1) with and without ischemic
stroke. All analyses are implemented in a single R Markdown document
(scRNA_Analyses_Candlishetal.Rmd).
The following sections gives an overview and explanation of the usage and different steps within the analytcal pipeline
This pipeline processes plate-based Smart-seq2 single-cell RNA-seq data from mouse brain immune cells (primarily microglia) across two sequencing batches. The goal is to characterize microglial heterogeneity and transcriptional responses in the context of amyloid pathology (APPPS1 genotype) and acute ischemic stroke, alone and in combination.
| Group | Genotype | Treatment |
|---|---|---|
| WT_Ctrl | Wild-type | Control |
| WT_Stroke | Wild-type | Ischemic stroke |
| APPPS1_Ctrl | APPPS1 (amyloid model) | Control |
| APPPS1_Stroke | APPPS1 (amyloid model) | Ischemic stroke |
Data were collected across two sequencing batches and integrated prior to downstream analysis.
All analyses run in R. The following packages are required:
Core single-cell analysis: Seurat,
SingleCellExperiment, scCustomize,
harmony, SingleR, celldex,
slingshot, TSCAN
Differential expression & statistics:
DESeq2, limma, EnhancedVolcano,
compareGroups, chisq.posthoc.test,
mclust
Gene co-expression: WGCNA
Annotation & enrichment:
AnnotationHub, org.Mm.eg.db,
gprofiler2 (via custom getGOresults()
wrapper)
Data wrangling & visualization:
tidyverse, data.table, openxlsx,
ggpubr, ggrepel, viridis,
plotly, clustree, DT,
knitr
A custom function library (./code/custom_functions.R)
provides helper functions including display_tab(),
getGOresults(), GOplot(),
gg_qqplot(), and convertpvaltostars().
Mouse gene annotations are fetched from Ensembl via
AnnotationHub, retrieving the most recent
EnsDb for Mus musculus. These are used throughout
the pipeline to map Ensembl gene IDs to gene symbols. Cell cycle gene
lists (S-phase and G2/M-phase) are downloaded from the Harvard
Bioinformatics Core’s TinyAtlas repository.
Batch 2 (second sequencing run): - Sample metadata
(plate, well, genotype, treatment, sex, age, mouse ID) are read from an
Excel file and expanded from well-range notation (e.g.,
A1–P12) into individual cell-level entries. - Kallisto
count matrices are loaded per plate, Ensembl IDs are mapped to gene
symbols, and duplicate gene names are resolved by summing counts. Each
plate is assembled into a Seurat object.
Batch 1 (first sequencing run): - Count matrices from two Kallisto alignment runs are merged, gene IDs are converted to symbols, and metadata are joined from per-plate sample tables. A single combined Seurat object is created and split by plate.
Both batches are then merged into a single Seurat object for joint processing.
Mitochondrial read percentage (percent.mt) is calculated
for all cells. Cells are retained if they pass all of the following
thresholds:
Genes are additionally filtered to retain only those expressed in ≥ 0.5% of cells. Pre- and post-filter QC distributions are visualized as violin plots grouped by sequencing plate.
The merged object is normalized (NormalizeData),
variable features are identified, and the data are scaled and reduced
via PCA. An initial unintegrated UMAP is generated to visualize batch
effects.
Integration is performed using CCA (Canonical Correlation
Analysis) via IntegrateLayers with 30 PCs,
followed by graph-based clustering and UMAP embedding. The integrated
UMAP is inspected for plate effects and expression of canonical immune
markers (e.g., Aif1, Ptprc).
A second integration pass using Harmony
(RunHarmony, grouping by sequencing run) is applied after
isolating microglia (Step 6) for the final microglia-specific
analysis.
Automated cell type annotation is performed using
SingleR against the MouseRNAseqData
reference from celldex. Fine-grained cell type labels are
assigned to all cells and visualized on the integrated UMAP, split by
plate and experimental condition.
Microglia isolation: Cells annotated as microglia by
SingleR (Immgen_sc_labels_agc containing “Microglia”) are
extracted and carried forward for all subsequent analyses.
Cell cycle scoring (S-phase and G2/M scores) is applied to the microglia subset using the downloaded cell cycle gene lists.
The microglia subset is re-normalized, re-scaled, and re-embedded (PCA → Harmony → CCA integration → UMAP) specifically to resolve microglial subpopulations.
Resolution optimization: Clustering is run across a
range of resolutions (0.2–1.2) and stability is assessed using the
Adjusted Rand Index (ARI) between successive
resolutions. The optimal resolution (selected as 0.6) is applied to
produce the final clusters, labeled Microglia0 through
MicrogliaX.
UMAP visualizations are produced grouped by cluster, plate,
condition, mouse ID, sex, and homeostatic marker expression. A cluster
tree (clustree) is saved to document resolution
stability.
To facilitate comparison with prior published analyses, cluster labels from the initial Batch 1–only analysis are mapped onto the integrated UMAP.
Module scores are computed with AddModuleScore for the
following microglial states using published gene signatures:
| State | Key Genes |
|---|---|
| DAM (Disease-Associated Microglia) | Itgax, Cst7, Apoe,
Trem2, Lpl, Axl,
Clec7a, … |
| DAM stage 1 (DAM1) | Tyrobp, Ctsb, Ctsd,
Apoe, B2m, Fth1,
Lyz2 |
| DAM stage 2 (DAM2) | Trem2, Axl, Cst7,
Ctsl, Lpl, Cd9,
Csf1, … |
| Homeostatic | P2ry12, P2ry13, Tmem119,
Cx3cr1, Selplg, Cd33 |
| Interferon-responsive | Ifit2, Ifit3, Ifitm3,
Irf7, Oasl2 |
| Cycling | Top2a, Mcm2, Tubb5,
Mki67, Cdk1 |
| Activation-responsive | Cd74, H2-Ab1, H2-Aa,
Ctsb, Ctsd |
| Axon-tract associated | Spp1, Gpnmb, Igf1,
Lgals3, Fabp5, … |
| BAM (Border-Associated Macrophages) | Cd36, Cd38, Lyve1,
Cd163, Cd169 |
| CAM (CNS-Associated Macrophages) | Emilin2, Pf4, Hp,
F5, Mki67 |
Additionally, spatial gene signatures derived from companion spatial transcriptomic data (peri-ictus, distal, and control regions) are scored and binarized using density-derived cutoffs. Each continuous score is converted to a binary classification, and proportions are compared across conditions and clusters using chi-square tests with post-hoc testing.
A combined Homeostatic vs. DAM classification assigns each cell to one of four mutually exclusive states: Homeostatic, DAM, Both, or Neither.
Cluster marker discovery:
FindAllMarkers (Wilcoxon rank-sum, log2FC threshold = 1,
positive markers only) is run across all microglia clusters. GO term
enrichment is computed on per-cluster marker gene sets using a custom
gprofiler2-based wrapper.
Condition-level DEG: Pairwise DEG comparisons are performed between experimental conditions (Wilcoxon, no fold-change pre-filter) for the following contrasts applied within each microglia cluster:
A reusable compareCluster() function wraps each
comparison, handling: FindMarkers execution with error
logging, volcano plots (EnhancedVolcano) combined with Q-Q
plots, and GO enrichment with per-direction gene sets (up-regulated,
down-regulated, all significant).
Cross-condition scatter plots: Genes significantly
regulated in both WT_Stroke and APPPS1_Stroke (vs. respective controls)
are plotted as log2FC scatter plots to identify shared and
condition-specific transcriptional responses. Selected genes of interest
(e.g., S100a6, Apoc4, C4b,
Fn1) are labelled.
Dot plots: Scaled mean expression and percentage-expressed are visualized as bubble dot plots across cluster × condition for DAM marker genes and spatial gene sets.
All DEG results and GO term tables are exported to Excel
(.xlsx).
Trajectory inference is performed with Slingshot on
the integrated UMAP embedding, using Microglia0 as the root
(start cluster) and Microglia3, Microglia4,
and Microglia5 as terminal states. Multiple lineages are
inferred with omega-scaling to prevent trajectory extension.
Pseudotime values for each lineage are transferred back to the Seurat object metadata and overlaid on UMAPs.
Trajectory-associated genes are identified using
TSCAN (testPseudotime) for each lineage
independently. The top 10 up- and down-regulated genes along each
trajectory are visualized as expression-along-pseudotime plots, with
cells colored by cluster label.
Results (ranked gene tables with log-fold changes and FDR values) are exported per lineage to Excel.
Weighted Gene Co-expression Network Analysis (WGCNA) is applied to the microglia expression matrix to identify co-regulated gene modules.
Pre-processing: The expression matrix is filtered to
genes expressed in > 10% of cells and genes with annotated mouse
Entrez IDs (org.Mm.eg.db). Network topology analysis
(pickSoftThreshold) is used to select the soft-thresholding
power (set to 3).
Network construction: A consensus blockwise network
is built with blockwiseConsensusModules (unsigned network,
Pearson correlation, maximum block size 30,000 genes). Modules are
assigned color labels.
Module-trait associations: Module eigengenes (MEs)
are extracted and their association with genotype × treatment
interaction is tested using linear models
(lm(value ~ Genotype * Treatment)). ANOVA is used to test
cluster-level associations. Results are visualized as labeled heatmaps
(mean ± SD eigengene per condition/cluster) and violin plots, with
significance annotated by p-value stars.
Gene-level heatmaps: For each module, heatmaps of the top 50 most highly expressed member genes are generated grouped by both condition and cluster.
GO enrichment: GO term enrichment is computed for
non-grey modules (modules with significant cluster or condition
associations) using gprofiler2.
The following additional visualizations are produced throughout the pipeline:
MG_score_*) across clusters.barplot_genes(), allowing visualization of any gene split
by condition, cluster, or both.All results are written to ./output/Res_202602/. Key
output categories include:
| Category | Format |
|---|---|
| UMAP / DimPlot visualizations | |
| Feature plots, ridge plots, violin plots | |
| Volcano plots + Q-Q plots (per DEG comparison) | |
| GO term bar charts (per comparison/module) | |
| WGCNA heatmaps and eigengene plots | |
| Pseudotime trajectory plots and gene expression plots | |
| Cluster marker gene tables | .xlsx |
| DEG result tables (per comparison) | .xlsx |
| GO enrichment tables | .xlsx |
| WGCNA linear model statistics | .xlsx |
| Cell count / cluster proportion tables | .xlsx |
| Processed Seurat object | .rds |
A global random seed (set.seed(1324567)) is set at the
start of the document. The WGCNA section is computationally intensive
(can take several hours); results are cached to .rds files
and reloaded on subsequent runs via force.recalc flags.
Intermediate Seurat objects are saved and can be loaded directly to
resume from mid-pipeline checkpoints.