Chapter 10 — The Talapas Analysis Pipeline

Author

Single Cell RNA-seq Workshop

Note

Where this chapter sits. This is the conceptual reading for the second of the two Friday-morning Talapas modules. Companion to Lecture 10 — Running the Analysis Pipeline on Talapas and Tutorial 10 — Talapas Pipeline. Prerequisite: Chapter 9 — VS Code & SLURM Basics on Talapas, which covers first-time login, VS Code Remote-SSH setup, Lmod modules, Talapas file systems, reproducible environments, partitions, and sbatch. Where Chapter 9 taught the mechanics of connecting to and operating the cluster, this chapter explains the batch-mode architecture of the scRNA-seq analysis pipeline you run on it.

10.1 From notebooks to scripts: the fourth leg of 1:1:1:1

Every module in this workshop is taught through four parallel artifacts — a lecture, a reading, a laptop tutorial, and a Talapas script — all covering the same workflow at different altitudes. The laptop tutorials (notebooks 01–08) are learning documents: rich prose, “Think about it” prompts, inline plots, and explanation woven through every step. They are written to be read and understood at a desk, one cell at a time.

The Talapas scripts are the fourth leg. They run the same analysis, on the same dataset, but stripped down to runnable R with only minimal comments that mark which notebook each step parallels. This division of labor is deliberate. A teaching notebook and a production script have opposite goals: the notebook maximizes explanation, the script maximizes reproducibility and unattended execution. Rather than compromise either, the workshop keeps them separate — the explanation lives in the notebook, the execution lives in the script — and the comment markers in the script point you back to the notebook whenever you need the “why.”

The payoff is that moving to the cluster introduces no new analysis to learn. You already understand QC, clustering, annotation, integration, pseudobulk DE, functional enrichment, and differential abundance from the laptop series. Chapter 10 is purely about how that same understanding is packaged so it can run as a chain of SLURM jobs on Talapas.

The practical difference between a notebook and a script shows up in two concrete ways. First, a notebook runs cell by cell in RStudio; you see each plot and table as it is produced, and you can pause, inspect, and change course. A script runs from start to finish as a single process: either it completes successfully and writes its output, or it fails. There is no mid-run intervention. This is why the scripts print explicit progress messages (“==== RUNNING 01_qc_preprocessing.R ====”) and end with cat("Wrote ...") — those messages are your only window into what happened, and reading the .out log is how you confirm success. Second, a script on a shared cluster must be self-contained: it cannot depend on variables in your current R session, on packages installed in an interactive session that wasn’t loaded by module load, or on data that happen to be in your current working directory. Everything it needs — the package list, the data path, the parameters — must be explicit and inside the script itself.

10.2 The dataset: why ifnb is not a toy

The scripts run the ifnb dataset1, the same one used throughout the laptop tutorials. “ifnb” is the familiar name for these Kang et al. PBMCs; the pipeline actually loads them via muscData::Kang18_8vs8() — a Bioconductor ExperimentHub object (the full ~29,000-cell build, ~24,000 annotated singlets across 8 donors) that carries the per-cell donor IDs the older SeuratData ifnb lacks. It is worth dwelling on why this dataset was chosen, because the choice is what makes the cluster run a real analysis rather than a demonstration.

ifnb carries two kinds of structure that most single-cell demo datasets lack:

A real biological condition — stim, with two levels: CTRL (control) and STIM (cells stimulated with interferon-β) — this is a genuine perturbation with a known biological consequence. And real biological replicates — ind, identifying 8 lupus patient donors, each measured in both conditions, providing true replicates at the level that matters for statistics: the donor.

Contrast this with a single dissociated-tumor run, which gives you one sample, no condition, and no replicates. With such data you can cluster and annotate cells, but you cannot ask whether anything changes between conditions, because there is no condition and nothing to compare. ifnb has both axes by design, so the downstream steps — integration (05), pseudobulk differential expression (06), and differential abundance (08) — operate on genuine sample and donor structure. Integration corrects the CTRL/STIM batch axis; the DE and DA tests recover a real interferon-stimulated gene (ISG) response rather than returning a null result. The cluster run reproduces the workshop’s actual biological findings, at scale, on the real machine.

10.3 Pipeline architecture: a chain of object hand-offs

The eight core scripts are not eight independent programs. They form a pipeline — a chain in which each step’s output is the next step’s input. The connective tissue is the filesystem: each script reads an object from disk, does its work, and writes a new object that the following script will read. This pattern has a formal name in software engineering — it is a dataflow pipeline, where each stage is connected to the next by a named artifact rather than by a direct function call. The advantage is that stages can be developed, debugged, and re-run independently; the disadvantage is that you bear responsibility for verifying each artifact before passing it forward.

Table 1: Core pipeline scripts — input and output of each step. The chain reads left to right: each script reads the .rds written by the preceding step and writes a new one. Steps 06 and 08 both consume ifnb_integrated.rds independently, so they can run in parallel once step 05 completes.
Script Parallels notebook Reads Writes
01_qc_preprocessing.R Tutorial 01 — QC & Preprocessing muscData::Kang18_8vs8() ifnb_preprocessed.rds
02_dimreduction_clustering.R Tutorial 02 — DimReduction & Clustering ifnb_preprocessed.rds ifnb_clustered.rds
03_markers_annotation.R Tutorial 03 — Markers & Annotation ifnb_clustered.rds ifnb_annotated.rds
04_reference_annotation.R Tutorial 04 — Reference Annotation ifnb_annotated.rds ifnb_annotated_final.rds
05_integration.R Tutorial 05 — Integration ifnb_annotated_final.rds ifnb_integrated.rds
06_pseudobulk_de.R Tutorial 06 — Pseudobulk DE ifnb_integrated.rds ifnb_pseudobulk_de.csv
07_functional_analysis.R Tutorial 07 — Functional Analysis ifnb_pseudobulk_de.csv functional/*.csv
08_differential_abundance.R Tutorial 08 — Differential Abundance ifnb_integrated.rds ifnb_milo_da.csv

Reading Table 1 top to bottom traces the entire workflow: raw ifnb enters at step 01, becomes a QC’d object, then a clustered object, then an annotated object, then an integrated object — and from the integrated object the pipeline branches to pseudobulk DE (06, which 07 then interprets functionally) and to differential abundance (08). The only external input to the whole chain is ifnb itself, fetched once by script 01. Everything else is an internal hand-off.

Because the steps are linked by their inputs and outputs, they must run in order. Step 02 cannot run until step 01 has written ifnb_preprocessed.rds; step 06 cannot run until step 05 has written ifnb_integrated.rds. A missing or corrupt hand-off breaks the chain at that point. Notice also that steps 06 and 08 both read from ifnb_integrated.rds — they have no dependency on each other, so once step 05 finishes they can run simultaneously if you are using the per-step wrapper (run_rscript.sbatch). The run_core_pipeline.sbatch wrapper runs all eight sequentially, so it cannot exploit this parallelism; the gain from submitting 06 and 08 in parallel is modest for a single-dataset run but becomes meaningful when processing many datasets.

Note🔑 Key concept — why .rds checkpoints make a pipeline robust

Each .rds file is not just a saved object — it is a checkpoint: a self-contained, on-disk snapshot of the analysis state at a specific stage. This pattern decouples the stages from each other. If step 06 (pseudobulk DE) fails because of a memory spike, you do not lose steps 01–05. You fix the memory request, resubmit step 06, and it reads the already-existing ifnb_integrated.rds as if nothing happened. Without checkpoints, a pipeline failure anywhere forces a complete re-run from scratch. The cost of checkpointing is disk space (each Seurat object is 100–500 MB) and the discipline of verifying each checkpoint before passing it forward — a truncated, corrupt .rds from an OOM-killed job looks like a file but is unreadable. The verification habit (check seff for State: COMPLETED, then look for the Wrote ... line in the .out log) is what keeps the chain clean.

10.4 Why .rds checkpoints?

The hand-off format is R’s native serialization, .rds (with .csv for the tabular results of steps 06–08). Saving a fully-formed Seurat object to a .rds file after each step is more than a convenience; it is a deliberate design pattern with several payoffs.

Each .rds is a checkpoint. A step’s output is a complete, self-contained object on disk. If step 06 fails, you do not re-run steps 01 through 05 — you fix step 06 and point it at the already-saved ifnb_integrated.rds. Steps can be right-sized independently. Because each is its own SLURM job reading from disk, you can give step 06 (pseudobulk DE) and step 08 (differential abundance) the extra memory they need without inflating the cheaper steps. Right-sizing per step (the skill from Chapter 9) only works because the pipeline is decomposed this way. The contract between steps is inspectable. The interface from one step to the next is a named file with a known structure. You can open it, examine it, or hand it to a different downstream analysis. This is exactly how reproducible pipelines are built — explicit artifacts, not hidden in-memory state.

The cost of this pattern is the discipline it demands: you must confirm each step truly succeeded before launching the next. A job that runs out of memory or time mid-write leaves behind a truncated, corrupt .rds — a file that exists, so the next step happily tries to read it, but is incomplete, so the read fails in some confusing downstream way. The fix is to verify completion (covered in §10.5c) before chaining onward.

This explicit-checkpoint pattern is a hand-rolled version of what dedicated workflow managers automate: Snakemake2 and Nextflow3 track the file-level dependency graph for you, re-running only the steps whose inputs changed and parallelizing independent branches. For an eight-step teaching pipeline, shell wrappers are clear and sufficient; once a pipeline grows or must run on many datasets, moving to one of these managers — combined with per-tool environments from Bioconda4 and containers such as Singularity/Apptainer5 for full reproducibility — is the natural next step.

Note🔑 Key concept — the distinction between a wrapper script and a workflow manager

The run_rscript.sbatch wrapper and run_core_pipeline.sbatch scripts in the workshop are shell-level orchestration: they submit or loop over the R scripts, but they do not track dependencies, detect partial failures, or skip already-completed steps. If step 03 fails and you resubmit run_core_pipeline.sbatch, it re-runs steps 01 and 02 unnecessarily. A dedicated workflow manager like Snakemake or Nextflow represents the pipeline as an explicit dependency graph — it knows that step 03 reads from step 02’s output file, and it re-runs step 03 only when that file is newer than step 03’s output or missing. For a small teaching pipeline this difference is a convenience, not a necessity. For a production pipeline processing dozens of samples, the ability to automatically skip completed steps and restart only from the point of failure saves hours of compute time and is worth the learning investment.

Note🔑 Key concept — experimental design requirements for differential analysis

The pseudobulk DE (step 06) and differential abundance (step 08) steps require two design features that not all datasets have: a condition (something that differs between groups of samples) and biological replicates (multiple independent samples per condition). The ifnb dataset provides both — stim is the condition (CTRL vs IFN-β) and ind identifies 8 independent lupus donors each measured in both conditions. Without replicates, pseudobulk DE collapses to a single-replicate comparison, and no valid statistical test exists for that design. Without a condition, there is nothing to compare. When you adapt the pipeline to your own data (§10.8), check these two requirements first: if either is absent, steps 06 and 08 will run without error but will return meaningless or null results, not a statistical failure message. The design check belongs in your experimental planning, not in your computational analysis.

10.5 Reproducibility: pre-caching and DATA_DIR

Two practical knobs make the pipeline reproducible and portable on a cluster.

Pre-caching ifnb. Script 01 downloads ifnb (about 25 MB) through Bioconductor’s ExperimentHub the first time it runs. On Talapas this collides with a quirk from Chapter 9: compute nodes may have no outbound internet. A download that works fine on your laptop fails inside the batch job, because the compute node cannot reach the internet to fetch the data. The fix is to populate the cache once, on a login node (which does have internet), before you submit any jobs:

Rscript -e 'invisible(muscData::Kang18_8vs8())'

After this, the data sit in your local ExperimentHub cache, and script 01 finds them there with no network access required.

DATA_DIR — controlling where the hand-off files land. By default the scripts write their .rds/.csv hand-offs into ../data/. The DATA_DIR environment variable overrides that destination, which matters for two reasons. First, the intermediate Seurat objects can be large, and Chapter 9’s etiquette says large working data belongs in /scratch/<PIRG>/, not in /home/. Second, making the output location a variable rather than a hard-coded path is itself a reproducibility practice — the same scripts run unchanged whether the data go to your project directory, scratch, or somewhere else entirely:

export DATA_DIR=/scratch/<PIRG>/<duckid>/ifnb/data

10.5b Design consideration: one allocation vs per-step jobs

The two wrapper scripts embody opposite trade-offs. run_rscript.sbatch gives fine control — each step is its own SLURM job, can request its own memory and time, and can be re-run individually if something fails. The cost is that you must monitor the queue and wait between submissions. run_core_pipeline.sbatch is hands-off — submit once and come back — but the whole chain occupies a single large allocation, so the memory ceiling must accommodate the most memory-hungry step (pseudobulk DE and differential abundance, which each benefit from --mem=96G) even when earlier steps need far less. For a first run or for learning, run_core_pipeline.sbatch is simplest. For production runs on large or evolving datasets, the per-step approach wastes less shared compute and lets you re-run only the failing step.

10.5c Reading the output: confirming each step succeeded

Before launching the next script in the chain, confirm the preceding step genuinely completed. Two checks in combination:

  1. Look at the .out log for the Wrote ... line each script prints when it successfully writes its hand-off object to disk.
  2. Run seff <jobid> and verify State: COMPLETED — not OUT_OF_MEMORY or TIMEOUT, either of which would have left a truncated, unreadable .rds on disk.

A job killed mid-write by an OOM or timeout leaves a file that exists and appears to have a non-zero size, but the .rds is corrupt and the following readRDS() call will fail with a confusing error. The Wrote ... message in the log is the true signal that the object was fully written.

10.6 Running the pipeline: the two wrappers

You do not write a separate .sbatch file for each of the eight scripts. The workshop ships two wrappers that handle the SLURM boilerplate.

run_rscript.sbatch — one step at a time. This is a generic wrapper whose essential line is just Rscript "${SCRIPT}". You pass it the R script to run and override the job name (and, where needed, the memory) per submission. This is the right tool when you want fine control — running steps individually, watching each finish, and giving the memory-hungry steps more RAM:

mkdir -p logs   # one-time: SLURM will not create the log dir for you (Chapter 9)

sbatch --job-name=qc            run_rscript.sbatch 01_qc_preprocessing.R
sbatch --job-name=cluster       run_rscript.sbatch 02_dimreduction_clustering.R
sbatch --job-name=markers       run_rscript.sbatch 03_markers_annotation.R
sbatch --job-name=refannot      run_rscript.sbatch 04_reference_annotation.R
sbatch --job-name=integrate     run_rscript.sbatch 05_integration.R
sbatch --job-name=pb_de  --mem=96G run_rscript.sbatch 06_pseudobulk_de.R
sbatch --job-name=fa     --mem=32G run_rscript.sbatch 07_functional_analysis.R
sbatch --job-name=da     --mem=96G run_rscript.sbatch 08_differential_abundance.R

Because the steps are chained, you launch them one at a time and confirm each before starting the next. Concretely: look in the .out log for the Wrote ... line that the script prints when it has finished writing its object, and run seff <jobid> to confirm State: COMPLETED — not OUT_OF_MEMORY or TIMEOUT, either of which would have left a corrupt hand-off (§10.4).

run_core_pipeline.sbatch — the whole series in one allocation. When you simply want to reproduce the entire core analysis end-to-end, this wrapper runs all eight scripts in order inside a single SLURM allocation (roughly 12 hours, 96 GB):

STEPS=(01_qc_preprocessing.R 02_dimreduction_clustering.R … 08_differential_abundance.R)
for step in "${STEPS[@]}"; do
  echo "==== RUNNING ${step} ===="
  Rscript "${step}"
done

You submit it once with sbatch run_core_pipeline.sbatch (after mkdir -p logs and pre-caching ifnb), then confirm State: COMPLETED with seff. This is the simplest path to “run the whole thing.” The trade-off versus the one-at-a-time approach is granularity: a single fat allocation cannot right-size memory per step, so use the per-step wrapper when you need that control and the pipeline wrapper when you just want the full run.

10.7 The bonus scripts (13–17)

Beyond the core chain, five bonus scripts parallel the downstream bonus laptop tutorials. Unlike the core, these are largely standalone — they run after the core pipeline rather than in a fixed sequence with each other, and several bring their own datasets. Table 2 lists each script with its inputs and outputs.

Table 2: Bonus pipeline scripts (13–17). Unlike the core chain, these are largely standalone — each can run independently once its input data exist. The run_bonus_pipeline.sbatch wrapper runs all five in sequence but continues past failures so a missing dataset for one step does not abort the others.
Script Parallels notebook Reads Writes
13_wgcna.R Module 13 — WGCNA GSE152418_*_RawCounts.txt wgcna_*.csv
14_trajectory_cellcomm.R Module 14 — Trajectory & Cell–Cell Comm ifnb_annotated.rds ifnb_slingshot_pseudotime.csv
15_scatac_signac.R Module 15 — scATAC-seq (Signac) PBMC 10k scATAC files pbmc_atac_clustered.rds
16_spatial.R Module 16 — Spatial Transcriptomics stxBrain (SeuratData) brain_spatial_integrated.rds
17_fair_metadata.R Module 17 — FAIR & Metadata ifnb_annotated_final.rds ifnb_for_submission.h5ad

Scripts 13 (WGCNA) and 15 (scATAC) use their own datasets, downloaded into ../data/; 14 and 17 run on the annotated ifnb object the core pipeline produces; 16 uses the stxBrain spatial dataset. Because they do not chain on each other, you can run whichever bonus analyses interest you in any order, once the relevant input exists.

Just as the core eight have run_core_pipeline.sbatch, the bonus five have a parallel wrapper, run_bonus_pipeline.sbatch, that runs 13–17 in one allocation. Because they’re standalone (each needs its own dataset), that wrapper does not abort if one step fails — it logs the failure and continues to the next — so a missing dataset for, say, WGCNA won’t stop the spatial run. Submit it with sbatch run_bonus_pipeline.sbatch after mkdir -p logs.

One caveat, consistent with the notebooks: steps that are inherently interactive or Python-based — Monocle3 root selection, scVelo, CellChat plotting, the cellxgene-schema validator — cannot be fully automated in a batch script and are kept as commented sketches inside the scripts, exactly as they appear in the notebooks. The bonus scripts mirror the runnable core of each notebook, not its interactive parts. The one exception is RNA velocity (Module 14), which is packaged as a runnable Python script, 14_rna_velocity_scvelo.py, submitted via sbatch run_python.sbatch 14_rna_velocity_scvelo.py.

10.7b Per-step resource guide

Each core script has different CPU, memory, and wall-clock needs. Table 3 gives the practical defaults for the ifnb dataset and notes where to scale up for larger inputs. Use seff after each first run to verify your actual usage and trim accordingly.

Understanding why different steps have different needs makes the table more useful than memorizing it. Steps 01–04 are data-loading, normalization, clustering, and annotation — these are primarily memory-proportional-to-cell-count operations, and for ifnb (~24,000 cells) 32 GB is more than adequate. Step 05 (integration) is the first step that benefits from more cores: Harmony and anchor-based integration both parallelize, and the object being manipulated is the fully-annotated Seurat object across all conditions, making the memory footprint larger. Steps 06 and 08 are the memory peaks: pseudobulk DE builds one count matrix per cell type per donor (DESeq2’s model matrix can be large), and Milo’s neighborhood graph for differential abundance stores adjacency information across all cells. Both benefit from --mem=96G rather than the generic 32 GB default. Step 07 (functional analysis with fgsea) is computationally cheap by comparison — it reads a CSV, not a Seurat object — and 32 GB is generous.

Table 3: Per-step SLURM resource defaults for the ifnb dataset (~24,000 cells, 8 donors). Scale --mem proportionally for larger datasets; the pseudobulk DE (06) and differential abundance (08) steps spike highest. Always run seff <jobid> after first execution to right-size repeated runs.
Script --mem --cpus-per-task --time Notes
01_qc_preprocessing.R 32G 4 1:00:00 Fetches ifnb on first run; pre-cache on login node
02_dimreduction_clustering.R 32G 4 1:00:00 PCA is the slow step; parallelizes with future
03_markers_annotation.R 32G 4 2:00:00 FindAllMarkers is CPU-bound; more cores help
04_reference_annotation.R 32G 4 1:00:00 SingleR reference download needs internet (pre-cache)
05_integration.R 64G 8 4:00:00 Harmony/anchor integration is memory-intensive
06_pseudobulk_de.R 96G 4 2:00:00 DESeq2 in-memory model matrix drives memory spike
07_functional_analysis.R 32G 4 1:00:00 fgsea is fast; MSigDB download needs pre-caching
08_differential_abundance.R 96G 8 3:00:00 Milo neighborhood graph is both CPU and memory hungry

10.8 Adapting the pipeline to your own data

The whole point of learning this pipeline on ifnb is to run it on your data afterward, and the chained architecture makes that swap clean. The only place the pipeline touches the external world is script 01, which loads ifnb from Bioconductor. To run the analysis on your own counts:

Replace the loader in script 01. Swap the muscData::Kang18_8vs8() call for code that builds a Seurat object from your own counts matrix — for example, reading a Cell Ranger output directory or a 10x .h5 file — and carry your condition and donor information into the object’s metadata.

Keep the rest of the chain. Steps 02–08 read from the previous step’s .rds, so once script 01 writes a Seurat object in the expected shape, the downstream steps run unchanged. Mind the design: the DE and DA steps (05/06/08) depend on having a real condition and real replicates in the metadata. If your data lacks a condition or has no biological replication, those steps will not return meaningful results — the limitation is your experimental design, not the scripts. This is the same lesson §10.2 drew from why ifnb was chosen. Right-size for your scale: a dataset much larger than ifnb will need more memory and time; Table 3 gives the ifnb defaults as a starting point, and seff from Chapter 9 is how you tune each step once you have a first-run baseline.

Note⚙️ Key parameter — DATA_DIR (pipeline output directory)

Default: ../data/ (relative to the scripts directory). All eight pipeline scripts write their .rds and .csv hand-off files to the path in the DATA_DIR environment variable. If DATA_DIR is unset, it defaults to ../data/. For cluster runs, set it to a scratch location before submission: export DATA_DIR=/scratch/<PIRG>/<duckid>/ifnb/data. This keeps large intermediate Seurat objects (100–500 MB each) out of your home directory and on the high-throughput scratch filesystem. Because it is a shell variable, the same eight scripts run unchanged whether the data go to your project directory, scratch, or an external storage path — no edits to the scripts themselves are required.

Note⚙️ Key parameter — --mem per step (pipeline-specific defaults)

Default in run_rscript.sbatch: 32 GB. Most pipeline steps (01–04, 07) run comfortably at --mem=32G on the ifnb dataset. Steps 05 (integration) should be raised to --mem=64G; steps 06 and 08 (pseudobulk DE and differential abundance) should be raised to --mem=96G. These overrides are passed inline: sbatch --job-name=pb_de --mem=96G run_rscript.sbatch 06_pseudobulk_de.R. If you are adapting the pipeline to a dataset larger than ifnb (~24,000 cells), scale --mem proportionally — roughly linearly with cell count for steps 01–05, and more aggressively for 06 and 08 because DESeq2 and Milo both build in-memory data structures that scale with the number of cell types × donors × cells.

10.9 Where this leads

You now have the full picture: the cluster mechanics from Chapter 9, and in this chapter the architecture of the analysis that runs on it — a chain of .rds checkpoints, two wrappers to submit it, the pre-caching and DATA_DIR knobs that make it reproducible, and the ifnb design that gives the whole run real interferon biology to recover.

If you want to go upstream of script 01 — to start from raw FASTQs and generate the counts matrix yourself with Cell Ranger — see the bonus Cell Ranger modules and the Full Talapas RNA-seq run from raw FASTQs roadmap. Otherwise, you are ready to reproduce the entire single-cell workflow on Talapas and then point it at your own experiment.

Going further

References

1.
Kang, H. M. et al. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology 36, 89–94 (2018).
2.
Mölder, F. et al. Sustainable data analysis with Snakemake. F1000Research 10, 33 (2021).
3.
Di Tommaso, P. et al. Nextflow enables reproducible computational workflows. Nature Biotechnology 35, 316–319 (2017).
4.
Grüning, B. et al. Bioconda: Sustainable and comprehensive software distribution for the life sciences. Nature Methods 15, 475–476 (2018).
5.
Kurtzer, G. M., Sochat, V. & Bauer, M. W. Singularity: Scientific containers for mobility of compute. PLoS ONE 12, e0177459 (2017).
Back to top