3. Resource Management

Why Resources Matter

FastQC and fastp are fast and lightweight. You can run many of them simultaneously without issue. STAR alignment is different: it loads the entire genome index into memory and saturates CPU cores efficiently. If Snakemake launches three STAR jobs at once on a machine with 8 cores and 32 GB RAM, each job would compete for the same resources and run slower than if they were scheduled intelligently.

The threads: and resources: directives let you declare what each rule needs. Snakemake uses these declarations to schedule jobs without exceeding your system's capacity.

Adding `rule star_align`

Add the following two configuration variables at the top of your Snakefile, directly below the READS line:

GENOME_DIR = "genome/star_index"
GTF        = "genome/Saccharomyces_cerevisiae.R64-1-1.gtf"

Then add rule star_align below rule fastp:

rule star_align:
    input:
        r1=   "results/trimmed/{sample}_R1.fastq.gz",
        r2=   "results/trimmed/{sample}_R2.fastq.gz",
        index=GENOME_DIR
    output:
        bam="results/aligned/{sample}.Aligned.sortedByCoord.out.bam",
        log="results/aligned/{sample}.Log.final.out"
    threads: 8
    resources:
        mem_mb=16000
    shell:
        """
        mkdir -p results/aligned
        STAR \
            --runThreadN {threads} \
            --genomeDir {input.index} \
            --readFilesIn {input.r1} {input.r2} \
            --readFilesCommand zcat \
            --outSAMtype BAM SortedByCoordinate \
            --outFileNamePrefix results/aligned/{wildcards.sample}. \
            --outSAMattributes NH HI AS NM MD
        """

Key details:

threads: 8 — Snakemake passes this value to {threads} in the shell command, so STAR receives it via --runThreadN. When you run snakemake --cores 16, Snakemake knows each STAR job needs 8 cores and will run at most 2 in parallel.
resources: mem_mb=16000 — declares 16 GB of memory. For S. cerevisiae, STAR uses far less, but this pattern is critical when working with human genomes (~30 GB).
{wildcards.sample} — when you need the wildcard's value inside the shell command (not in an input/output path), access it via {wildcards.sample}. This is necessary here because STAR constructs output filenames from --outFileNamePrefix.

Adding `rule featurecounts`

featureCounts takes all three sorted BAM files at once and produces a single count matrix. This is a many-to-one aggregation step. It uses expand() in its input: rather than a wildcard, because it is not run per-sample.

rule featurecounts:
    input:
        bams=expand("results/aligned/{sample}.Aligned.sortedByCoord.out.bam",
                    sample=SAMPLES),
        gtf=GTF
    output:
        counts="results/counts/all_samples.txt"
    threads: 4
    shell:
        """
        mkdir -p results/counts
        featureCounts \
            -T {threads} \
            -p \
            -a {input.gtf} \
            -o {output.counts} \
            {input.bams}
        """

The -p flag tells featureCounts that reads are paired-end. The {input.bams} placeholder expands to all three BAM paths, space-separated.

Adding `rule multiqc`

MultiQC is the final step. It scans the results/ directory for log files generated by FastQC, fastp, STAR, and featureCounts, and compiles them into a single interactive HTML report.

rule multiqc:
    input:
        expand("results/fastqc/{sample}_{read}_fastqc.zip",
               sample=SAMPLES, read=READS),
        expand("results/trimmed/{sample}_fastp.json", sample=SAMPLES),
        expand("results/aligned/{sample}.Log.final.out", sample=SAMPLES),
        "results/counts/all_samples.txt"
    output:
        "results/multiqc/multiqc_report.html"
    shell:
        "multiqc results/ --outdir results/multiqc/ --force"

The input: list here serves an important purpose: it forces Snakemake to complete all upstream steps before launching MultiQC. MultiQC itself does not need these files passed on the command line. It finds them by scanning the directory, but listing them as input: ensures the dependency chain is enforced.

Update `rule all`

Simplify rule all to a single final target. Snakemake will trace backward through the full dependency chain automatically:

rule all:
    input:
        "results/multiqc/multiqc_report.html"

Your Complete Snakefile So Far

SAMPLES    = ["wt_rep1", "wt_rep2", "wt_rep3", "mut_rep1", "mut_rep2", "mut_rep3"]
READS      = ["R1", "R2"]
GENOME_DIR = "genome/star_index"
GTF        = "genome/Saccharomyces_cerevisiae.R64-1-1.gtf"

rule all:
    input:
        "results/multiqc/multiqc_report.html"

rule fastqc:
    input:
        "data/{sample}_{read}.fastq.gz"
    output:
        html="results/fastqc/{sample}_{read}_fastqc.html",
        zip= "results/fastqc/{sample}_{read}_fastqc.zip"
    shell:
        "fastqc {input} --outdir results/fastqc/"

rule fastp:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz"
    output:
        r1=  "results/trimmed/{sample}_R1.fastq.gz",
        r2=  "results/trimmed/{sample}_R2.fastq.gz",
        json="results/trimmed/{sample}_fastp.json",
        html="results/trimmed/{sample}_fastp.html"
    shell:
        """
        fastp \
            --in1 {input.r1}  --in2 {input.r2} \
            --out1 {output.r1} --out2 {output.r2} \
            --json {output.json} --html {output.html}
        """

rule star_align:
    input:
        r1=   "results/trimmed/{sample}_R1.fastq.gz",
        r2=   "results/trimmed/{sample}_R2.fastq.gz",
        index=GENOME_DIR
    output:
        bam="results/aligned/{sample}.Aligned.sortedByCoord.out.bam",
        log="results/aligned/{sample}.Log.final.out"
    threads: 8
    resources:
        mem_mb=16000
    shell:
        """
        mkdir -p results/aligned
        STAR \
            --runThreadN {threads} \
            --genomeDir {input.index} \
            --readFilesIn {input.r1} {input.r2} \
            --readFilesCommand zcat \
            --outSAMtype BAM SortedByCoordinate \
            --outFileNamePrefix results/aligned/{wildcards.sample}. \
            --outSAMattributes NH HI AS NM MD
        """

rule featurecounts:
    input:
        bams=expand("results/aligned/{sample}.Aligned.sortedByCoord.out.bam",
                    sample=SAMPLES),
        gtf=GTF
    output:
        counts="results/counts/all_samples.txt"
    threads: 4
    shell:
        """
        mkdir -p results/counts
        featureCounts \
            -T {threads} \
            -p \
            -a {input.gtf} \
            -o {output.counts} \
            {input.bams}
        """

rule multiqc:
    input:
        expand("results/fastqc/{sample}_{read}_fastqc.zip",
               sample=SAMPLES, read=READS),
        expand("results/trimmed/{sample}_fastp.json", sample=SAMPLES),
        expand("results/aligned/{sample}.Log.final.out", sample=SAMPLES),
        "results/counts/all_samples.txt"
    output:
        "results/multiqc/multiqc_report.html"
    shell:
        "multiqc results/ --outdir results/multiqc/ --force"

Running the Full Pipeline

Do a dry-run first to verify that all 14 planned jobs look correct (6 FastQC + 3 fastp + 3 STAR + 1 featureCounts + 1 MultiQC):

snakemake -n

Then run:

snakemake --cores 16

With 16 cores and threads: 8 on rule star_align, Snakemake can align two samples in parallel while simultaneously running FastQC and fastp jobs on the remaining cores.

Use all available cores

On a machine where you have exclusive access, skip the mental arithmetic:

snakemake --cores all

Snakemake detects the number of available CPUs and uses them all, respecting each rule's threads: declaration.

Where we are

You now have a complete, parallelised 5-step RNA-seq pipeline in ~40 lines. The threads: and resources: directives do more than just control local parallelism. In the next section they become SLURM job specifications automatically. Continue to HPC Integration.

3. Resource Management

Why Resources Matter

Adding rule star_align

Adding rule featurecounts

Adding rule multiqc

Update rule all