RNA-Seq Quantification Pipeline

Version 1.2.0

A reproducible Snakemake pipeline for RNA-Seq read processing and gene-level quantification. Handles both paired-end and single-end libraries with automatic or user-specified strandness detection.

Repository Structure

.
├── config/
│   ├── config.yaml          # Pipeline parameters & reference paths
│   └── samplesheet.csv      # Sample manifest (edit this)
├── workflow/
│   ├── Snakefile             # Main entry point (common code + includes)
│   ├── rules/
│   │   ├── qc.smk           # FastQC (raw & trimmed)
│   │   ├── trim.smk         # fastp (adapter/quality trimming)
│   │   ├── filter.smk       # BBDuk (rRNA filtering)
│   │   ├── align.smk        # HISAT2 alignment
│   │   └── count.smk        # Strandness resolution + featureCounts
│   └── scripts/
│       └── featurecounts.R   # R script for Rsubread::featureCounts
├── environment.yaml          # Conda environment specification
├── LICENSE
└── README.md

Pipeline Overview

Raw FASTQ ─→ FastQC ─→ fastp (trim) ─→ FastQC ─→ BBDuk (rRNA filter) ─→ HISAT2 (align) ─→ featureCounts
                                                                                              ↑
                                                                          RSeQC infer_experiment.py
                                                                          (auto-detect strandness)

Step	Tool	Output directory
1. Raw QC	FastQC	results/01_fastqc_raw_reports/
2. Adapter/quality trimming	fastp	results/02_fastp_output/
3. Post-trim QC	FastQC	results/03_fastqc_trimmed_reports/
4. rRNA removal	BBDuk (BBTools)	results/04_rrna_filter/
5. Alignment	HISAT2 + samtools	results/05_sorted_bams/
6. Strandness resolution	RSeQC / samplesheet	results/06_strandness/
7. Gene counting	Rsubread featureCounts	results/07_counts/

Quick Start

1. Clone the repository

git clone https://github.com/<your-username>/snakemake-rnaseq-pipeline.git
cd snakemake-rnaseq-pipeline

2. Create the conda environment

conda env create -f environment.yaml
conda activate rnaseq_pipeline

3. Prepare reference files

The pipeline requires:

• HISAT2 genome index (e.g., GRCh38)
• GTF annotation file
• BED gene model file (for strandness inference via RSeQC)
• rRNA reference FASTA (e.g., humanRrna_mod.fasta or BBMap's ribokmers.fa.gz)

Download or generate these and note their paths.

4. Configure the pipeline

Edit config/config.yaml — replace all <REF_DIR> placeholders with your actual paths:

hisat2_index:    "/path/to/GRCh38_hisat2_index"
gene_bed_file:   "/path/to/gene_model.bed"
gtf_annotation:  "/path/to/annotation.gtf"
bbduk_rrna_ref:  "/path/to/humanRrna_mod.fasta"

5. Create your samplesheet

Edit config/samplesheet.csv following the template format:

sample,fastq1,fastq2,strandness
sample_PE_rep1,/path/to/R1.fastq.gz,/path/to/R2.fastq.gz,2
sample_SE_rep1,/path/to/reads.fastq.gz,,

Column	Description
sample	Unique sample identifier
fastq1	Absolute path to R1 FASTQ (required)
fastq2	Absolute path to R2 FASTQ (empty for single-end)
strandness	0 = unstranded, 1 = sense, 2 = antisense, empty = auto-infer

6. Run the pipeline

Run from the project root directory:

snakemake --cores 12

Dry run (check what will execute without running):

snakemake -n

Note: Snakemake automatically discovers workflow/Snakefile when run from the project root. All output paths are relative to this working directory.

Strandness Handling

The pipeline supports three modes per sample:

strandness value	HISAT2 flag (PE/SE)	featureCounts strandSpecific
0 (unstranded)	(none)	0
1 (sense)	--rna-strandness FR / F	1
2 (antisense)	--rna-strandness RF / R	2
(empty)	(none) during alignment; inferred post-alignment	Auto-detected via RSeQC

When the strandness column is left empty, the pipeline:

1. Aligns reads without a strandness flag
2. Runs infer_experiment.py (RSeQC) on the resulting BAM
3. Classifies as sense/antisense/unstranded using a 0.75 threshold
4. Uses the resolved code for featureCounts quantification

Output

The primary output is results/07_counts/counts.txt — a tab-separated gene×sample count matrix ready for differential expression analysis (DESeq2, edgeR, limma, etc.).

Additional outputs:

• results/07_counts/counts.txt.summary — featureCounts assignment summary
• results/07_counts/counts.rds — R serialized featureCounts object
• results/logs/ — per-rule log files for debugging

Configuration

All parameters are in config/config.yaml:

Parameter	Default	Description
threads	6	Number of threads per rule
min_read_length	36	Minimum read length after trimming (fastp)
bbduk_k	31	K-mer size for rRNA filtering
bbduk_hdist	1	Hamming distance for BBDuk

Requirements

All tools are specified with exact versions in environment.yaml:

• Snakemake ≥ 9.11.2
• fastp 1.0.1
• FastQC 0.12.1
• BBTools (BBDuk) 39.06
• HISAT2 2.2.1
• samtools 1.22.1
• RSeQC 5.0.4
• pandas 2.3.2
• R 4.4.1 + Rsubread 2.20.0

Troubleshooting

Problem	Solution
ERROR: config['...'] still contains <REF_DIR>	Edit config/config.yaml — replace all <REF_DIR> placeholders
infer_experiment.py fails	Ensure RSeQC is installed and gene_bed_file path is correct
Out of memory on BBDuk	Increase -Xmx in the bbduk_* rules or reduce bbduk_k
HISAT2 index not found	Run hisat2-build genome.fa genome_index to create the index

License

MIT

Citation

If you use this pipeline, please cite:

( manuscript in preparation)