A Snakemake 8 workflow for performing differential expression analyses (DEA) of processed (multimodal) scRNA-seq data powered by the R package Seurat's functions FindMarkers and FindAllMarkers.

• Stephan Reichl
• Christoph Bock

This project wouldn't be possible without the following software and their dependencies:

This is a template for the Methods section of a scientific publication and is intended to serve as a starting point. Only retain paragraphs relevant to your analysis. References [ref] to the respective publications are curated in the software table above. Versions (ver) have to be read out from the respective conda environment specifications (workflow/envs/*.yaml file) or post-execution in the result directory (dea_seurat/envs/*.yaml). Parameters that have to be adapted depending on the data or workflow configurations are denoted in squared brackets e.g., [X].

The outlined analyses were performed using the R package Seurat (ver) [ref] unless stated otherwise.

Differential Expression Analysis (DEA). DEA was performed on the assay [X] and data slot [X] with Seurat's [FindMarkers|FindAllMarkers] function using the statistical test [X] with the parameters log2(fold change) threshold of [X] and minimal percentage of expression [X]. The results were filtered for relevant features by adjusted p-value of [X], absolute log2(fold change) of [X] and minimum percentage of expression [X].

Visualization. All filtered result statistics, i.e., number of statistically significant results split by positive (up) and negative (down) effect-sizes, were separately visualized with stacked bar plots using ggplot (ver) [ref]. To visually summarize results of the same analysis the filtered log2(fold change) values of features that were found to be at least in one comparison statistically significantly differentially expressed were visualized in a hierarchically clustered heatmap using pheatmap (ver) [ref]. Volcano plots were generated for each analysis using EnhancedVolcano (ver) [ref] with adjusted p-value threshold of [X] and log2(fold change) threshold of [X] as visual cut-offs for the y- and x-axis, respectively.

The analysis and visualizations described here were performed using a publicly available Snakemake [ver] (ref) workflow [10.5281/zenodo.10689139].

The workflow performs the following steps to produce the outlined results (dea_seurat/{analysis}/).

• Differential Expression Analysis (DEA)
  • using Seurat's FindMarkers or FindAllMarkers depending on the configuration (results.csv). This step is parallelized using the R package future.
  • feature list per comparison group and direction (up/down) for downstream analysis (e.g., enrichment analysis) (feature_lists/{group}_{up/down}_features.txt)
  • (optional) feature score tables (with two columns: "feature" and "score") per comparison group using {score_formula} for downstream analyses (e.g., preranked enrichment analysis) (feature_lists/{group}_featureScores.csv).
  • results filtered according to the configured thresholds:
    • statistical significance (adjusted p-value).
    • effect-size (log2 fold change).
    • expression (minimum percentage of expression) in one of the comparison groups.
  • filtered result statistics: number of statistically significant results split by positive (up) and negative (down) change (stats.csv).
• Visualizations (dea_seurat/{analysis}/plots)
  • filtered result statistics: number of features and direction as bar plot (stats.png).
  • volanco plots per comparison group with effect size on the x-axis and raw p-value(*_rawp)/adjusted p-value (_adjp) on the y-axis (volcano/{feature_list}/{group}.png).
    • highlighting features according to configured cut-offs for statistical significance (pCutoff) and effect size (FCcutoff).
    • (optional) highlighting features according to configured feature lists.
  • hierarchically clustered heatmap of effect sizes (LFC) per comparison (features x comparisons) indicating statistical significance with a star '*' (heatmap/{feature_list}.png).
    • using all filtered features (FILTERED)
    • (optional) using configured feature lists
    • in case of more than 100 features the row labels and significance indicators (*) are removed

Here are some tips for the usage of this workflow:

• Perform your first run with loose filtering options/cut-offs and set the same for filtering and plotting to see if further filtering is even necessar or useful.
• Try one small/simple analysis first before running all desired analyses.

Detailed specifications can be found here ./config/README.md

We selected a scRNA-seq data set consisting of 15 CRC samples from Lee et al (2020) Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nature Genetics. Downloaded from the Weizmann Institute - Curated Cancer Cell Atlas (3CA) - Colorectal Cancer section.

• samples/patients: 15
• cells: 21657
• features (genes): 22276
• preprocessed using the compatible MrBiomics module for scRNA-seq data processing & visualization
• We performed a 1 vs rest analysis using the cell type annotation ("ALL").
• total runtime on HPC w/ SLURM (32GB RAM; only DEA with 8 cores otherwise 1 core): <25 minutes for 17 jobs in total

A comparison of the cell type marker expression split by cell types visualized as a dot plot with the DEA results as hierarchically clustered heatmap of the effect sizes.

data source/authors	Workflow Output

We provide metadata, annotation and configuration files for this data set in ./test. The processed and prepared Seurat RDS object has to be downloaded from Zenodo by following the instructions below.

# download Zenodo records using zenodo_get

# install zenodo_get v1.3.4
conda install -c conda-forge zenodo_get=1.3.4

# download the prepare Seurat RDS object
zenodo_get --record 10688824 --output-dir=test/data/Lee2020NatGenet/

• GitHub Repository
• GitHub Page
• Zenodo Repository
• Snakemake Workflow Catalog Entry

• Recommended compatible MrBiomics modules for:
  • upstream analysis
    • scRNA-seq Data Processing & Visualization for processing (multimodal) single-cell transcriptome data.
  • downstream analyses
    • Unsupervised Analysis to understand and visualize similarities and variations between cells/samples, including dimensionality reduction and cluster analysis. Useful for all tabular data including single-cell and bulk sequencing data.
  • Enrichment Analysis for biomedical interpretation of (differential) analysis results using prior knowledge.

The following publications successfully used this module for their analyses.

• FirstAuthors et al. (202X) Journal Name - Paper Title.
• ...