GPTomics/bio-small-rna-seq-smrna-preprocessing

Version Compatibility

Reference examples tested with: cutadapt 4.4+, fastp 0.23+, seqkit 2.6+, umi_tools 1.1+, matplotlib 3.8+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Small RNA Preprocessing

"Preprocess my small RNA-seq reads" -> Remove the 3' adapter, remove any UMI or 4N degenerate bases, size-select to the target class window, and collapse identical reads to a counted FASTA before quantification or discovery.

• CLI: cutadapt -a ADAPTER -m 18 -M 30 --discard-untrimmed then class-specific UMI/4N handling and seqkit rmdup -s

The governing principle: the insert IS its ends, so adapter handling decides everything

In small RNA-seq the molecule is shorter than the read (insert ~18-32 nt, read 50-75 nt), so the 3' adapter is sequenced through on EVERY real insert. A read with no adapter is therefore not a complete small RNA (the insert was too long, or it is an adapter dimer or junk), which inverts the genomic-DNA intuition: here --discard-untrimmed is the correct default, not an aggressive one. Because the insert is defined by its exact 5' and 3' ends, ligation bias never averages out the way fragmentation does in mRNA-seq: T4 RNA ligase captures some miRNA ends 10-100x more efficiently than others, so absolute, cross-miRNA abundance WITHIN a sample is not trustworthy (only the same miRNA compared ACROSS samples, where the per-sequence bias cancels, is reliable; Giraldez 2018). Preprocessing cannot fix ligation bias, but mishandling adapters, UMIs, or size windows manufactures artifacts on top of it.

The read-length histogram after trimming is the primary QC readout, not an afterthought: a sharp peak at 21-23 nt is a healthy miRNA library; a 26-32 nt peak is piRNA (expected in germline, suspicious in soma/plasma); a broad 30+ nt smear with no 22 nt peak is degradation or tRNA/rRNA-fragment contamination or failed size selection; a spike near insert length 0 is adapter dimer eating flowcell capacity. Read the histogram before trusting any downstream count. This degradation heuristic assumes a standard ligation library: for T4-PNK / PANDORA-seq / phospho-RNA-seq preps (which deliberately capture 5'-OH and cyclic-phosphate tRFs and rRFs) the heuristic INVERTS - broad ~18-35 nt tRF/rRF peaks are the expected signal, not contamination.

Decision: how to handle the 3' end depends on the kit

| Kit | 3' adapter | Degenerate / UMI design | Preprocessing consequence | |-----|-----------|--------------------------|----------------------------| | Illumina TruSeq | TGGAATTCTCGGGTGCCAAGG | invariant ends (high ligation bias) | trim adapter only; do NOT PCR-dedup | | NEBNext | AGATCGGAAGAGCACACGTCT | invariant ends | trim adapter only; do NOT PCR-dedup | | NEXTflex (Bioo/PerkinElmer) | TGGAATTCTCGGGTGCCAAGG | 4 random nt on each adapter end (debiasing spacer) | trim adapter, then STRIP 4 nt from each insert end (-u 4 -u -4); the 4N is NOT a UMI, discard it | | QIAseq miRNA | AACTGTAGGCACCATCAAT | true 12-nt UMI 3' of the adapter | EXTRACT the UMI (keep it), align, then UMI-dedup; never position-dedup | | SMARTer / CATS (template-switching) | no ligation adapter | adds a 3' poly-A/tail, not a 4N or UMI | trim the 3' poly-tail, NOT a ligation adapter; different bias profile; ultra-low input | | RealSeq (circularization) | single adapter, one ligation | sidesteps the two-junction ligation bias | low-input; one-ligation chemistry, not two |

The single most damaging error in this table is conflating the NEXTflex 4N debiasing spacer (only 4^4=256 combinations, must be DISCARDED) with the QIAseq 12-nt UMI (must be KEPT and used to separate PCR duplicates from biological duplicates). Using the 4N as a pseudo-UMI saturates instantly and undercounts abundant miRNAs.

Adapter trimming with cutadapt

# Standard ligation-based small-RNA library (TruSeq adapter shown)
cutadapt \
    -a TGGAATTCTCGGGTGCCAAGG \
    -m 18 \
    -M 30 \
    -q 20 \
    --discard-untrimmed \
    -j 8 \
    -o trimmed.fastq.gz \
    input.fastq.gz

# -a: 3' adapter (cutadapt finds it even when only a prefix is sequenced)
# -m 18 / -M 30: keep the small-RNA window; -m drops adapter dimers (trim to ~0)
# -q 20: light 3' quality trim, applied BEFORE adapter removal (cutadapt orders it internally)
# --discard-untrimmed: a read with no adapter is not a complete small RNA

Class-specific size windows

# miRNA-focused window (mature miRNAs cluster at 21-23 nt)
cutadapt -a TGGAATTCTCGGGTGCCAAGG -m 18 -M 26 --discard-untrimmed -o mirna.fastq.gz input.fastq.gz

# piRNA / tRNA-half window (widen -M; do not clip the very class of interest)
cutadapt -a TGGAATTCTCGGGTGCCAAGG -m 24 -M 35 --discard-untrimmed -o pirna.fastq.gz input.fastq.gz

Removing 4N degenerate bases (NEXTflex / high-definition adapters)

# ORDER MATTERS: trim the adapter FIRST, then strip the 4 random nt from each insert end.
# Stripping a fixed 4 nt before adapter removal would corrupt the adapter search.
cutadapt -a TGGAATTCTCGGGTGCCAAGG -m 18 -M 30 --discard-untrimmed -o adapter_trimmed.fastq.gz input.fastq.gz
cutadapt -u 4 -u -4 -o final.fastq.gz adapter_trimmed.fastq.gz

# -u 4: remove 4 nt from the 5' end; -u -4: remove 4 nt from the 3' end (negative = 3')

Extracting and using a true UMI (QIAseq)

# The 12-nt UMI sits immediately 3' of the QIAGEN adapter. Capture it into the read name,
# align, then collapse reads sharing sequence+position+UMI (PCR duplicates) but keep reads
# that differ in UMI (distinct biological molecules). Position-only dedup is WRONG for small RNA.
umi_tools extract --extract-method=regex \
    --bc-pattern='.+(?P<discard_1>AACTGTAGGCACCATCAAT)(?P<umi_1>.{12}).*' \
    -I input.fastq.gz -S umi_extracted.fastq.gz
# ... adapter-trim, map ...
umi_tools dedup --method=directional -I aligned.bam -S deduped.bam
# directional models 1-edit UMI sequencing errors; raw unique-UMI counting overcounts

Using fastp as an alternative

fastp \
    --in1 input.fastq.gz \
    --out1 trimmed.fastq.gz \
    --adapter_sequence TGGAATTCTCGGGTGCCAAGG \
    --length_required 18 \
    --length_limit 30 \
    --json report.json --html report.html

# --length_limit caps the small-RNA window; do NOT use fastp --dedup on small RNA
# (it is sequence-based and deletes real biological duplicates)

Collapse identical reads to a counted FASTA

# seqkit rmdup -s only DEDUPLICATES identical sequences; it does NOT append the _xN
# count that miRDeep2 needs. Use it to shrink the file, but generate the counted FASTA
# with the awk/Python helper below. For a UMI library, collapse on sequence+UMI (or skip
# collapsing) so the UMI survives dedup.
seqkit rmdup -s trimmed.fastq.gz -o dedup.fasta

import gzip
from collections import Counter

def collapse_reads(fastq_path, lo=18, hi=30):
    counts = Counter()
    with gzip.open(fastq_path, 'rt') as f:
        while True:
            header = f.readline()
            if not header:
                break
            seq = f.readline().strip()
            f.readline()
            f.readline()
            if lo <= len(seq) <= hi:
                counts[seq] += 1
    return counts

def write_collapsed_fasta(counts, output_path):
    # miRDeep2 reads the _xN suffix as the read count; preserve it
    with open(output_path, 'w') as f:
        for i, (seq, count) in enumerate(counts.most_common()):
            f.write(f'>seq_{i}_x{count}\n{seq}\n')

QC and contamination gate with miRTrace

# Run miRTrace BEFORE quantifying. Beyond length/complexity, it reports the RNA-class
# composition (miRNA vs rRNA/tRNA/artifact) AND fingerprints clade-specific miRNAs to
# detect cross-species / reagent / sample-swap contamination (found in >7% of public
# datasets) that a good genome mapping rate hides. It has kit presets via --protocol.
mirtrace qc --species hsa --protocol illumina -o mirtrace_out *.fastq.gz
# Read: a miRNA-dominant composition is healthy; rRNA/tRNA-dominant means poor size
# selection, degraded input, or low real miRNA; a foreign-clade signal flags contamination.

For plasma/serum specifically, hemolysis is the dominant QC: red blood cells are loaded with miR-451a, so even slight hemolysis floods the sample with erythroid miRNAs and corrupts the circulating profile. Flag it with the miR-451a (RBC-enriched, rises with hemolysis) vs miR-23a-3p (hemolysis-insensitive) relationship - an elevated miR-451a fraction (or delta-Cq(miR-23a-3p - miR-451a) > ~7 by qPCR) marks a hemolyzed sample. Exclude or model hemolyzed samples before differential analysis (see differential-mirna). Use exogenous spike-ins (cel-miR-39) for low-biomass technical normalization.

Read the length distribution as QC

import gzip
from collections import Counter
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

def plot_length_distribution(fastq_path, out_png):
    lengths = Counter()
    with gzip.open(fastq_path, 'rt') as f:
        for i, line in enumerate(f):
            if i % 4 == 1:
                lengths[len(line.strip())] += 1
    xs = sorted(lengths)
    plt.bar(xs, [lengths[x] for x in xs])
    plt.axvspan(21, 23, color='green', alpha=0.15)  # healthy miRNA peak
    plt.xlabel('read length (nt)')
    plt.ylabel('count')
    plt.savefig(out_png)

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | Almost nothing maps; reads ~50-75 nt | 3' adapter never removed (wrong sequence or step skipped) | Set the kit's exact adapter; reads must shrink to ~18-30 nt after trimming | | Length histogram peaks at ~8 random nt over the real peak | 4N spacer not stripped after adapter removal | Add cutadapt -u 4 -u -4 as a second pass | | Abundant miRNAs look flat / undercounted | Position-based PCR dedup on non-UMI data, or 4N used as a UMI | Do not dedup without a real UMI; for QIAseq use umi_tools dedup | | Broad 30+ nt smear, no 22 nt peak | Degraded input / tRNA-rRNA fragments / failed size selection | Inspect RNA quality (DV200, not RIN); rerun size selection; expect mostly non-miRNA classes | | Huge spike at insert length ~0 | Adapter dimers (no-insert ligation), common at low input | -m 18 discards them; report the dimer fraction as a library-quality flag | | Cross-sample counts incomparable | Libraries built with different kits/protocols (bias is protocol-specific) | Never merge or compare counts across kits; rebuild with one protocol | | High mapping rate but odd composition / foreign reads | Cross-species or reagent contamination that mapping rate hides | Run miRTrace clade fingerprinting; exclude or investigate contaminated samples | | Good phospho/PANDORA library flagged as "degraded" | Standard length-histogram heuristic applied to a 5'-OH/cP-capture prep | Expect broad ~18-35 nt tRF/rRF peaks for these preps; the heuristic inverts | | Plasma profile dominated by a few miRNAs across all samples | Hemolysis: red-cell miR-451a contamination | Flag with miR-451a:miR-23a-3p; exclude/model hemolyzed samples; spike-in normalize |

Related Skills

• mirdeep2-analysis - Novel miRNA discovery; consumes collapsed reads
• mirge3-analysis - Fast known-miRNA + isomiR quantification; has its own trimming
• trf-pirna-profiling - tRF and piRNA profiling, where wider size windows and 5'-OH/cP end chemistry matter
• read-qc/adapter-trimming - General adapter trimming concepts and tool behavior
• read-qc/umi-processing - UMI extraction and deduplication mechanics

References

• Martin M. 2011. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17:10-12. doi:10.14806/ej.17.1.200
• Chen S, Zhou Y, Chen Y, Gu J. 2018. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:i884-i890. doi:10.1093/bioinformatics/bty560
• Giraldez MD, Spengler RM, Etheridge A, et al. 2018. Comprehensive multi-center assessment of small RNA-seq methods for quantitative miRNA profiling. Nat Biotechnol 36:746-757. doi:10.1038/nbt.4183
• Smith T, Heger A, Sudbery I. 2017. UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res 27:491-499. doi:10.1101/gr.209601.116
• Sorefan K, Pais H, Hall AE, et al. 2012. Reducing ligation bias of small RNAs in libraries for next generation sequencing. Silence 3:4. doi:10.1186/1758-907X-3-4
• Kang W, Eldfjell Y, Fromm B, et al. 2018. miRTrace reveals the organismal origins of microRNA sequencing data. Genome Biol 19:213. doi:10.1186/s13059-018-1588-9
• Shi J, Zhang Y, Tan D, et al. 2021. PANDORA-seq expands the repertoire of regulatory small RNAs by overcoming RNA modifications. Nat Cell Biol 23:424-436. doi:10.1038/s41556-021-00652-7