NEBULA Differential Expression Pipeline

This pipeline performs cell-type–specific differential expression analysis using NEBULA on Seurat objects. It is implemented in Snakemake for reproducibility, scalability, and memory efficiency.

It is designed to be flexible regarding the variables of interest for the differential expression analysis (e.g., compare cases vs controls, males vs females, APOEe4 positive vs APOEe4 negative). It also allows for the inclusion of covariates in the model.

It was designed to be a downstream analysis after cell type annotation and quality control of single-cell RNA-seq data.

The pipeline can handle multiple cell types in one run, and it splits the analysis into gene batches to reduce memory load.

Features

• Works with Seurat objects stored as .rds or .qs
• Supports baseline-centered or all-pairwise comparisons
• Splits analyses into gene batches to reduce memory load
• Runs the differential expression analysis for all genes in a seurat object or an user-provided list of genes
• Optional: automatically installs NEBULA from GitHub inside a conda env
• Produces combined results tables per cell type × comparison. Multiple testing correction is applied to alls genes tested in a given cell type and comparison.

Inputs

• Required: Seurat objects stored in a specified input directory. Each cell type should be stored in a separate Seurat object file. The filenames should match the cell type names specified in the config file.

• Optional: A list of genes to test. If not provided, all genes in the Seurat object will be tested.

Usage

1. Clone this repository:

   git clone https://github.com/HarariLab/nebula_de.git

2. Install the conda environment specified in config/config.yaml (or let the pipeline do it for you if you set install_nebula: true in the config file).

bash setup_environment.sh

3. Modify the config/config.yaml file to specify your parameters, input/output directories, and variables of interest.

• Some important parameters to set in the config file:
  • input_dir: Directory containing the Seurat objects
  • output_dir: Directory to save results
  • gene_split_factor: Split factor to apply to the number of genes in each cell type. For example, if a cell type has 2000 genes and gene_split_factor is set to 4, the analysis will be split into 4 batches of 500 genes each. Set to 1 to disable splitting.
  • cell_types: List of cell types to analyze (should match Seurat object filenames)
  • target_variable: Metadata column for the main variable of interest. Should match the column name in the Seurat object metadata.
  • target_levels: The values of the target variable to be tested. For example, if target_variable is diagnosis, and you want to compare SAD vs Control and ADAD vs Control, set target_levels: ["SAD", "ADAD" "Control"].
  • baseline_level: (Optional) Baseline level for comparisons. In the previous example, if you set baseline_level: "Control", the pipeline will set the Control level as reference level. This sets the direction of the log fold changes (upregulated genes in SAD and ADAD will have positive logFC values). In this case, the pipeline will also create the non-baseline pairwise comparison SAD vs ADAD. But the order of this comparison will be lexographically determined. If you want to create all pairwise comparisons without a baseline, set baseline_level: null or leave it empty. In this case, the comparisons will also be lexographically determined.
  • covariates: (Optional) List of covariates to include in the model. Should match the column names in the Seurat object metadata.
  • id_column: Column in the metadata for individual IDs (used as random effect in NEBULA)
  • nebula_params: Parameters for NEBULA (e.g., min_cells, min_expression, fdr_threshold)
  • stratify_by / stratify_levels: (Optional) Run the DE analysis separately within each level of a metadata variable, in addition to the unstratified run. Cells are subset to stratify_by == level before the comparison is applied, so each stratum tests the target_variable contrast within that subgroup only. stratify_levels must match the metadata values exactly; see Outputs for how they map onto directory names.
  • interactions: (Optional) List of interaction terms to add to the model — one per pipeline run. Each entry is a single interaction such as "comparison_group:Gender" or "comparison_group:Gender:APOE". Use the literal token comparison_group to refer to the target-variable contrast. For each interaction, all main effects and lower-order interactions among the participating variables are added automatically (principle of marginality). Cannot be combined with stratify_by — set one or the other.

4. Run the pipeline with Snakemake:

Pipeline wrote and tested with Snakemake v9.3.2.

Two profiles were included: local and slurm. These can be configured as desired. The local profile is intended for local execution (e.g. DORA) and the slurm is for slurm scheduler (e.g. OSC).

# On the root directory of the cloned repository, choose a profile:

# Local execution (no scheduler):
snakemake --workflow-profile profiles/local

# SLURM cluster:
snakemake --workflow-profile profiles/slurm

Each profile (profiles/local/config.yaml and profiles/slurm/config.yaml) sets the executor, conda usage, parallelism, and per-rule resources. Edit the profile file to adjust cores, memory, runtime, or the SLURM account.

Gene filtering

Two filters are applied before fitting the model:

1. Pre-NEBULA filter (configurable via nebula_params.min_cells and nebula_params.min_expression): a gene is kept if at least min_cells cells have raw count > min_expression. With the template defaults (min_cells: 10, min_expression: 0.1) this is effectively "≥ 10 cells with at least 1 read."
2. NEBULA-internal filter (hardcoded defaults mincp = 5, cpc = 0.005): a gene must additionally have ≥ 5 cells with a non-zero count and a mean count per cell > 0.005. The mean threshold scales with the number of cells, so for large cell-type matrices this is the stricter cutoff.

Genes that pass (1) but fail (2) appear in the NEBULA log as Remove N genes having low expression. — they're dropped silently from the results.

Outputs

The results will be saved in the specified output_dir under the path final_results/{stratum}/{comparison}/ (or final_results/{stratum}/{interaction}/{comparison}/ when interactions is configured). When interactions is set, one directory per configured interaction is created; when unset, the interaction segment is omitted from the path entirely.

{stratum} names both the variable and the level it came from, so results from runs with different variables of interest sit side by side instead of overwriting each other. The unstratified run is named after the target_variable being tested; spaces and other path-hostile characters are stripped from the level names:

Config	Stratum directories
target_variable: Sex, stratify_by: Case_sub_status, stratify_levels: ["sporadic", "CT old"]	all_Sex/, Case_sub_status_sporadic/, Case_sub_status_CTold/
target_variable: Case_sub_status, stratify_by: Sex, stratify_levels: ["Female", "Male"]	all_Case_sub_status/, Sex_Female/, Sex_Male/
target_variable: Sex, stratify_by: null	all_Sex/

Snakemake prints the full directory → cell-subset mapping at the start of every run.

For each cell type, comparison, stratum and interaction, three files are written:

• *_final_de_results.csv — combined results with the following columns:
  • gene: Gene name
  • log2FC*: Log2FC values for the variables included in the model (target variable and covariates). Nebula returns the logFC (in ln scale), so we convert it to log2 scale.
  • p_*: P-values for every model term (main effects, covariates, interaction if configured).
  • neg_log10_p_*: -log10 transformed p-values for the variables included in the model.
  • abs_log2FC*: Absolute log2FC values.
  • padj_* / fdr_*: Adjusted p-values for every term.
  • significant: Boolean flag for the main contrast (target-variable effect).
  • significant_interaction: Boolean flag for the interaction term (only meaningful when interactions is configured).
• *_significant_genes.txt — genes significant for the main contrast.
• *_significant_genes_interaction.txt — genes significant for the interaction term (empty for runs without an interaction, or when the interaction does not involve the target variable).