Building a Small RNA-Seq Analysis Pipeline You Can Audit

Small RNA sequencing is straightforward to run and easy to misinterpret. This guide walks through a complete pipeline and explains the reasoning behind each choice. By the end you will have five things for each sample:

• a trimmed and collapsed FASTA of unique inserts with read counts
• a hierarchical annotation that assigns every read to miRNA, tRNA fragment, piRNA, snoRNA, rRNA, or unmapped
• a miRBase- or MirGeneDB-keyed count matrix
• an isomiR-level GFF3 from mirtop
• a DESeq2 result table with size factors you chose deliberately

The hardware requirements are modest. A Linux machine with 16 GB RAM is enough, because small RNA alignment is cheap and the genomes are the only large thing involved. You also need conda or a container runtime, the FASTQ files, and one piece of metadata that people routinely fail to record: which library prep kit was used. Everything below assumes single-end reads of 50–75 cycles from human samples.

One caveat comes before anything else. Nothing here is a clinical test. Small RNA levels in blood shift with hemolysis, platelet count, and time of day. They also shift with exercise and sample handling. Interpreting any of it against your health is a conversation with a clinician rather than a script.

1. Identify the library chemistry before you touch a read

The first task is to establish which chemistry produced the data, because everything downstream depends on it. Specifically, it depends on the 3’ adapter sequence and on whether the chemistry includes a unique molecular identifier (UMI). A UMI is a short random barcode attached to each original molecule before amplification. Get the adapter wrong and you will either trim nothing at all or trim into the insert itself.

The common kits each put a different sequence on the 3’ end of a read:

• Illumina TruSeq Small RNA: 3’ adapter TGGAATTCTCGGGTGCCAAGG
• NEBNext Small RNA: 3’ adapter AGATCGGAAGAGCACACGTCT
• QIAseq miRNA UMI: 3’ adapter AACTGTAGGCACCATCAAT, followed by a 12 nt UMI
• NEXTflex (4N randomized ends): TruSeq adapter, plus 4 random nucleotides on each side of the insert that must be clipped after trimming

If the kit was not recorded, you can infer it from the data itself. Take 100k reads and count the most common 12-mers starting at position 25, which sits just past the typical insert:

zcat sample.fastq.gz | head -400000 | awk 'NR%4==2 {print substr($0,25,12)}' \
  | sort | uniq -c | sort -rn | head

A real small RNA library will show one 12-mer at 20–60% of reads, and that sequence is your adapter’s 5’ end. If the top 12-mer appears at position 1 rather than 25, you are looking at adapter dimer. That is the single most common way these libraries fail.

2. Raw QC that tells you something

Standard QC reports are worth generating, but only part of what they show is informative for this assay. FastQC is a reasonable choice, though the per-base quality plots are not what matters here. The two numbers that matter are the insert length distribution after trimming and the adapter-dimer fraction.

fastqc -t 8 -o qc/ *.fastq.gz
multiqc qc/ -o qc/

Go straight to the “Adapter Content” plot. A healthy miRNA library has adapter starting around cycle 22–24 in most reads, because the modal insert is 22 nt. Adapter starting at cycle 1–3 is dimer. Adapter that never appears means one of two things. Either the inserts were long, in which case you sequenced something other than small RNA, or the adapter sequence you supplied is wrong.

3. Trim, filter by length, and handle UMIs

Trimming is the step where most of the signal is either preserved or lost, so the parameters deserve some thought. Use cutadapt rather than Trimmomatic. Small RNA trimming is a 3’ adapter problem with short inserts, and cutadapt’s error model and its --discard-untrimmed flag are built for exactly that. If a read contains no adapter at all, the insert was longer than the read length and therefore cannot be a mature miRNA.

For TruSeq or NEBNext libraries, the call looks like this:

cutadapt \
  -a TGGAATTCTCGGGTGCCAAGG \
  -e 0.1 -O 7 \
  -m 16 -M 45 \
  --max-n 0 \
  -q 20 \
  --discard-untrimmed \
  -j 8 \
  -o trimmed/sample.fq.gz \
  --json=trimmed/sample.cutadapt.json \
  raw/sample.fastq.gz

A few of those settings carry most of the weight. -O 7 requires a 7 nt adapter overlap, which cuts spurious trimming of reads that happen to end in a chance match. -m 16 keeps tRNA halves’ shorter cousins out while retaining the miRNA range. You can drop it to 15 if you care about the shortest tRFs. -M 45 is generous enough for tRNA halves at 30–35 nt and for Y RNA fragments.

QIAseq libraries need two steps, because the UMI sits 3’ of the adapter and has to be moved out of the read before alignment:

cutadapt -a AACTGTAGGCACCATCAAT -e 0.1 -O 7 -m 28 --discard-untrimmed \
  -j 8 -o tmp/sample.umi.fq.gz raw/sample.fastq.gz

umi_tools extract \
  --extract-method=regex \
  --bc-pattern='.+(?P<umi_1>.{12})$' \
  -I tmp/sample.umi.fq.gz -S trimmed/sample.fq.gz

Note the -m 28 in the first call. That is the 16 nt minimum insert plus the 12 nt UMI. After extraction you deduplicate post-alignment with umi_tools dedup --method=unique. UMIs matter more here than in mRNA-seq, because the input is often picograms of RNA and the PCR cycle count is high. If you have them, use them. If you do not, accept that your counts carry amplification bias and do not chase two-fold differences.

For NEXTflex 4N libraries, add -u 4 -u -4 after adapter trimming to clip the randomized ends.

Whatever the chemistry, record the trimming statistics. The percentage of reads passing, together with the length histogram, is your first real QC gate:

zcat trimmed/sample.fq.gz | awk 'NR%4==2 {print length($0)}' \
  | sort -n | uniq -c > qc/sample.lengths.txt

The shape of that histogram tells you what you have. A plasma or serum library should show a peak at 22 nt. A tissue library shows the same peak plus a shoulder at 30–34 nt from tRNA halves. A flat distribution with no 22 nt peak points to degraded RNA or a failed size selection.

4. Collapse to unique sequences

Small RNA data is massively redundant, with the same handful of sequences observed thousands of times. Collapsing identical reads before alignment cuts runtime by an order of magnitude and makes the intermediate files readable.

seqcluster collapse -f trimmed/sample.fq.gz -m 1 -o collapsed/

The same thing can be done without seqcluster:

seqkit fq2fa trimmed/sample.fq.gz | seqkit rmdup -s -D dups.txt > collapsed/sample.fa

We prefer seqcluster collapse for two reasons. It writes the count into the FASTA header in the format that mirtop and the isomiRs R package expect. It also keeps a quality summary for each unique sequence.

5. Annotate hierarchically rather than with a single alignment

Annotation is where quick pipelines most often go wrong, so it is worth structuring deliberately. The common mistake is to align against miRBase hairpins and call the job done. In plasma, miRNAs are frequently a minority of the library. The rest is made up of rRNA and tRNA fragments, Y RNA fragments, and degradation products of mRNA. If you do not count them you cannot tell a real biological shift from a change in the degradation background. Extracellular RNA libraries in particular carry a heavy non-miRNA load and are sensitive to isolation method, which is the central reproducibility problem in that field ¹.

The remedy is to align sequentially with Bowtie 1. These reads are 16–45 nt, and Bowtie 1’s end-to-end short-read mode is the right tool for that length range. Bowtie 2 and STAR are built for longer reads.

# 1. contaminants first: rRNA, tRNA, snRNA, snoRNA, Mt_rRNA from Ensembl ncRNA fasta
bowtie -f -v 1 -k 1 --best --strata --norc \
  -p 8 --un unmapped/step1.fa \
  -x idx/contam collapsed/sample.fa \
  -S align/contam.sam 2> logs/contam.log

# 2. mature miRNAs
bowtie -f -v 1 -k 20 --best --strata --norc \
  -p 8 --un unmapped/step2.fa \
  -x idx/mature_hsa unmapped/step1.fa \
  -S align/mirna.sam 2> logs/mirna.log

# 3. whole genome, everything left
bowtie -f -v 1 -k 50 --best --strata \
  -p 8 --un unmapped/final.fa \
  -x idx/GRCh38 unmapped/step2.fa \
  -S align/genome.sam 2> logs/genome.log

The flags encode specific decisions. -v 1 allows one mismatch across the whole read, which accommodates the single 3’ non-templated addition that is common in miRNAs. It does so without letting genuinely different sequences collapse together. --norc appears on the miRNA step because mature miRNA sequences are already strand-resolved. -k 50 --best --strata on the genome step reflects the fact that miRNAs live in repeated families, and you want the multi-mapping structure visible rather than silently collapsed to one locus.

Build the contaminant index from the Ensembl ncrna.fa filtered to rRNA, tRNA, snRNA, snoRNA, and misc_RNA. Add GtRNAdb mature tRNAs with the CCA tail included. Keep the per-class counts as you go. The ratio of miRNA reads to rRNA fragment reads is one of the more useful sample QC metrics you will have. It is stable within a person and across draws when handling is consistent.

6. Quantify miRNAs and isomiRs

There are two useful ways to turn alignments into numbers, and we run both because they answer different questions.

For fast counts, use miRDeep2’s quantifier.pl against miRBase hairpin and mature FASTAs.

quantifier.pl -p hairpin_hsa.fa -m mature_hsa.fa \
  -r collapsed/sample.fa -t hsa -y sample -d -W

For correct isomiR handling, use mirtop. It converts an alignment against hairpins into a GFF3 with explicit 5’ and 3’ trimming, non-templated additions, and internal SNVs annotated for each variant.

mirtop gff --sps hsa --hairpin hairpin_hsa.fa --gtf hsa.gff3 \
  -o mirtop/ align/mirna_hairpin.bam

mirtop counts --gff mirtop/mirtop.gff --hairpin hairpin_hsa.fa \
  --gtf hsa.gff3 --sps hsa -o mirtop/

For the annotation itself, we recommend MirGeneDB rather than the full miRBase set as your primary reference. miRBase contains a long tail of entries that do not satisfy the structural criteria for a miRNA. Including them inflates your multiple testing burden with loci that are mostly degradation noise. Report against miRBase IDs as well, because that is what the literature uses.

The isomiR layer is worth the extra step. 3’ non-templated uridylation and adenylation shift with cellular state and with handling. Treating a miRNA as a single number can therefore hide a real change in the isomiR distribution. Integrated small RNA, mRNA, and protein studies have shown that miRNA-level regulatory structure only becomes interpretable when the small RNA layer is quantified carefully and then joined to target expression ².

7. Everything that is not a miRNA

The non-miRNA fraction deserves its own pass rather than being treated as leftover. Run unitas or sports1.1 on the trimmed reads as a parallel track. Both classify against tRNA, rRNA, snoRNA, piRNA, and mRNA-derived fragments, and both give you a class composition table for each sample. For tRNA fragments specifically, MINTmap gives exact and unambiguously assigned tRF counts. It also makes the distinction between exclusive and ambiguous tRF sequences explicit.

The unmapped fraction is worth keeping too. If you want to know what is in it, assemble and BLAST it rather than discarding it. Small RNA sequencing of a biofluid captures non-host sequence. The same read population that yields host miRNAs can be used to reconstruct viral genomes and characterize the small RNA response to them ³. Cross-kingdom small RNA signal is real enough that at least some studies now treat it as a first-class layer rather than contamination ⁴.

8. Normalization and differential expression

With a count matrix in hand, the analysis moves into R, and the main risk shifts from alignment to normalization. Load the count matrix into DESeq2 and filter first, keeping features with at least 10 counts in at least the smaller group size.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(cts, coldata, ~ sex + condition)
keep <- rowSums(counts(dds) >= 10) >= 4
dds <- dds[keep, ]
dds <- estimateSizeFactors(dds, type = "poscounts")
dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "condition_B_vs_A", type = "apeglm")

The one non-default choice there is type = "poscounts". It is needed because small RNA matrices are sparse, and the standard median-of-ratios estimator fails when too many features have a zero in any sample.

The composition problem is worse here than in mRNA-seq. In plasma, a handful of miRNAs (miR-486-5p, miR-451a, miR-16-5p, miR-92a-3p) can take 40–70% of all miRNA reads. If one of them moves, median-of-ratios normalization moves everything else in the opposite direction. There are two defenses and both are cheap. Recompute size factors after excluding the top 5 features and check that your calls survive. Also report CPM within the miRNA-mapped subset alongside the DESeq2 output. Pipelines that bundle normalization, differential expression, and target prediction into one run make this easy to skip. That is precisely why you should look at the size factors by hand ⁵.

Include sex as a covariate by default, as the model formula above does. Sex-dependent effects in omics data are pervasive and systematically under-modeled. Tooling now exists specifically to test for them rather than regress them away ⁶.

9. If you have matched mRNA and protein

The reason to run small RNA-seq at all, rather than a targeted panel, is that it joins to the rest of your molecular data. The standard join is anticorrelation between a miRNA and its predicted targets. Use TargetScan conserved-site predictions plus miRDB and require agreement between them. Then test for enrichment of your differential miRNAs’ target sets among down-regulated transcripts, rather than scoring pairs one at a time.

This kind of design has a track record. Multi-omic mouse work has shown small RNA layers integrating cleanly with mRNA and protein when the same tissue and timepoint are profiled together ⁷. In developing human heart, the same design recovered a coherent miRNA network across all three layers rather than a list of individually significant features ².

10. Reuse public data as your control cohort

A single subject gives you no sense of what is typical, even with many timepoints. Public miRNA-seq from GEO and SRA supplies a distributional reference. Reanalysis of public small RNA data through a uniform pipeline is a well-established way to get there ⁸.

The important detail is where you join the two datasets. Reprocess the raw FASTQs through your own pipeline, and never merge someone else’s count matrix with yours. The adapter, the trimming parameters, and the annotation version all shift the numbers more than most biological effects do.

Common problems

Most failures in this assay come from a short list of recurring causes, and nearly all of them are easier to catch early than to correct later.

Adapter dimer dominates the library. The top sequence is 20–40 nt of pure adapter, and you lose most of your depth. Nothing downstream fixes it. Check the fraction with a grep for the adapter at position 1. If it is over 30%, resequence with a better size selection rather than analyzing what is left.

Hemolysis. Red cell lysis floods plasma with miR-451a and miR-486-5p and makes every ratio meaningless. The established qPCR metric is the Cq difference between miR-23a-3p and miR-451a, with a gap above about 5 cycles flagged as hemolyzed. In sequencing data, the practical proxy is the miR-451a fraction of miRNA-mapped reads. Track it for each sample and flag outliers relative to your own baseline. Exclude affected samples rather than correcting them.

Platelet contamination. Slow or single-spin plasma processing carries platelets, which contribute their own miRNA cargo. Use a double spin (1,900 g for 10 min, then 16,000 g for 10 min) and record that you did.

Ligation bias. T4 RNA ligase strongly prefers certain 5’ and 3’ end sequences. Absolute abundances between different miRNAs within a sample are therefore not comparable. Comparisons across samples for the same miRNA are fine. Do not read “miR-X is 50× more abundant than miR-Y” out of a count table.

Zeros that are not zeros. A miRNA at 3 counts in one sample and 0 in another is sampling noise at typical plasma depths. Set a detection floor and hold to it.

Annotation drift. miRBase 21 and miRBase 22 differ in both sequence and naming for a nontrivial number of entries. Pin the version in your pipeline config and record it in the output filenames.

GUI pipelines that hide parameters. Point-and-click frontends over standard tools are useful for getting a first result and for people who do not want to write shell ⁹. The cost is that the trimming overlap, the mismatch allowance, and the normalization method are the parameters that determine your answer. You should be able to state all three from memory for your own runs.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Rebecca T. Miceli, Tzu‐Yi Chen, Yohei Nose, et al. Extracellular vesicles, RNA sequencing, and bioinformatic analyses: Challenges, solutions, and recommendations. Journal of Extracellular Vesicles, 2024. https://doi.org/10.1002/jev2.70005 ↩

2. Adar Aharon-Yariv, Yaxu Wang, Abdalla Ahmed, et al. Integrated small RNA, mRNA and protein omics reveal a miRNA network orchestrating metabolic maturation of the developing human heart. BMC Genomics, 2023. https://doi.org/10.1186/s12864-023-09801-8 ↩ ↩²

3. Mikhail M. Pooggin. Small RNA-Omics for Plant Virus Identification, Virome Reconstruction, and Antiviral Defense Characterization. Frontiers in Microbiology, 2018. https://doi.org/10.3389/fmicb.2018.02779 ↩

4. Gholamhossein Badeli, Kami Kaboosi, Alireza Mohebbi, et al. Cross-Kingdom Multi-Omics Harmonization Uncovers Coordinated Host Defense and Vector Small RNA Regulatory Networks in Begomovirus Transmission. 2026. https://doi.org/10.64898/2026.08.14.744792 ↩

5. Giorgio Giurato, Maria Rosaria De Filippo, Antonio Rinaldi, et al. iMir: An integrated pipeline for high-throughput analysis of small non-coding RNA data obtained by smallRNA-Seq. BMC Bioinformatics, 2013. https://doi.org/10.1186/1471-2105-14-362 ↩

6. Sophie Le Bars, Mohamed Soudy, Enrico Glaab. XYomics: detecting sex-dependent molecular mechanisms in omics data. Nucleic Acids Research, 2026. https://doi.org/10.1093/nar/gkag759 ↩

7. Mira Pavkovic, Lorena Pantano, Cory V. Gerlach, et al. Multi omics analysis of fibrotic kidneys in two mouse models. Scientific Data, 2019. https://doi.org/10.1038/s41597-019-0095-5 ↩

8. Jun-ichi Satoh, Yoshihiro Kino, Shumpei Niida. MicroRNA-Seq Data Analysis Pipeline to Identify Blood Biomarkers for Alzheimer’s Disease from Public Data. Biomarker Insights, 2015. https://doi.org/10.4137/bmi.s25132 ↩

9. Anjana Anilkumar Sithara, Devi Priyanka Maripuri, Keerthika Moorthy, et al. iCOMIC: a graphical interface-driven bioinformatics pipeline for analyzing cancer omics data. NAR Genomics and Bioinformatics, 2022. https://doi.org/10.1093/nargab/lqac053 ↩