Bio Small Rna Seq Trf Pirna Profiling

Version Compatibility

Reference examples tested with: MINTmap 2.0+, unitas 1.7+, SPORTS1.0, proTRAC 2.4+, Python 3.10+ (numpy 1.26+, pandas 2.2+)

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.

tRF and piRNA Profiling

"Profile the tRFs and piRNAs in my small RNA-seq" -> Annotate every small-RNA class, quantify tRNA-derived fragments at locus resolution, and test whether a piRNA population is real and active.

• CLI: MINTmap (tRFs), unitas / SPORTS1.0 (all classes), proTRAC (piRNA clusters)

The governing principle: small RNA-seq is a size cut over 8+ classes, and the kit decides what is captured

A small-RNA library is a ~18-40 nt size selection that pools miRNAs, piRNAs, endo-siRNAs, tRNA-derived fragments, rRNA-derived fragments, snoRNA-derived RNAs, and Y-RNA fragments - treating the output as "a miRNA dataset" is the field's most common error. Two facts dominate the analysis. First, end chemistry decides capture: standard TruSeq-style ligation requires a 5'-monophosphate and a 3'-OH, so 5'-OH and 2',3'-cyclic-phosphate species (angiogenin-cleaved tRNA halves, many tRFs and rRFs) are SILENTLY ABSENT, not lowly expressed - their absence in a TruSeq library is an assay artifact until proven otherwise, and capturing them needs T4 PNK pre-treatment (cP-RNA-seq / PANDORA-seq). Second, detection is not function: any abundant structured RNA sheds breakdown products into the 18-40 nt window, so high read count proves nothing. Functionality must be EARNED by precise reproducible ends, strand bias, phasing, the piRNA ping-pong signature, or AGO/PIWI loading - and database membership (piRBase, MINTbase) is annotation, not proof.

Multimapping is the other defining hazard: tRNA and piRNA loci are highly redundant (many genomic copies), so a read often cannot be assigned to one locus. This is why tRF tools report EXCLUSIVE versus AMBIGUOUS counts, and why piRNAs are quantified at the cluster/family level rather than per sequence.

Two end-chemistry and biogenesis facts change piRNA conclusions specifically. piRNAs (and plant miRNAs) carry a 3' 2'-O-methyl (HENMT1) that suppresses standard 3'-adapter ligation, so they are systematically UNDER-counted - low piRNA yield can be a 3'-end-chemistry artifact, not low abundance. And the ping-pong signature evidences the SECONDARY (slicer-driven, transposon) pathway only: primary piRNAs are PHASED (1U, Zucchini-dependent trail biogenesis), not ping-pong, and adult mammalian testis is >95% pachytene piRNAs that are repeat-depleted and largely non-transposon. A flat ping-pong z-score therefore does NOT mean "no piRNAs" - test phasing as well.

Decision: which tool for which class

tRF quantification with MINTmap

# MINTmap maps trimmed reads against a tRNA-space lookup and emits two tables:
# EXCLUSIVE tRFs (reads that map only within tRNA space) and AMBIGUOUS tRFs.
# Trust exclusive counts; ambiguous reads are shared with non-tRNA loci.
MINTmap -f trimmed.fastq -p sample_out
# Outputs: sample_out-MINTmap_v2-exclusive-tRFs.expression.txt
#          sample_out-MINTmap_v2-ambiguous-tRFs.expression.txt
# tRF type (tRF-5/tRF-3/tRF-1/i-tRF/tRNA-half) and the source tRNA are reported per row.

The tRF subtype carries a biogenesis tell: tRF-1 comes from the pre-tRNA 3' trailer (RNase Z/ELAC2, ending at the Pol III terminator), tRF-3 includes the post-transcriptional CCA (a marker of mature-tRNA origin), and tRNA halves are angiogenin-cleaved and stress-induced. Mitochondrially-encoded tRFs (mse-tRFs) are lost or misassigned if reads are mapped only to the nuclear genome.

All-class annotation with unitas

# Hierarchical annotation: each read assigned to the first matching class.
# Reading the class composition is the first interpretation step.
unitas -input trimmed.fastq -species human
# Output: a UNITAS folder with per-class read fractions (miRNA / piRNA / tRF / rRF / snoRNA / ...)

piRNA cluster detection with proTRAC

# Map reads (e.g. with sRNAmapper/bowtie), then call clusters probabilistically.
proTRAC_2.4.4.pl -genome genome.fa -map reads.map -format SAM
# A real primary-piRNA cluster shows strand asymmetry, 1U bias, and phased 3' ends.

Test the ping-pong signature (is this an active piRNA pathway?)

Goal: Decide whether a putative piRNA population shows the slicer-driven ping-pong amplification signature.

Approach: For sense/antisense read pairs, count 5'-5' overlaps; an active pathway shows a sharp excess at exactly 10 nt (with 1U on primary and 10A on secondary piRNAs).

import numpy as np
from collections import defaultdict

def ping_pong_zscore(plus_5p, minus_5p, max_overlap=30):
    # plus_5p / minus_5p: dict mapping genomic 5' coordinate -> read count, per strand.
    # A sense read at position i and an antisense read whose 5' end sits at i+overlap-1
    # overlap by 'overlap' nt at their 5' ends. Score the overlap histogram; a 10-nt
    # spike (z >> 0) is the ping-pong signature, evidence of an active piRNA pathway.
    hist = np.zeros(max_overlap + 1)
    for pos, n in plus_5p.items():
        for overlap in range(1, max_overlap + 1):
            partner = pos + overlap - 1
            if partner in minus_5p:
                hist[overlap] += n * minus_5p[partner]
    others = np.concatenate([hist[1:10], hist[11:]])
    z10 = (hist[10] - others.mean()) / (others.std() + 1e-9)
    return hist, z10


def phasing_zscore(same_strand_5p, period=27, max_dist=60):
    # Primary piRNAs are produced head-to-tail, so adjacent SAME-strand 5' ends are
    # spaced ~one piRNA length apart. Score the 5'-to-5' distance histogram: a peak at
    # the modal piRNA length (~26-28 nt) is the phasing signal (the primary-pathway
    # complement to ping-pong; proTRAC reports it natively). Test BOTH, not just ping-pong.
    pos = sorted(same_strand_5p)
    hist = np.zeros(max_dist + 1)
    for a in pos:
        for d in range(1, max_dist + 1):
            if (a + d) in same_strand_5p:
                hist[d] += same_strand_5p[a] * same_strand_5p[a + d]
    others = np.delete(hist[1:], period - 1)
    zp = (hist[period] - others.mean()) / (others.std() + 1e-9)
    return hist, zp