nf-core/rnaseq

When writing tests, follow these snapshotting conventions:

• Non-deterministic outputs (e.g., STAR, Salmon, Kallisto, RSEM, HISAT2 indices; qualimap reports): Snapshot by file-name-only using getSnapshot() filtering.
• Deterministic text outputs: Snapshot using md5.
• Verbose JSON test output: Snapshot the .md5() of the result rather than inlining the JSON.
• Avoid: Do not snapshot timestamps or paths containing hash directories.

If a pipeline-level snapshot is unstable (e.g., due to timestamps, varying log order, or multiqc_data paths), add the file glob patterns to tests/.nftignore (or tests/.nftignore_rustqc for the RustQC variant) to exclude them from nf-test snapshots.

Generating indices (STAR, Salmon, etc.) is computationally expensive. To avoid redundant computation across runs:

1. Use --save_reference to retain generated indices.
2. Store these indices in a central location.
3. Reuse them in subsequent runs by specifying the path directly (e.g., --star_index '/path/to/STAR/index/') or by including them in a custom --genome catalogue.

dupRadar relates duplication rates to gene expression levels to identify technical duplication.

• Good Sample: High duplication should only occur for highly expressed genes.
• Bad Sample: High duplication occurs across all genes, regardless of transcription level.

Once differentially expressed genes are identified, use Over Representation Analysis (ORA) to assign biological meaning to the results.

ORA identifies biological pathways, molecular functions, or cellular processes that are enriched within your DE gene list, typically using the Gene Ontology (GO) database. The output is often visualized as an enrichment plot showing terms like 'clathrin-coated pit' or 'endoplasmic reticulum lumen', which helps link gene-level changes to specific cellular structures or biological processes.

StringTie assembles RNA-Seq alignments into transcripts. The pipeline supports two modes:

1. Reference-guided (Default): Performs a single pass of quantification against the provided GTF.
2. De novo workflow: Triggered by setting --stringtie_ignore_gtf. This follows a three-stage process:
   • Per-sample de novo assembly (<SAMPLE>.denovo.transcripts.gtf).
   • Merging per-sample assemblies into a unified annotation (stringtie_merge.gtf).
   • Re-quantifying all samples against the merged annotation. This mode is recommended to ensure expression values are directly comparable across samples.

The pipeline quantifies RNA-sequenced reads relative to genes/transcripts in the genome and normalizes the data.

Important Limitations:

• The pipeline does not perform statistical comparisons between samples (it does not calculate FDR or P-values).

Downstream Usage:

• The output files (gene expression matrices, etc.) are designed to be used in statistical environments like R or Julia.
• For differential abundance testing, it is recommended to use the nf-core/differentialabundance pipeline.

The choice of statistical test depends on your experimental design:

• Wald Test: Best for simple, two-condition comparisons (e.g., treated vs. untreated). It calculates a z-statistic by dividing the log2 fold change by its standard error and evaluates it against a standard normal distribution.
• Likelihood Ratio Test (LRT): Best for complex, multi-factor designs. It compares a 'full model' (including the factor of interest) against a 'reduced model' (excluding it). This is used to identify genes that change in any direction across multiple sample classes.

The pipeline supports several quantification strategies depending on your choice of aligner:

• STAR + Salmon/RSEM: Use --aligner star_salmon or --aligner star_rsem for standard genomic alignment and quantification.
• Pseudo-alignment: For faster quantification without full alignment, use --pseudo_aligner salmon or --pseudo_aligner kallisto.

Note: Quantification is not performed when using --aligner hisat2 because it lacks an appropriate option for accurate expression estimates in this pipeline. HISAT2 is available for users who prefer its specific alignment and QC outputs.

The pipeline can scan reads for contamination from other species using either Sylph or Kraken2 (with optional Bracken refinement).

• Default behavior: Screening runs on aligner-unmapped reads (reads rejected by STAR or HISAT2).
• Custom input: To inspect reads entering the alignment stage (e.g., to evaluate BBSplit filtering), set --contaminant_screening_input trimmed.
• Sylph: Faster and more memory-efficient with lower false-positive rates, but only provides overall abundance estimates and cannot assign a specific percentage of reads as unclassified.
• Kraken2: Highly dependent on the database used. It is crucial to ensure the host genome/transcriptome is included in the database to avoid false positives. It is recommended to always use Bracken alongside Kraken2 to refine abundance estimates.

The read_distribution.py script calculates how mapped reads are distributed across genomic features (e.g., CDS exon, 5’UTR exon, 3’ UTR exon, Intron, Intergenic).

For standard RNA-seq, you typically expect a high proportion of exonic reads (CDS_Exons). A high proportion of intronic reads may indicate DNA contamination, though it can be expected in total RNA preparations.