You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
# RNA-seq Data Processing, Quantification & Annotation Pipeline
10
-
A [Snakemake 8](https://snakemake.readthedocs.io/en/stable/) workflow for end-to-end processing, quantification, and annotation of gene expression for RNA-seq experiments, starting from single- or paired-end reads within raw/unaligned/unmapped [uBAM](https://gatk.broadinstitute.org/hc/en-us/articles/360035532132-uBAM-Unmapped-BAM-Format) files, including a comprehensive MultiQC report.
10
+
A [Snakemake](https://snakemake.readthedocs.io/en/stable/) workflow for end-to-end processing, quantification, and annotation of gene expression for RNA-seq experiments, starting from single- or paired-end reads within raw/unaligned/unmapped [uBAM](https://gatk.broadinstitute.org/hc/en-us/articles/360035532132-uBAM-Unmapped-BAM-Format) files, including a comprehensive MultiQC report.
11
11
12
-
> [!NOTE]
12
+
> [!NOTE]
13
13
> This workflow adheres to the module specifications of [MrBiomics](https://github.com/epigen/MrBiomics), an effort to augment research by modularizing (biomedical) data science. For more details, instructions, and modules check out the project's repository.
14
14
>
15
15
> ⭐️ **Star and share modules you find valuable** 📤 - help others discover them, and guide our future work!
16
16
17
17
> [!IMPORTANT]
18
-
> **If you use this workflow in a publication, please don't forget to give credit to the authors by citing it using this DOI [10.5281/zenodo.15119355](https://doi.org/10.5281/zenodo.15119355) and acknowledging the [rna-seq-star-deseq2 workflow](https://github.com/snakemake-workflows/rna-seq-star-deseq2) DOI [10.5281/zenodo.4737358](https://doi.org/10.5281/zenodo.4737358) from which some structure and code were adapted.**
18
+
> **If you use this workflow in a publication, please don't forget to give credit to the authors by citing it using this DOI [10.5281/zenodo.15119355](https://doi.org/10.5281/zenodo.15119355) and acknowledge the [rna-seq-star-deseq2 workflow](https://github.com/snakemake-workflows/rna-seq-star-deseq2) DOI [10.5281/zenodo.4737358](https://doi.org/10.5281/zenodo.4737358) from which some structure and code were adapted.**
@@ -47,15 +52,15 @@ This is a template for the Methods section of a scientific publication and is in
47
52
48
53
**Quantification.** Gene expression quantification was performed on the filtered and trimmed reads. The STAR aligner (ver) [ref] was utilized in `--quantMode GeneCounts` mode to count reads overlapping annotated genes based on the Ensembl [ver `config: ref: release`] gene annotation for the [`config: ref: species`] genome (build [`config: ref: build`]) and using parameters [`config["star_args"]`]. Library strandedness (`[none/yes/reverse]`, specified per sample in `config/annotation.csv`) was accounted for during counting. Counts for individual samples were aggregated into a single gene-by-sample count matrix using a custom Python script utilizing the pandas library (ver) [ref]. Quality control metrics from various tools, including fastp (ver) [ref], RSeQC (ver) [ref] and STAR (ver) [ref], were aggregated using MultiQC (ver) [ref].
49
54
50
-
**Annotation.** Gene annotations from Ensembl were retrieved using the R package biomaRt (ver) [ref]. Annotations included Ensembl gene ID, gene symbol (`external_gene_name`), gene biotype, and description. Additionally, exon-based GC content (`exon_gc`) and cumulative exon length (`exon_length`) were calculated for each gene using a custom R function adapted from EDASeq (ver) [ref], leveraging biomaRt (ver) [ref] to fetch exon coordinates and sequences. This exon-based approach was chosen as sequencing reads in poly(A)-selected libraries primarily derive from exonic regions, making these metrics more appropriate for downstream bias correction (e.g., Conditional Quantile Normalization - CQN) than whole-gene metrics. A sample annotation file was generated, integrating input annotations with QC metrics.
55
+
**Annotation.** Gene annotations from Ensembl were retrieved using the R package biomaRt (ver) [ref]. Annotations included Ensembl gene ID, gene symbol (`external_gene_name`), gene biotype, and description. Additionally, exon-based GC content (`exon_gc`) and cumulative exon length (`exon_length`) were calculated for each gene based on the genome GFT and FASTA files using a custom R function, leveraging GenomicRanges (ver) [ref] and rtracklayer (ver) [ref]. This exon-based approach was chosen because only exonic reads are quantified during alignment and contribute to the count matrix (also sequencing reads in poly(A)-selected libraries primarily derive from exonic regions), making these metrics more appropriate for downstream bias correction (e.g., Conditional Quantile Normalization - CQN) than whole-gene metrics (i.e., including introns). A sample annotation file was generated, integrating input annotations with QC metrics.
51
56
52
-
The processing and quantification described here was performed using a publicly available Snakemake [ver] (ref) workflow [[10.5281/zenodo.15119355](https://doi.org/10.5281/zenodo.15119355)], which adopted code from anotehr workflow (ref) [10.5281/zenodo.4737358](https://doi.org/10.5281/zenodo.4737358).
57
+
The processing and quantification described here was performed using a publicly available Snakemake [ver] (ref) workflow [[10.5281/zenodo.15119355](https://doi.org/10.5281/zenodo.15119355)], which adapted code from another workflow (ref) [10.5281/zenodo.4737358](https://doi.org/10.5281/zenodo.4737358).
53
58
54
59
# 🚀 Features
55
60
This workflow offers several key advantages for RNA-seq analysis over existing pipelines:
56
61
57
62
***MrBiomics Module:** Designed for modularity and seamless integration with other [MrBiomics](https://github.com/epigen/MrBiomics/) analysis workflows (e.g., filtering/normalization, differential expression, unsupervised analysis) and includes example analysis recipes.
58
-
***Exon-Centric Annotation:** Calculates *exon-based* GC content and length, providing more accurate metrics for downstream bias correction (like CQN) for typical poly(A)-selected RNA-seq libraries.
63
+
***Exon-Centric Annotation:** Calculates *exon-based* GC content and length, providing more accurate metrics for downstream bias correction (like CQN).
59
64
***Robust Input Handling:** Starts directly from raw uBAM files and includes automated checks to ensure consistency between annotated read types (single/paired) and BAM file content.
60
65
***Efficiency Focused:** Optimized for performance and disk space using data streaming between processing steps and temporary intermediate files.
61
66
***User-Friendly:** Offers comprehensive documentation, clear configuration, standard Snakemake usage, and practical usage tips including detailed QC guidelines.
@@ -73,22 +78,23 @@ The workflow performs the following steps that produce the outlined results:
73
78
- De-interleaves the filtered FASTQ stream into separate compressed R1 and R2 files for paired-end data, or compresses directly for single-end data using shell commands and `pigz`.
74
79
> [!NOTE]
75
80
> `fastp` adapter auto-detection is disabled because we use STDIN mode (i.e., stream the data through pipes) to be disk space efficient.
81
+
> We do not deduplicate aligned reads. ["...the computational removal of duplicates does improve neither accuracy nor precision and can actually worsen the power and the False Discovery Rate (FDR) for differential gene expression."](https://www.nature.com/articles/srep25533)
76
82
-**Quantification:**
77
83
- Uses STAR `GeneCounts` to quantify reads per gene based on the specified Ensembl reference genome and annotation (`star/{sample}/`).
78
84
- Handles unstranded, forward-stranded, and reverse-stranded library protocols based on the `strandedness` column.
79
85
- Aggregates counts into a single matrix (`counts/counts.csv`).
- Retrieves gene annotations (Ensembl ID, gene symbol, biotype, description) from Ensembl using `biomaRt`.
83
-
- Calculates **exon-based** GC content and cumulative exon length for each gene, suitable for poly(A) selected libraries.
84
-
- Outputs a sample annotation table containing sample-wise general MultiQC statistics (`counts/sample_annotation.csv`).
89
+
- Calculates **exon-based** GC content and cumulative exon length for each gene using the genome's GFT and FASTA files.
90
+
- Outputs a sample annotation table containing sample-wise general MultiQC statistics (`counts/sample_annotation.csv`).
85
91
> [!NOTE]
86
92
> Gene annotation can take a while since it depends on the availability of external data sources accessed via `biomaRt`.
87
93
>
88
-
> GC-content and length are **exon-based**: In poly(A)‑selected libraries (such as Illumina TruSeq, Smart-seq or QuantSeq), the sequencing reads mainly come from exonic regions. Therefore, potential correction for GC bias and gene length should ideally use exon‑level GC content and effective exon length rather than whole‑gene metrics that include introns.
94
+
> GC-content and length are **exon-based**, because we only use exonic reads during the count matrix generation. Furthermore, in poly(A)‑selected libraries (such as Illumina TruSeq, Smart-seq or QuantSeq), the sequencing reads mainly come from exonic regions. Therefore, potential correction for GC bias and gene length should use exon‑level GC content and effective exon length rather than whole‑gene metrics that include introns.
89
95
-**QC & Reporting:**
90
-
- Employs RSeQC tools to generate key quality metrics like strand specificity and read distribution across genomic features (`rseqc/`).
91
-
- Aggregates QC metrics from fastp, STAR and RSeQC into a single report using MultiQC (`report/multiqc_report.html`) with [AI summaries](https://seqera.io/blog/ai-summaries-multiqc/).
96
+
- Employs `RSeQC` tools to generate key quality metrics like strand specificity and read distribution across genomic features (`rseqc/`).
97
+
- Aggregates QC metrics from `fastp`, `STAR` and `RSeQC` into a single report using `MultiQC` (`report/multiqc_report.html`) with [AI summaries](https://seqera.io/blog/ai-summaries-multiqc/).
92
98
93
99
---
94
100
@@ -111,15 +117,14 @@ The workflow produces the following directory structure:
111
117
> [!IMPORTANT]
112
118
> Resources are downloaded automatically to `resources/{config::project_name}/rnaseq_pipeline/)`, are large (>`3GB`) and have to be manually removed if not needed anymore.
113
119
114
-
115
120
# 🛠️ Usage
116
121
117
122
Here are some tips for the usage of this workflow:
118
123
119
124
- Configure and run the workflow first for a few samples.
120
125
- Once everything works, run it for all samples.
121
-
- If you are unsure about the memory requirements for alignment (the most memory intensive step) provide the Snakemake parameter `--retries X`, where `X` denotes the number of retries, and with every retry the memory is increased to `attempts * config:mem`.
122
-
- To save disk space the intermediate gzipped FASTQ files are marked as temporary using Snakemake's `temp()` directive. To remove them upon successful completion you have to include the `--delete-temp-output` flag in your Snakemake command.
126
+
- If you are unsure about the memory requirements for alignment (the most memory intensive step) provide the Snakemake parameter `--retries X`, where `X` denotes the number of retries, and with every retry the memory is increased to `attempts * config::mem`.
127
+
- To save disk space the intermediate gzipped `FASTQ` files are marked as temporary using Snakemake's `temp()` directive. To remove them upon successful completion you have to run Snakemake with the `--delete-temp-output` flag.
123
128
124
129
This workflow is written with Snakemake and its usage is described in the [Snakemake Workflow Catalog](https://snakemake.github.io/snakemake-workflow-catalog?usage=epigen/rnaseq_pipeline).
125
130
@@ -133,9 +138,10 @@ Explore detailed examples showcasing module usage in comprehensive end-to-end an
133
138
# 🔍 Quality Control
134
139
Below are some guidelines for the manual quality control of each sample using the generated `MultiQC` report, but keep in mind that every experiment/dataset is different. Thresholds are general suggestions and may vary based on experiment type, organism, and library prep.
-**Alignment Scores & Gene Counts (STAR):** High proportion of uniquely mapped reads assigned to exonic regions (e.g., >60-70% for poly(A) mRNA-seq). Low rates could suggest gDNA contamination or high intronic reads. In case of many intronic reads or non-poly(A) mRNA-seq protocols do not use exon-based gene annotations (gc-content and length).
138
-
-**Read Quality (fastp):** High average quality scores across reads after trimming. Ensure effective adapter removal.
141
+
-**Read Depth (STAR)**: Count of `(Uniquely) Mapped Reads` >10M is minimum, >20M reads is optimal for differential expression analysis.
-**Alignment Scores & Gene Counts (STAR):** High proportion of uniquely mapped reads assigned to exonic regions (e.g., >60-70% for poly(A) mRNA-seq). Low rates could suggest gDNA contamination or high intronic reads. In case of many intronic reads or non-poly(A) mRNA-seq protocols do not use exon-based gene annotations (gc-content and length) downstream.
144
+
-**Read Quality (fastp):**`% > Q30` (=Percentage of bases with Phred score > 30, after filtering/trimming) > 90%. High average quality scores across reads after trimming. Ensure effective adapter removal.
139
145
-**Library Complexity (fastp/RSeQC):** % Duplication Rate should not be excessively high (highly variable, interpret in context of expression). Very high rates might indicate low input or PCR issues.
140
146
-**Strand Specificity (RSeQC):** For stranded protocols, >90-95% reads should match the expected strand.
141
147
- Inspect [**Genome Browser Tracks**](https://github.com/epigen/genome_tracks/) using UCSC Genome Browser (online) or IGV (local)
0 commit comments