# Building a Robust RNA-seq QC Pipeline with Nextflow: A Template By the end of this guide you will have a working Nextflow RNA-seq QC pipeline that takes a samplesheet of gzipped FASTQ files and returns raw-read quality reports, adapter and quality trimming, post-trim reports, optional alignment statistics from STAR, and a single aggregated MultiQC report for the whole project. You will be able to run it on a laptop or a cluster, resume it after a failure without redoing finished work, read every QC metric it produces, and fix the errors that show up most often in real runs. Have these on hand before you start: - A Unix shell (Linux, macOS, or Windows Subsystem for Linux) - Java 17 or later - Miniforge or Miniconda, which supply the software environments - At least two paired-end RNA-seq samples in gzipped [FASTQ format](/knowledge/bioinformatics/fastq-file-format-phred-scores-qc) - For the optional alignment stage, a reference genome FASTA, a GTF annotation file, and a machine with roughly 30 GB of RAM FASTQ files, the plain-text format that stores sequencing reads as four lines per read (header, sequence, separator, per-base quality string), are the only required input. Everything in this pipeline starts from them. ## Why Nextflow for RNA-seq QC Nextflow is a workflow manager. You describe each analysis step as a process, connect processes with channels, and Nextflow handles execution order, parallelism, caching, and retries. A Nextflow RNA-seq QC pipeline gives you three things a shell script cannot: 1. **Accountability.** Every task records its inputs, command, and outputs in a content-addressed work directory, so a crashed run resumes where it stopped. 2. **Pinned software.** Conda environments or containers lock tool versions, so the same command produces the same environment on a laptop and on a cluster. 3. **Portability.** The same code runs under a local executor or a scheduler such as Slurm without changing the analysis logic. Robustness here means the pipeline survives failures, records what it did, and reproduces the same result when rerun. The published literature backs this design. The nf-rnaSeqCount project built an RNA-seq pipeline on Nextflow with Singularity containerization specifically to create a robust, portable, and reproducible workflow that automates analysis from QC through alignment, and its authors note that such pipelines reduce analysis time while increasing the accuracy and reproducibility of transcriptome research [1]. The same paper frames the motivation directly: large omics datasets bring computational challenges and potential pitfalls, and improvised analysis is where those pitfalls multiply [1]. ## The RNA-seq errors your QC pipeline must catch QC exists because RNA-seq errors cluster into recognizable patterns. Catch them at the read and library level, before they reach differential expression where they masquerade as biology. **Instrument-level problems** - **Base quality decay.** Quality scores typically drop toward the ends of reads. A Phred score Q is defined as Q = -10 × log10(P), where P is the probability that the base call is wrong. Q30 means 99.9 percent accuracy, one error per thousand bases. Q20 means one error per hundred. A library whose tail quality sits below Q20 contributes noisy bases to every downstream count. - **Overclustering.** Too dense a loading on the flow cell degrades quality across whole runs, visible as uniformly low scores rather than a clean tail. - **Index misassignment.** Reads carry the wrong sample barcode at low but nonzero rates, which matters most in low-input experiments. **Library-level problems** - **Adapter read-through.** When the inserted DNA fragment is shorter than the read length, the sequencer reads into the adapter sequence. An adapter is the short synthetic DNA oligo ligated to fragment ends during library prep. Untrimmed adapters poison alignment. - **Duplication.** PCR amplification and optical duplicates inflate read counts for the same original fragment. High duplication in a low-input library means the effective depth is far lower than the raw read count suggests. - **rRNA carryover.** [Ribosomal RNA](/knowledge/bioinformatics/ribosomal-rna-structure-taxonomic-profiling) reads indicate failed poly(A) selection or ribodepletion, and they waste [sequencing depth](/knowledge/diagnostics/molecular/how-to-calculate-sequencing-depth-and-coverage-for-your-ngs-run). - **Positional bias.** Degraded RNA shifts coverage toward transcript ends. Some 3'-end chemistries produce strong 3' bias by design, which is expected behavior rather than damage. - **Strandedness mismatch.** Strandedness describes which reads correspond to which RNA strand, determined by the library kit. Counting features with the wrong strandedness setting produces a count matrix dominated by antisense signal. **Sample-level problems** - **Swaps and outliers.** Mislabelled files, cross-contamination, and failed libraries appear as samples that refuse to cluster with their replicates. - **Wrong reference.** Aligning against the wrong genome build silently lowers mapping rates and distorts gene-level counts. QC detects these patterns. It does not repair them. Trimming removes adapters and low-quality tails, but no tool recovers a degraded library or fixes a swap. ## What the template does ![Flow of the RNA-seq QC pipeline: samplesheet parsing, raw FastQC, fastp trimming, trimmed FastQC, optional STAR alignment and MultiQC aggregation.](/article-images/illustrations/nextflow-rna-seq-qc-pipeline-qc-stages.webp) The pipeline runs seven stages in order: 1. Parse the samplesheet into a channel of sample metadata plus read pairs. 2. Run FastQC on raw reads. 3. Trim adapters and low-quality tails with fastp. 4. Run FastQC on trimmed reads. 5. Optionally align with STAR and summarize mapping with samtools flagstat. 6. Aggregate every report with MultiQC, a tool that parses outputs from dozens of bioinformatics programs into one HTML page. 7. Publish results and record software versions per sample. ## Prerequisites and installation Check Java first: ``` java -version ``` Expected output (yours will differ slightly): ``` openjdk version "17.0.12" 2024-07-16 OpenJDK Runtime Environment Temurin-17.0.12+7 ``` Install Nextflow, the launcher is a single self-updating script: ``` curl -s https://get.nextflow.io | bash sudo mv nextflow /usr/local/bin/ nextflow -version ``` Expected output: ``` N E X T F L O W version 24.10.0 ``` The full line includes a build number that varies by release. Any 24.x or later version runs this template. For software environments, install Miniforge from the conda-forge GitHub project (github.com/conda-forge/miniforge). Conda is a package manager that installs pinned tool versions into isolated environments, and recent Miniforge installers include Mamba, a faster dependency solver. This template sets `conda.useMamba = true` so Nextflow builds each process environment in seconds rather than minutes. ## Project layout and samplesheet format A Nextflow pipeline has a small, predictable format: a main script, a configuration file, and a folder of assets. Create the skeleton: ``` mkdir -p rna-seq-qc/assets cd rna-seq-qc touch main.nf nextflow.config samplesheet.csv ``` The files that matter: - `main.nf` holds the workflow and all processes. - `nextflow.config` holds parameters, resources, and profiles. - `samplesheet.csv` lists your input files. - `assets/` holds reference files and test data. - `work/` and `results/` appear at runtime. `work/` stores each task's private sandbox and can grow to hundreds of gigabytes. `results/` holds only what you publish. Run `nextflow clean -f` when you no longer need the cache and want the disk space back. The samplesheet is a plain CSV. The header names are fixed by the template: ``` sample,fastq_1,fastq_2 CTRL_rep1,reads/CTRL_rep1_R1.fastq.gz,reads/CTRL_rep1_R2.fastq.gz CTRL_rep2,reads/CTRL_rep2_R1.fastq.gz,reads/CTRL_rep2_R2.fastq.gz TREAT_rep1,reads/TREAT_rep1_R1.fastq.gz,reads/TREAT_rep1_R2.fastq.gz TREAT_rep2,reads/TREAT_rep2_R1.fastq.gz,reads/TREAT_rep2_R2.fastq.gz ``` Rules that prevent most input failures: keep the header names exactly as shown, write paths relative to the directory where you launch the pipeline, list R1 before R2, gzip every file, make sample IDs unique, and save the file with Unix (LF) line endings. ## Building the pipeline in main.nf ### The workflow block and channels A channel is a stream of data flowing between processes. The workflow block builds the input channel and wires the processes together: ```groovy #!/usr/bin/env nextflow nextflow.enable.dsl = 2 workflow { reads_ch = Channel.fromPath(params.input, checkIfExists: true) .splitCsv(header: true) .map { row -> def meta = [id: row.sample, stage: 'raw'] tuple(meta, [file(row.fastq_1, checkIfExists: true), file(row.fastq_2, checkIfExists: true)]) } FASTQC(reads_ch) FASTP(reads_ch) FASTQC(FASTP.out.reads) if (params.star_index) { STAR_ALIGN(reads_ch, Channel.value(file(params.star_index, checkIfExists: true))) SAMTOOLS_FLAGSTAT(STAR_ALIGN.out.bam) } qc_inputs = FASTQC.out.zips.mix(FASTP.out.json) if (params.star_index) { qc_inputs = qc_inputs.mix(STAR_ALIGN.out.log) } MULTIQC(qc_inputs.collect()) } ``` Key points: - `nextflow.enable.dsl = 2` selects DSL2, Nextflow's current scripting model that organizes code into processes and workflows. Recent Nextflow versions default to DSL2, and the line keeps the script compatible with older installs. - The `meta` map (sample ID and stage) travels with the reads through every process, which is how reports stay attached to the right sample. - `checkIfExists: true` makes the pipeline fail within seconds on a typo in the samplesheet instead of failing deep inside a task. - FASTQC is called twice, once on raw reads and once on trimmed reads. Both invocations publish into the same output channel, so `FASTQC.out.zips` contains all of them. - The alignment stage runs only when `params.star_index` is set. - `.collect()` gathers many channel items into one list, which is what MultiQC needs to see everything at once. If Nextflow prints a warning that DSL 1 is deprecated, you are running an old version or the enable line is missing. Add the line or upgrade. ### FASTQC ```groovy process FASTQC { tag "${meta.id}_${meta.stage}" label 'process_low' conda 'bioconda::fastqc=0.11.9' publishDir "${params.outdir}/fastqc", mode: 'copy' input: tuple val(meta), path(reads) output: tuple val(meta), path("fastqc/*_fastqc.zip"), emit: zips path "fastqc/*_fastqc.html", emit: htmls script: """ fastqc --outdir fastqc --threads ${task.cpus} ${reads} fastqc --version > ${meta.id}_${meta.stage}_version.txt """ } ``` Each directive has a job: - `tag` labels the task in the console with the sample and stage, so you can see `CTRL_rep1_raw` versus `CTRL_rep1_trimmed` at a glance. - `label` attaches a resource block from the config (2 CPUs, 4 GB here). - `conda` pins the exact tool version. Nextflow builds that environment automatically under the conda profile. - `publishDir` copies outputs into `results/fastqc` while the raw files stay in `work/`. - The tuple input carries the metadata map plus both read files, and `${reads}` expands to both file paths in the script. - The version line records the tool version per sample, a habit that pays off when you need to prove what ran. ### FASTP fastp performs adapter removal and quality trimming in one fast pass and writes its own per-sample report: ```groovy process FASTP { tag "${meta.id}" label 'process_low' conda 'bioconda::fastp=0.23.2' publishDir "${params.outdir}/fastp", mode: 'copy' input: tuple val(meta), path(reads) output: tuple val(trimmed_meta), path("${meta.id}_trimmed_*.fastq.gz"), emit: reads path "fastp_${meta.id}.html", emit: html path "fastp_${meta.id}.json", emit: json script: def trimmed_meta = meta + [stage: 'trimmed'] """ fastp \ --in1 ${reads} \ --in2 ${reads[1]} \ --out1 ${meta.id}_trimmed_1.fastq.gz \ --out2 ${meta.id}_trimmed_2.fastq.gz \ --detect_adapter_for_pe \ --qualified_quality_phred 20 \ --length_required 36 \ --thread ${task.cpus} \ --html fastp_${meta.id}.html \ --json fastp_${meta.id}.json \ 2> fastp_${meta.id}.log """ } ``` Parameter choices: - `detect_adapter_for_pe` auto-detects the adapter sequences from the paired data, which avoids hard-coding the wrong kit adapter. - `qualified_quality_phred 20` sets Q20 as the bar for a high-quality base, and fastp's default rule discards reads where more than 40 percent of bases fall below it. - `length_required 36` discards trimmed reads shorter than 36 nucleotides, which removes adapter remnants that would otherwise fail alignment. - Report and JSON files carry the sample ID in their names. Without that, every sample writes `fastp.html` into the same folder and the outputs collide. - `trimmed_meta = meta + [stage: 'trimmed']` renames the stage so the second FastQC pass publishes and reports under the trimmed label. ### Optional alignment QC with STAR Alignment turns read-level QC into library-level QC. STAR needs a genome index, a precomputed data structure that lets the aligner find read positions quickly. Build it once with the genome FASTA and a GTF file, the standard annotation format listing genes, transcripts, and exon coordinates: ``` STAR --runMode genomeGenerate \ --genomeDir star_index \ --genomeFastaFiles genome.fa \ --sjdbGTFfile genes.gtf \ --genomeSAindexNbases 14 \ --runThreadN 8 ``` The suffix array parameter follows the documented STAR formula, min(14, log2(genome length)/2 - 1). A 3 Gb mammalian genome gives 14, a 100 Mb genome gives about 12. Building and loading a 3 Gb index needs on the order of 30 GB of RAM, which is why the alignment stage carries the high resource label. ```groovy process STAR_ALIGN { tag "${meta.id}" label 'process_high' conda 'bioconda::star=2.7.10a' publishDir "${params.outdir}/star", mode: 'copy', pattern: '*Log.final.out' input: tuple val(meta), path(reads) path index output: tuple val(meta), path("${meta.id}Aligned.sortedByCoord.out.bam"), emit: bam path "${meta.id}Log.final.out", emit: log script: """ STAR \ --genomeDir ${index} \ --readFilesIn ${reads} ${reads[1]} \ --readFilesCommand zcat \ --runThreadN ${task.cpus} \ --outFileNamePrefix ${meta.id} \ --outSAMtype BAM SortedByCoordinate \ --outBAMsortingThreadN ${task.cpus} \ --quantMode GeneCounts """ } process SAMTOOLS_FLAGSTAT { tag "${meta.id}" label 'process_low' conda 'bioconda::samtools=1.17' publishDir "${params.outdir}/flagstat", mode: 'copy' input: tuple val(meta), path(bam) output: path "${meta.id}.flagstat", emit: flagstat script: """ samtools flagstat ${bam} > ${meta.id}.flagstat """ } ``` Notes on the choices: - `readFilesCommand zcat` tells STAR the inputs are gzipped. Forgetting it is one of the most common RNA-seq errors in hand-written STAR commands. - The `publishDir` pattern publishes only the small log file. The sorted BAM, which can exceed 10 GB per sample, stays in `work/` where downstream processes can still use it. - `quantMode GeneCounts` produces `ReadsPerGene.out.tab`, a per-gene count table. Count matrices are the standard handoff to differential expression: nf-rnaSeqCount was built for exactly this deliverable, mapping raw RNA-seq reads to a reference genome, quantifying genomic features, and passing the resulting count matrix to tools such as DESeq2 and edgeR [1]. STAR's GeneCounts output is a quick version of that matrix. For strand-aware counting you control explicitly, run featureCounts from the subread package on the BAM instead. ### MULTIQC ```groovy process MULTIQC { label 'process_low' conda 'bioconda::multiqc=1.14' publishDir "${params.outdir}/multiqc", mode: 'copy' input: path reports output: path "multiqc_report.html", emit: report path "multiqc_data", emit: data script: """ multiqc . -n multiqc_report.html -f """ } ``` MultiQC scans the staged files recursively and recognizes FastQC zip archives, fastp JSON reports, and STAR log files natively. Because it only aggregates what you feed it, a report missing a tool usually means that tool's output never reached the process input. ## nextflow.config explained ```groovy params { input = 'samplesheet.csv' outdir = 'results' star_index = null } manifest { name = 'rna-seq-qc-template' description = 'RNA-seq QC pipeline template built with Nextflow DSL2' version = '1.0.0' mainScript = 'main.nf' } process { errorStrategy = 'retry' maxRetries = 2 withLabel: process_low { cpus = 2 memory = 4.GB time = '2 h' } withLabel: process_high { cpus = 8 memory = 32.GB time = '8 h' } } profiles { conda { conda.enabled = true conda.useMamba = true } slurm { process.executor = 'slurm' conda.enabled = true conda.useMamba = true } test { params.input = "$projectDir/assets/test_samplesheet.csv" } } ``` - `errorStrategy = 'retry'` with `maxRetries = 2` reruns a failed task twice before the pipeline stops, which absorbs transient cluster faults. - The two labels separate light tasks (FastQC, fastp, MultiQC) from the memory-hungry alignment stage. - Profiles switch behavior with one flag: `-profile conda` on a laptop, `-profile slurm` on a cluster, `-profile test` for a built-in test samplesheet. Key options and formulas in one place: | Parameter or flag | Set where | What it controls | Sensible value | |, - |, - |, - |, - | | `input` | CLI or config | Path to the samplesheet CSV | `samplesheet.csv` | | `outdir` | CLI or config | Published results folder | `results` | | `star_index` | config | Directory with a built STAR index | `null` keeps alignment off | | `qualified_quality_phred` | fastp | Base quality counted as high quality | 20 (fastp default is 15) | | `length_required` | fastp | Minimum read length after trimming | 36 (fastp default is 15) | | `genomeSAindexNbases` | STAR index build | Suffix array bucket size | min(14, log2(genome length)/2 - 1) | | `cpus`, `memory` | process labels | Resources per task | 2 CPUs, 4 GB (low), 8 CPUs, 32 GB (high) | ## Running the Nextflow RNA-seq QC pipeline From the project directory: ``` nextflow run main.nf --input samplesheet.csv --outdir results -profile conda ``` The first launch builds the Conda environments, then the tasks run. Expected console output: ``` N E X T F L O W ~ version 24.10.0 Launching `main.nf` [soggy_curie] DSL2 - revision: 8f3a1b2c4d executor > local (13) [1a/2b3c4d] process > FASTQC (TREAT_rep2_trimmed) [100%] 8 of 8 ✔ [5e/6f7a8b] process > FASTP (TREAT_rep2) [100%] 4 of 4 ✔ [9c/0d1e2f] process > MULTIQC [100%] 1 of 1 ✔ Completed at 11-NOV-2026 14:32:11 Duration: 6m 12s ``` The session name and revision hash vary per run. FASTQC shows 8 of 8 because the process ran on both raw and trimmed reads for four samples. Published results look like this: ``` results/ fastqc/ CTRL_rep1_R1_fastqc.html CTRL_rep1_R1_fastqc.zip CTRL_rep1_trimmed_1_fastqc.html ... fastp/ fastp_CTRL_rep1.html fastp_CTRL_rep1.json ... multiqc/ multiqc_report.html multiqc_data/ multiqc_general_stats.txt ... ``` Rerun with `-resume` after changing anything (a parameter, a sample, one process script): ``` nextflow run main.nf --input samplesheet.csv --outdir results -profile conda -resume ``` Nextflow hashes each task's code, inputs, and software environment. Anything unchanged prints as cached and finishes in seconds. Only genuinely affected tasks rerun. For run accounting, add trace and timeline output: ``` nextflow run main.nf --input samplesheet.csv -profile conda -with-trace trace.tsv -with-timeline timeline.html ``` `trace.tsv` lists every task with its real CPU time, peak memory, and status, which is how you right-size the resource labels later. ## Reading the QC results Open `multiqc_report.html` first. The General Statistics table puts one row per sample with total sequences, duplication, and per-tool metrics side by side. Then drill into the per-sample reports. From a STAR log, expect a summary like: ``` Uniquely mapped reads % | 87.4 Unmapped reads % | 4.1 Multiple loci % | 8.5 Mismatch rate per base, % | 0.2 ``` From flagstat, a shortened example: ``` 2000000 + 0 in total (QC-passed reads + QC-failed reads) 1894000 + 0 mapped (94.70% : N/A) 1800000 + 0 properly paired (91.84% : N/A) ``` The two reports answer different questions. STAR splits mapped reads into unique versus multi-mapping, while flagstat counts everything that aligned, so flagstat's mapped percentage will sit above STAR's uniquely mapped percentage for the same BAM. For strandedness, run RSeQC (installable from Bioconda) on a sorted BAM with a BED12 gene model built from the same GTF: ``` infer_experiment.py -r genes.bed -i CTRL_rep1Aligned.sortedByCoord.out.bam ``` Example output and interpretation: ``` Fraction of reads explained by "1++,1--,2+-,2--": 0.03 Fraction of reads explained by "1+-,1-+,2++,2--": 0.96 ``` A fraction near 1 on one line means a stranded library (this example is reverse-stranded, typical of dUTP-type kits). Both lines near 0.5 means unstranded. Set your counting tool accordingly, because featureCounts uses `-s 0` for unstranded, `-s 1` or `-s 2` for the two stranded orientations. The metrics worth watching, in one table: | Metric | Where to check | Healthy pattern | Common causes when abnormal | |, - |, - |, - |, - | | Per-base quality | FastQC per-base plot | Most positions at or above Q30 | Chemistry decay, overclustering | | Adapter content | FastQC adapter plot, fastp report | Near zero after trimming | Short inserts, wrong adapter set | | Duplication | FastQC, fastp, MultiQC | Low to moderate | Low input, over-sequencing | | Uniquely mapped | STAR log | Most reads map uniquely in whole-transcriptome libraries | Wrong reference, contamination, over-trimming | | rRNA fraction | Alignment statistics | Low for poly(A) or ribodepleted libraries | Selection or depletion failure | | Gene body coverage | RSeQC geneBody_coverage | Roughly uniform (3'-end libraries excepted) | Degradation, or expected 3' chemistry | | Strandedness | RSeQC infer_experiment | Matches the kit | Wrong counting settings | | Replicate clustering | MultiQC PCA or heatmap | Replicates group by condition | Swaps, outliers, batch effects | ## Adapting the template to your library type **Stranded libraries.** Once `infer_experiment.py` confirms orientation, use featureCounts with the matching `-s` setting rather than relying on defaults. A strandedness mismatch is silent at the pipeline level and catastrophic at the differential expression level. **3'-end and oligo(dT) libraries.** These chemistries concentrate reads near transcript 3' ends by design, so geneBody_coverage plots will show strong 3' skew that is not degradation. Specialized pipelines handle this space well. PeakPrime, for example, is a Nextflow pipeline for 3'-end RNA-seq that identifies high-coverage regions with MACS2, designs strand-appropriate primers with Primer3, screens specificity with Bowtie2, and emits QC summaries and plots as first-class outputs [2]. Its validation numbers show what targeted enrichment can and cannot do: targeted libraries preserved gene-level expression correlations with random-primed 3'-end RNA-seq (Pearson r ≈ 0.88 across conditions) while raising the fraction of reads assigned to target genes from under 1 percent to roughly 25 to 42 percent and improving detection of the lowest-abundance transcripts [2]. The QC lesson generalizes: judge a library against the expectations of its chemistry, not against whole-transcriptome rules. **Low-input and small-RNA libraries.** Expect high duplication as a property of the protocol rather than a defect, and treat duplicate removal as mandatory rather than optional. ## Bioinformatics troubleshooting: exact errors and fixes ![Decision tree for common Nextflow RNA-seq QC failures: Java version, FASTQ file checks, MultiQC results and resume cache.](/article-images/illustrations/nextflow-rna-seq-qc-pipeline-troubleshoot-path.webp) ### UnsupportedClassVersionError at launch The JVM throws `UnsupportedClassVersionError` naming a class file version when you run `nextflow`. Your Java is older than the release requires. Install Temurin 17 or later, confirm with `java -version`, and if several JDKs coexist, point `JAVA_HOME` at the new one. ### Input file name collision Nextflow aborts with an input file name collision and names the duplicated file. Two inputs staged into one task share a filename, usually because sample IDs repeat in the samplesheet or two channels deliver same-named files. Make output names unique per sample (as the template does with `${meta.id}` prefixes), or namespace one input with `stageAs`: ```groovy input: tuple val(meta), path(reads, stageAs: 'reads/*') ``` ### Process requirement exceeds available memory Nextflow prints "Process requirement exceeds available memory" before the task starts. The label requests more than the machine or queue provides. Lower the label memory in `nextflow.config`, or route the heavy stage to a large-memory queue. The 32 GB label on STAR exists because the mammalian index needs it, so on a laptop either skip the alignment stage or run it on a beefier node. ### MultiQC finds no analysis results MultiQC finishes and reports that it found no analysis results. It scanned only files it does not recognize. Look inside the MULTIQC task directory under `work/` (find the hash with `nextflow log`) and confirm the FastQC zips and fastp JSON files actually landed there. If they did, run `multiqc` manually on that folder to see the parse warnings. ### fastp or gzip fails on a FASTQ file A task fails and the log shows gzip reporting invalid compressed data. The file is truncated, or a download saved an HTML error page with a `.fastq.gz` name. Test with `gzip -t reads/sample_R1.fastq.gz` and `file reads/sample_R1.fastq.gz` (a healthy file reports gzip data). Re-transfer the file and verify its checksum against the source. ### STAR exits with a fatal error on a small genome STAR exits with a fatal error telling you to decrease `genomeSAindexNbases`, or the process dies while loading the index. The first case means the suffix array parameter is too large for the genome size, so recompute it with the formula min(14, log2(genome length)/2 - 1). The second case means the node lacks the RAM to hold the index, so move the stage to a large-memory machine. ### The samplesheet parses but files are not found Tasks fail immediately with file not found, or the input channel comes up empty. The usual culprits: Windows CRLF line endings leave a carriage return at the end of every path, the header names do not match `sample`, `fastq_1`, `fastq_2`, or the paths are relative to the wrong directory. Convert the file with `dos2unix samplesheet.csv`, fix the header, and launch from the directory that contains `reads/`. ### Resume does not skip finished tasks `-resume` reruns everything. Nextflow invalidates the cache when the process script text, the input files, or the software directives change, and it never reuses a task whose files you edited inside `work/`. Keep process scripts byte-identical between runs, never touch `work/` by hand, and check `nextflow log` to compare run hashes. Remember that changing a resource label also invalidates every task carrying that label. ### Conda environment creation fails or crawls Set `conda.useMamba = true`, pin exact package versions in the `conda` directives, and delete a corrupted `work/conda` directory if a build died halfway. On shared systems, a warm Conda package cache makes the second pipeline launch nearly instant. ## Verifying that the pipeline is correct A completed run is not a validated run. Check these before trusting the numbers: 1. **Sample completeness.** MultiQC's General Statistics table should list every sample from the samplesheet. A missing row means a task failed or an output was never staged. 2. **Cross-report agreement.** For each sample, flagstat's mapped percentage should sit modestly above STAR's uniquely mapped percentage. A huge gap means excessive multi-mapping. 3. **Count sanity.** `ReadsPerGene.out.tab` from STAR has four columns: gene ID, unstranded counts, forward-stranded counts, reverse-stranded counts. Confirm the row count roughly matches your annotation's gene count, and confirm a known housekeeping gene carries plausible nonzero counts. 4. **Caching behavior.** Run twice with `-resume`. The second run should print cached status for every task. If it does not, something is nondeterministic, usually a timestamp written into an output name. 5. **Version records.** The per-sample version files in `results/fastqc/` should show the versions you pinned. 6. **Resource reality.** Open `trace.tsv` from a `-with-trace` run and compare requested versus used memory. Labels that request 32 GB for tasks using 3 GB waste cluster hours. Sensible next steps once the template passes these checks: add a featureCounts process for strand-aware counts, wrap the RSeQC commands into processes so strandedness and gene body coverage land in MultiQC automatically, put the pipeline under Git with a test profile in continuous integration, and graduate production studies to nf-core/rnaseq, which implements the same architecture at community-maintained scale. ## Common Mistakes and Limitations 1. **Treating FastQC warnings as failures.** FastQC's flag thresholds were tuned for small-genome sequencing. RNA-seq libraries legitimately trip the duplication and per-sequence quality warnings. Read the plots, not the pass or fail banner. 2. **Over-trimming.** Aggressive quality and length cutoffs delete short but mappable reads. Compare raw versus trimmed mapping rates for a few samples before locking in parameters. 3. **Ignoring strandedness until after differential expression.** By then the damage is a garbage count matrix. Run `infer_experiment.py` on one sample per library prep, early. 4. **Deleting `work/` before validation.** The work directory is your audit trail and your resume cache. Clean it only after you have published, verified, and archived what you need. 5. **Equating exit code zero with correctness.** Tools that complete successfully can still be fed the wrong reference, the wrong strandedness, or the wrong samplesheet row. The checks in the verification section exist for this reason. 6. **Skipping version records.** Six months later, nobody remembers which fastp version produced a given report. The per-sample version files cost one line of code. 7. **Feeding MultiQC only flattering inputs.** A report that omits the raw-read FastQC because it looks bad hides the trim benefit and the original problem alike. Aggregate everything. 8. **Scope limits.** This template automates routine checks. It does not replace judgment about a specific dataset, which still needs an experienced analyst, and production-scale studies should adopt a community-maintained pipeline such as nf-core/rnaseq rather than growing a homegrown template indefinitely. ## Frequently Asked Questions ### How much compute does an RNA-seq QC pipeline need? A laptop with 8 GB of RAM handles the FastQC, fastp, and MultiQC stages for a small project. The optional STAR stage is the demanding part, needing on the order of 30 GB of RAM for a 3 Gb mammalian genome index, so run that stage on a large-memory machine or skip it. ### Can I run QC without aligning the reads? Yes. FastQC, fastp, and MultiQC catch most instrument-level and library-level problems without a reference. Alignment adds mapping rate, strandedness, rRNA fraction, and coverage-uniformity checks, which is why the template makes it optional rather than absent. ### What is the difference between FastQC and MultiQC? FastQC produces one report per FASTQ file, and MultiQC aggregates reports from FastQC, fastp, STAR, and many other tools into a single HTML page with one table per metric across all samples. You need both: FastQC generates the data, MultiQC makes it comparable. ### Why did MultiQC find no samples? MultiQC scanned a directory containing only files it does not recognize. Check the MULTIQC task directory under `work/` to confirm the FastQC zip archives and fastp JSON files were staged, then run `multiqc` manually on that folder to see its parse warnings. ### How do I adapt the template for single-end reads? Leave `fastq_2` empty in the samplesheet, set `single_end: true` in the meta map, and make the second input conditional in FASTP and STAR, for example building the second `-in2` argument only when `meta.single_end` is false. The FastQC and MultiQC stages need no changes. ### Which trimming settings should I use? Start with adapter auto-detection, a Q20 qualified-base threshold, and a 36 nucleotide minimum length, then compare raw versus trimmed mapping rates on a few samples. If trimming barely changes mapping, the settings are safe. If it lowers mapping, loosen them. ### How does this template relate to nf-core/rnaseq? nf-core/rnaseq is a production-grade community pipeline with far more stages, tests, and documentation than any template. Use this template to learn the architecture and prototype quickly, then move real studies to nf-core/rnaseq or a similarly maintained pipeline. ### How do I run the pipeline on an HPC cluster? Add a profile that sets `process.executor = 'slurm'` (or your scheduler's name), keep the conda or container profile, and launch with `nextflow run main.nf -profile slurm`. Nextflow submits each task as a scheduler job using the CPU, memory, and time values from the process labels. ## Sources 1. [nf-rnaSeqCount: A Nextflow pipeline for obtaining raw read counts from RNA-seq data.](https://pubmed.ncbi.nlm.nih.gov/35574063/) 2. [PeakPrime: a peak-guided primer design pipeline for target enrichment in 3'-end RNA-seq.](https://pubmed.ncbi.nlm.nih.gov/41919010/)