Bio Applied Mirna Seq Pipeline

miRNA-seq Processing and Analysis

When to Use

• Processing small RNA-seq FASTQ files (18-30 nt inserts, 3' adapter-ligated libraries)
• Quantifying mature miRNA expression against miRBase (mature.fa / GFF3)
• Running differential miRNA expression between conditions (tumor vs normal, treated vs control)
• Predicting miRNA-mRNA targets via seed-region matching (TargetScan-style) or looking up validated targets in miRTarBase
• Distinguishing isomiRs / arm switching (-5p vs -3p) in count data

Version Compatibility

• cutadapt ≥ 4.0, bowtie 1.3.x (short-read aligner, not bowtie2 — needed for ungapped 1-mismatch miRNA mapping)
• Subread/featureCounts ≥ 2.0, miRBase release 22.1 (hg38 coordinates)
• Python: pandas ≥ 2.0, scipy ≥ 1.10, statsmodels ≥ 0.14
• R: DESeq2 ≥ 1.42 (Bioconductor 3.18), on R ≥ 4.3

Prerequisites

• pip install pandas scipy statsmodels (Python DE path); BiocManager::install("DESeq2") (R path)
• cutadapt, bowtie, subread (featureCounts) on PATH
• miRBase mature.fa (or hairpin.fa) and species-specific miRBase GFF3 for the reference genome build
• Familiarity with bio-read-qc-adapter-trimming and bio-differential-expression-deseq2-basics

miRNA Naming Convention

• miR-21 / miR-21-5p: mature guide strand, 5' arm
• miR-21-3p: 3' arm (formerly annotated miR-21*)
• hsa-miR-21-5p: species prefix (human) + arm
• mir-21 (lowercase): the gene/precursor locus, not the mature product
• Seed region = positions 2-7 of the mature miRNA — the primary determinant of target recognition

Library Design

Standard poly-A-selected RNA-seq fails for miRNAs: they lack poly-A tails and are already too short to fragment further. Use 3'-adapter-ligation-based small RNA libraries (e.g. Illumina TruSeq Small RNA) with 50 bp single-end sequencing — sufficient since inserts are 18-25 nt.

Processing Pipeline

Goal: go from raw small-RNA FASTQ to a per-sample miRNA count matrix. Approach: trim the 3' adapter with size selection (16-28 nt keeps mature miRNAs and excludes untrimmed adapter dimers), align with Bowtie in ungapped mode with up to 1 mismatch, then count reads over miRBase mature-miRNA features.

# 1. Trim 3' adapter, keep only 16-28 nt inserts, discard reads with no adapter found
cutadapt \
  -a TGGAATTCTCGGGTGCCAAGG \
  -m 16 -M 28 \
  --discard-untrimmed \
  -j 8 -o trimmed.fastq.gz sample.fastq.gz

# 2. Build a Bowtie index from miRBase mature sequences and align
#    --norc: miRNA libraries are stranded, only match the sense strand
#    -v 1:   allow up to 1 mismatch (captures isomiRs with SNPs/editing)
#    -m 5:   discard reads mapping to >5 loci (multi-mapping precursors)
bowtie-build mature.fa mirbase_index
bowtie -x mirbase_index --norc -v 1 -m 5 -p 8 \
  -q trimmed.fastq.gz -S aligned.sam
samtools sort -o aligned.bam aligned.sam && samtools index aligned.bam

# 3. Quantify: count reads per mature miRNA feature (strand-specific, -s 1)
featureCounts -a hg38_mirbase_v22.gff3 -F GTF \
  -o raw_counts.txt -t miRNA -g Name -s 1 aligned.bam

QC Thresholds

Differential Expression

Goal: find miRNAs whose expression differs between two conditions from a raw count matrix. Approach: for a quick Python-side screen, CPM-normalize, log2-transform, and run a per-miRNA t-test with BH correction; for a publication-grade result, model raw counts directly with DESeq2's negative-binomial GLM (preferred — CPM/t-test ignores mean-variance structure at low counts).

import numpy as np
import pandas as pd
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests


def mirna_de_screen(counts_df: pd.DataFrame, group_a_prefix: str, group_b_prefix: str,
                     fc_threshold: float = 1.0, fdr: float = 0.05) -> pd.DataFrame:
    """Quick CPM + t-test differential expression screen for a miRNA count matrix.

    counts_df: samples as columns, miRNA IDs as index, raw read counts.
    group_a_prefix / group_b_prefix: column-name prefixes identifying each group
        (e.g. "Normal", "Tumor").
    Returns a DataFrame with log2FC, p-value, BH-adjusted p-value, and a
    significance flag. NOTE: for a final analysis prefer DESeq2 (see R block
    below) — it models raw counts with a negative-binomial GLM instead of
    treating log-CPM as normally distributed.
    """
    lib_sizes = counts_df.sum(axis=0)
    cpm = counts_df / lib_sizes * 1e6
    log2_cpm = np.log2(cpm + 1)

    group_a = [c for c in log2_cpm.columns if c.startswith(group_a_prefix)]
    group_b = [c for c in log2_cpm.columns if c.startswith(group_b_prefix)]

    rows = []
    for mirna in log2_cpm.index:
        a_vals = log2_cpm.loc[mirna, group_a].values
        b_vals = log2_cpm.loc[mirna, group_b].values
        log2fc = b_vals.mean() - a_vals.mean()
        _, pval = ttest_ind(b_vals, a_vals)
        rows.append({"miRNA": mirna, "log2FC": log2fc, "pvalue": pval})

    de_df = pd.DataFrame(rows)
    _, padj, _, _ = multipletests(de_df["pvalue"].fillna(1), method="fdr_bh")
    de_df["padj"] = padj
    de_df["significant"] = (de_df["padj"] < fdr) & (de_df["log2FC"].abs() > fc_threshold)
    return de_df.sort_values("padj")

# DESeq2 on the raw miRBase count matrix — preferred over CPM/t-test above.
library(DESeq2)

counts <- read.delim("raw_counts.txt", row.names = 1, comment.char = "#")
coldata <- data.frame(
  row.names = colnames(counts),
  condition = factor(rep(c("Normal", "Tumor"), each = ncol(counts) / 2))
)

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]      # filter low-count miRNAs
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Tumor", "Normal"), alpha = 0.05)
res <- res[order(res$padj), ]
sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
write.csv(as.data.frame(res), "mirna_de_results.csv")

Target Prediction

Goal: find candidate 3'UTR target sites for a mature miRNA's seed sequence. Approach: reverse-complement the seed (positions 2-7), scan the UTR for matches, and classify site strength by whether the flanking bases match the 8mer/7mer-m8/7mer-A1 rules.

def find_seed_sites(mirna_seq: str, utr_seq: str, seed_start: int = 1, seed_end: int = 7) -> list[dict]:
    """Find candidate seed-match sites for a miRNA in a 3'UTR sequence.

    mirna_seq: mature miRNA sequence, 5'->3', RNA alphabet (A/C/G/U).
    utr_seq: 3'UTR sequence (RNA or DNA alphabet).
    seed_start/seed_end: 0-indexed half-open window, default positions 2-7
        (i.e. index 1:7) per the canonical miRNA seed definition.
    Returns a list of {position, site_type, match} dicts, one per hit.
    Simplification: classifies 8mer/7mer variants using only the flanking-A
    heuristic (real TargetScan additionally requires exact Watson-Crick
    pairing at UTR position 8) — good enough for screening, not for
    publication-grade site calls.
    """
    seed = mirna_seq[seed_start:seed_end]
    complement = str.maketrans("ACGUTacgut", "UGCAAugcaa")
    seed_rc_dna = seed.translate(complement)[::-1].replace("U", "T")
    utr_dna = utr_seq.replace("U", "T").upper()

    hits = []
    for i in range(len(utr_dna) - len(seed_rc_dna) + 1):
        if utr_dna[i:i + len(seed_rc_dna)] != seed_rc_dna:
            continue
        has_a_upstream = i > 0 and utr_dna[i - 1] == "A"
        has_match_pos8 = i + len(seed_rc_dna) < len(utr_dna) and utr_dna[i + len(seed_rc_dna)] == "A"
        if has_a_upstream and has_match_pos8:
            site_type = "8mer"
        elif has_match_pos8:
            site_type = "7mer-m8"
        elif has_a_upstream:
            site_type = "7mer-A1"
        else:
            site_type = "6mer"
        hits.append({"position": i, "site_type": site_type, "match": utr_dna[i:i + len(seed_rc_dna)]})
    return hits

Databases

• TargetScan 8.0: context++ score (seed type, site accessibility, conservation)
• miRDB: machine-learning-based MirTarget score
• miRTarBase: experimentally validated targets (CLASH, luciferase reporter assays)

Pitfalls

• Adapter contamination: most reads must contain the adapter; --discard-untrimmed removes reads where none was found (usually contaminants, not real inserts)
• isomiR variation: 5'/3' end heterogeneity is biologically real; decide up front whether to collapse isomiRs into one miRNA count or keep them separate
• Wrong aligner mode: bowtie2 local alignment over-calls multi-mappers for 20-22 nt reads — use bowtie (v1) with -v 1 -m 5, not bowtie2, for mature miRNA mapping
• CPM/t-test on low counts: log-CPM t-tests are anti-conservative at low read depth; use DESeq2/edgeR's NB model for the final DE call
• Multiple testing: hundreds of miRNAs are tested at once — always apply BH/FDR correction, never raw p-values

See Also

• bio-read-qc-adapter-trimming — cutadapt trimming details and QC
• bio-differential-expression-deseq2-basics — DESeq2 model and contrasts in depth
• bio-rna-quantification-featurecounts-counting — featureCounts options and strandedness
• bio-small-rna-seq-target-prediction — dedicated seed-match/TargetScan target prediction skill