RNA-seq Snakemake Pipeline

A comprehensive RNA-seq analysis pipeline built with Snakemake that supports multiple data types including paired-end, single-end, and nanopore sequencing data.

This workflow powered the Meta-Analysis of VSG Expression reanalyzing 454 RNA-seq runs (168 single-end, 280 paired-end, and 6 nanopore datasets) from 31 publications, spanning 35 experimental factors across 78 experiments.

Features

• Multi-format support: Paired-end, single-end and nanopore data
• Flexible input sources: Local files and FTP downloads
• Quality control: Comprehensive QC with fastp, samtools, and Qualimap
• Run tag grouping: Combine multiple samples for joint analysis
• Checksum verification: Ensure data integrity
• Containerized: Singularity support
• Benchmarking: Detailed performance metrics
• Modular design: Rules are automatically included based on your data
• On-demand directories: Sample folders are created only when processing starts and cleaned up when finished

Pipeline Architecture

Input Data
    ↓
Acquisition & Checksum Verification
(local / ftp)
    ↓
Quality Filtering* (fastp)
    ↓
Alignment (bowtie2/minimap2)
    ↓
BAM Merging for technical replicates
    ↓
Mark Duplicates* (Picard)
    ↓
Quality Control (samtools + Qualimap*)
    ↓
Coverage Tracks (bamCoverage)
    ↓
Feature Counting (featureCounts)
    ↓
RPKM Calculation
    ↓
Results & Cleanup

* Step skipped for nanopore data

Data Type Processing:

• Paired/Single/Nanopore: Full pipeline with bowtie2
• Nanopore: Direct minimap2 alignment, skips QC/duplicate steps

Supported Data Types

1. Paired-end Illumina Data

• read_type: paired
• Description: Standard paired-end sequencing with separate R1 and R2 files
• Required fields: file_path_1 (R1), file_path_2 (R2)
• Processing: Quality filtering with fastp → Alignment with bowtie2 → Standard QC

2. Single-end Illumina Data

• read_type: single
• Description: Single-end sequencing with one file per sample
• Required fields: file_path_1 only
• Processing: Quality filtering with fastp → Alignment with bowtie2 → Standard QC

3. Nanopore Long-read Data

• read_type: nanopore
• Description: Oxford Nanopore long-read sequencing data
• Required fields: file_path_1 only
• Processing: Direct alignment with minimap2 → Coverage analysis (skips quality filtering and duplicate marking)
• Special considerations:
  • Uses minimap2 with splice-aware alignment (-ax splice)
  • Skips fastp quality filtering (nanopore reads handled differently)
  • Skips Picard MarkDuplicates (not applicable to long reads)
  • Skips Qualimap RNA-seq QC (optimized for short reads)

Quick Start

1. Installation

# Clone the repository
git clone <repository-url>
cd rna-seq-pipeline

# Install Snakemake (if not already installed)
conda install -c conda-forge -c bioconda snakemake==9.1.9

2. Prepare your sample sheet

Create a samples.csv file with your sample information:

sample_name,read_type,source_type,file_path_1,file_path_2,checksum_1,checksum_2,run_tag
sample1,paired,local,/data/sample1_R1.fastq.gz,/data/sample1_R2.fastq.gz,abc123,def456,experiment_A
sample2,single,ftp,ftp://server.com/sample2.fastq.gz,,xyz789,,experiment_B
sample4,nanopore,local,/data/sample4_nanopore.fastq.gz,,pqr678,,long_read_set

3. Configure the pipeline

Edit config.yaml to specify your reference genome and processing parameters:

# Core directories
processing_dir: "processing"
results_dir: "results"

# Sample information
samples_csv: "samples.csv"

# Reference genome (FASTA must end with .fa)
reference_fasta: "reference/genome/genome.fa"
gtf_file: "reference/genome/annotation.gtf"

# Processing parameters
cores_align: 8
feature_type: "CDS"

The workflow stages the FASTA inside the processing directory and builds the Bowtie2 index automatically, so you only need to provide the reference FASTA (with a .fa extension) and annotation file. The pipeline will exit with an error if the FASTA does not use the required .fa suffix.

Annotation Format Support: The pipeline accepts both GTF and GFF annotation files. When a GFF file is provided, it is implicitly converted to GTF format during processing.

gffread in.gff3 -T -o out.gtf

Users should ensure that the attributes and feature types specified in the configuration (e.g., feature_type, attribute_type) are present in the original file and will be preserved in the GFF-to-GTF conversion. This is particularly important for custom annotations or non-standard attribute names.

4. Run the pipeline

# Dry run to check the workflow
snakemake -n

# Run the pipeline
snakemake --cores 8

Running with Singularity

When running with Singularity on shared HPC filesystems, the container needs access to the full host paths used by the workflow. Use absolute paths in your configuration and bind the project directory when launching Snakemake.

snakemake --cores 40 \
  --configfile config_rit.yaml \
  --use-singularity \
  --singularity-args "-B /gpfs/uod-scale-01/cluster/majf_lab/mtinti --pwd /gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq"

In the config file, provide full paths for the core directories. For example:

# Core directories
processing_dir: "/tmp/3549775.1.rhel9.q/processing"
results_dir: "/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/results_singularity"
benchmark_dir: "/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/benchmarks"

Also use full absolute paths in the samples sheet (file_path_1 and file_path_2) so files are visible inside the container runtime. Example:

sample_name,read_type,source_type,file_path_1,file_path_2,checksum_1,checksum_2,run_tag
SAMPLE_rit1,paired,local,/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/indata/TbRIT-8166/V350168884_L02_B5GTBRjalrRAACA-4_1.fq.gz,/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/indata/TbRIT-8166/V350168884_L02_B5GTBRjalrRAACA-4_2.fq.gz,db36979b3164a2a29c9af85e6b3072fd,01aad38e21bc8ea452612b10037a58a7,TbRIT-8166
SAMPLE_rit2,paired,local,/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/indata/TbRIT-8166/V350194510_L02_B5GTBRjalrRAACA-4_1.fq.gz,/gpfs/uod-scale-01/cluster/majf_lab/mtinti/myRna-seq/indata/TbRIT-8166/V350194510_L02_B5GTBRjalrRAACA-4_2.fq.gz,876109f629db0fef6c93918137d467a0,6871417faaf83beae9c72051a12519ce,TbRIT-8166

You can pull the Singularity image from the Docker release by following the instructions in: mtinti/myRna-seq-docker.

# briefly
singularity pull rna_seq.sif docker://mtinti/rna_seq:latest

Testing the Pipeline

You can test the pipeline with the included test dataset:

snakemake --cores 4 --use-conda --config \
            processing_dir="tests/test_counts/test_out/processing" \
            results_dir="tests/test_counts/test_out/results" \
            benchmark_dir="tests/test_counts/test_out/benchmarks" \
            reference_fasta="tests/test_counts/genome/random_genome.fa" \
            gtf_file="tests/test_counts/genome/annotation.gtf" \
            samples_csv="test_samples_local.csv" \
            cleanup_processing='True' \
            cores_align=2 cores_coverage=2 cores_fastp=2 cores_featurecounts=2 \
            cores_flagstat=2 cores_nanopore_align=2 max_cores=2

Make sure the FASTA you reference ends with .fa; other extensions will cause the workflow to abort during validation.

Expected output from the test dataset:

Gene	Counts
gene1	10
gene2	8
gene3	2

Sample CSV Format

Your samples.csv file should include these columns:

Column	Required	Description
sample_name	Yes	Unique identifier for each sample
read_type	Yes	One of: paired, single, nanopore
source_type	Yes	local, ftp, or sra_paired
file_path_1	Yes	Path to first/only FASTQ file
file_path_2	Conditional	Path to second FASTQ file (required for paired only)
checksum_1	Optional	MD5 checksum for file_path_1
checksum_2	Optional	MD5 checksum for file_path_2
run_tag	Optional	Group samples for combined analysis

Example CSV entries:

sample_name,read_type,source_type,file_path_1,file_path_2,checksum_1,checksum_2,run_tag
illumina_pe_01,paired,local,/data/sample1_R1.fastq.gz,/data/sample1_R2.fastq.gz,abc123,def456,experiment_A
illumina_se_01,single,ftp,ftp://server.com/sample2.fastq.gz,,xyz789,,experiment_B
nanopore_01,nanopore,local,/data/sample4_nanopore.fastq.gz,,pqr678,,long_read_set

Using SRA Data

The pipeline can download reads directly from NCBI SRA. Set source_type to sra_paired and specify the accession in file_path_1:

sample_name,read_type,source_type,file_path_1,file_path_2,checksum_1,checksum_2
MySRA,paired,sra_paired,SRR12345678,,,

SRA reads are downloaded with fasterq-dump, reformatted, and gzipped before entering the standard workflow. Each download touches a .download_complete file to signal that the FASTQ files are fully written before validation and further processing.

Processing Differences by Data Type

Quality Control Steps

Step	Paired	Single	Nanopore
fastp filtering	✅ Yes	✅ Yes	❌ No
Alignment tool	bowtie2	bowtie2	minimap2
Mark duplicates	✅ Yes	✅ Yes	❌ No
Qualimap BAM QC	✅ Yes	✅ Yes	❌ No
Qualimap RNA-seq QC	✅ Yes	✅ Yes	❌ No
FeatureCounts QC	✅ Yes	✅ Yes	✅ Yes

Alignment Parameters

Illumina data (paired/single):

• Tool: bowtie2
• Options: --very-sensitive-local
• Read groups: Automatically added for Picard compatibility

Nanopore data:

• Tool: minimap2
• Preset: map-ont (configurable via nanopore_preset)
• Options: -ax splice -uf -k14 -G 10000 (splice-aware alignment)
• Read groups: Platform set to "ONT"

Configuration

Core Configuration (config.yaml)

# Core directories
processing_dir: "processing"
results_dir: "results"

# Sample information
samples_csv: "samples.csv"
selected_samples: []  # Empty list means all samples are processed

# Reference genome (FASTA must end with .fa)
reference_fasta: "reference/genome/genome.fa"
gtf_file: "reference/genome/annotation.gtf"

# Processing parameters
cores_align: 8
cores_coverage: 8
cores_featurecounts: 8
cores_fastp: 8
cores_default: 1
feature_type: "CDS"

# Nanopore-specific settings
nanopore_preset: "map-ont"        # Minimap2 preset
cores_nanopore_align: 8           # CPU cores for nanopore alignment
mem_nanopore_align: 16000         # Memory for nanopore alignment (MB)

# Mark Duplicates parameters
remove_duplicates: false  # Set to true to remove duplicates, false to only mark

# QC parameters
qualimap_memory: "10G"

# Coverage track parameters
coverage_bin_size: 10
coverage_normalize: "RPKM"
min_mapping_quality: 2

# File handling preferences
copy_fastq: false        # Whether to copy FASTQ files to results
copy_bam: false          # Whether to copy BAM files to results
copy_benchmarks: true    # Whether to copy benchmark files
copy_logs: true          # Whether to copy log files
cleanup_processing: false  # Remove processing directory and reference after completion

# Container support
singularity_image: "/path/to/container.sif"

Resource Configuration

# Resource limits
max_resources:
  network: 3    # Maximum concurrent FTP downloads
  io: 1000      # Maximum concurrent IO operations

Output Structure

results/
├── {sample_or_run_tag}/
│   ├── {sample}_all.bw                    # Coverage track (all reads)
│   ├── {sample}_unique.bw                 # Coverage track (unique reads)
│   ├── {sample}_counts_*_all.txt          # Feature counts (all reads)
│   ├── {sample}_counts_*_unique.txt       # Feature counts (unique reads)
│   ├── qc/                                # Quality control reports
│   │   ├── fastp/                         # fastp reports (if applicable)
│   │   ├── flagstat/                      # samtools flagstat
│   │   ├── stats/                         # samtools stats
│   │   ├── {sample}_qualimap_bam/         # Qualimap BAM QC
│   │   ├── {sample}_qualimap_rnaseq/      # Qualimap RNA-seq QC (if applicable)
│   │   ├── markduplicates/                # Picard metrics (if applicable)
│   │   └── feature_counts/                # featureCounts summaries
│   ├── benchmarks/                        # Performance benchmarks (optional)
│   └── logs/                              # Log files (optional)
└── copy_complete_all.txt                  # Pipeline completion marker

If cleanup_processing is enabled, the pipeline removes each sample's processing directory after the results are copied. The processing/reference directory is also removed at the end of the run.

Run Tags and Sample Grouping

The run_tag column allows you to group samples for combined analysis:

sample_name,read_type,source_type,file_path_1,run_tag
sample_A1,paired,local,/data/A1_R1.fastq.gz,experiment_1
sample_A2,paired,local,/data/A2_R1.fastq.gz,experiment_1
sample_B1,nanopore,local,/data/B1.fastq.gz,long_reads

Important notes:

• All samples within a run tag group must have the same read_type
• BAM files are merged by run tag before downstream analysis
• Final outputs use the run tag name, not individual sample names
• If no run_tag is specified, each sample is processed independently

Advanced Usage

Running Specific Samples

# Process only specific samples
snakemake --cores 8 --config selected_samples="['sample1','sample2']"

Container Usage

# With Singularity
snakemake --cores 8 --use-singularity

Resource Management

# Limit concurrent FTP downloads
snakemake --cores 8 --resources network=2

Resource Management

# Limit concurrent FTP downloads
snakemake --cores 8 --resources network=2

Standalone MultiQC

A dedicated MultiqcOnly.smk workflow is provided to generate a consolidated MultiQC report from existing pipeline outputs. This standalone workflow expects results organized in the same format as the main RNA-seq pipeline and can be executed with:

snakemake --snakefile MultiqcOnly.smk \
  --config results_dir=/path/to/pipeline/results/ \
  samples_csv=/path/to/samples.csv \
  --cores 1 --use-conda

This allows MultiQC to be run independently on previously generated results.