22G RNA Counting Pipeline

A Nextflow DSL2 pipeline for processing UMI-extracted small RNA-seq data to identify and quantify 22G RNAs in C. elegans.

Pipeline Overview

This pipeline processes UMI-extracted FASTQ files through the following steps:

1. Cutadapt Trimming - Trim adapters and filter by length (default: 22nt exactly)
2. 22G Filtering - Extract reads that are exactly 22nt and start with 'G'
3. Bowtie Alignment - Align filtered reads to genome with unique mapping (-m 1)
4. featureCounts - Generate gene-level count matrix for DESeq2 analysis

Requirements

• Nextflow >= 23.04.0
• Docker or Singularity (recommended)
• OR: cutadapt, gawk, bowtie, samtools, subread installed locally

Quick Start

1. Prepare Input Samplesheet

Create a CSV file with your samples (see samplesheet.csv example):

sample,fastq,strandedness
sample1,/path/to/sample1.umi_extract.fastq.gz,reverse
sample2,/path/to/sample2.umi_extract.fastq.gz,reverse

Strandedness options:

• forward - First-strand synthesis
• reverse - Second-strand synthesis (most common for small RNA-seq)
• unstranded - No strand information

2. Run the Pipeline

Basic usage:

nextflow run main.nf \
  --input samplesheet.csv \
  --outdir results \
  --fasta genome.fa \
  --gtf annotations.gtf

With adapter trimming:

nextflow run main.nf \
  --input samplesheet.csv \
  --outdir results \
  --fasta genome.fa \
  --gtf annotations.gtf \
  --adapter AGATCGGAAGAGC

With Docker:

nextflow run main.nf \
  -profile docker \
  --input samplesheet.csv \
  --outdir results \
  --fasta genome.fa \
  --gtf annotations.gtf

With pre-built Bowtie index:

nextflow run main.nf \
  --input samplesheet.csv \
  --outdir results \
  --fasta genome.fa \
  --gtf annotations.gtf \
  --bowtie_index /path/to/bowtie_index_dir

Parameters

Required Parameters

Parameter	Description
--input	Path to samplesheet CSV file
--outdir	Output directory for results
--fasta	Reference genome FASTA file
--gtf	Gene annotation GTF file

Optional Parameters

Parameter	Default	Description
--adapter	null	Adapter sequence to trim (e.g., 'AGATCGGAAGAGC')
--min_length	22	Minimum read length after trimming
--max_length	22	Maximum read length after trimming
--quality_cutoff	0	Quality trimming threshold (0 = disabled)
--bowtie_index	null	Pre-built Bowtie index directory
--feature_type	exon	Feature type for counting (GTF column 3)
--attribute_type	gene_id	Attribute for gene ID (GTF column 9)
--max_mismatches	2	Maximum mismatches for Bowtie alignment

Output Structure

results/
├── cutadapt/
│   ├── sample1.trimmed.fastq.gz
│   └── sample1.cutadapt.log
├── filtered_fastq/
│   ├── sample1.22g.fastq.gz
│   └── sample1.filter_stats.txt
├── bowtie_index/          # (if built from FASTA)
│   └── genome.*
├── alignments/
│   ├── sample1.bam
│   ├── sample1.bam.bai
│   └── sample1.align_stats.txt
├── counts/
│   ├── reverse.featureCounts.txt         # Count matrix
│   └── reverse.featureCounts.txt.summary # Alignment summary
└── pipeline_info/
    ├── timeline.html
    ├── report.html
    └── trace.txt

Input Data Preparation

Expected Input Format

This pipeline expects UMI-extracted FASTQ files, typically generated with:

umi_tools extract \
  -I raw_reads.fastq.gz \
  -S output.umi_extract.fastq.gz \
  --extract-method=string \
  --bc-pattern=NNNN

After UMI extraction, reads should have the UMI in the read name (e.g., @readname_ACGT).

Why Cutadapt is Needed

Even after UMI extraction, reads may:

• Contain adapter sequences at 3' end
• Vary in length (18-26nt instead of exactly 22nt)
• Have low-quality bases at ends

Cutadapt ensures only high-quality 22nt reads are retained for 22G RNA analysis.

Understanding the Output

Count Matrix (counts/*.featureCounts.txt)

The main output is a tab-delimited count matrix where:

• Rows = genes (from GTF annotation)
• Columns = samples
• Values = number of uniquely mapped 22G reads per gene

This file can be directly imported into DESeq2 for differential expression analysis:

library(DESeq2)
counts <- read.table("reverse.featureCounts.txt", header=TRUE, row.names=1, skip=1)
counts <- counts[, 6:ncol(counts)]  # Remove annotation columns
colnames(counts) <- gsub(".bam", "", colnames(counts))

Statistics Files

• cutadapt/*.cutadapt.log - Trimming statistics (reads in/out, adapter found)
• filtered_fastq/*.filter_stats.txt - 22G filtering stats (total vs. 22G reads)
• alignments/*.align_stats.txt - Bowtie alignment metrics (mapped/unmapped)
• counts/*.summary - featureCounts assignment summary (assigned/unassigned)

Troubleshooting

Low 22G Read Count

If very few reads pass the 22G filter:

1. Check read length distribution in UMI-extracted files
2. Verify reads start with 'G' (biological property of 22G RNAs)
3. Adjust --min_length and --max_length if needed

Low Mapping Rate

If few reads align to genome:

1. Verify genome FASTA matches your organism
2. Check if reads need different alignment parameters
3. Inspect reads in filtered_fastq/ for quality

featureCounts: No Counts

If all reads are "Unassigned_NoFeatures":

1. Verify GTF annotation matches genome coordinates
2. Check --strandedness setting (try 'reverse', 'forward', 'unstranded')
3. Confirm --feature_type and --attribute_type match your GTF format

Citation

If using this pipeline, please cite:

• Nextflow: Di Tommaso et al. (2017) Nat Biotechnol
• Cutadapt: Martin (2011) EMBnet.journal
• Bowtie: Langmead et al. (2009) Genome Biol
• Subread: Liao et al. (2014) Bioinformatics

License

MIT License