| name | bio-applied-lncrna-classification |
| description | Classify StringTie/gffcompare transcripts into lncRNA subtypes by class code/length/TPM, score coding potential with CPC2/CPAT, detect circRNAs via CIRI2 BSJ reads. Use for lncRNA annotation or circRNA calls. |
| tool_type | python |
| primary_tool | pandas |
Long Non-Coding RNA: Discovery and Classification
When to Use
- Filtering a StringTie + gffcompare assembly down to novel lncRNA candidates (lincRNA/antisense/intronic)
- Deciding whether a novel transcript is protein-coding or non-coding (CPC2, CPAT, PhyloCSF)
- Detecting and quantifying circRNAs from back-splice junction (BSJ) reads with CIRI2
- Building guilt-by-association co-expression modules to infer lncRNA function
- Explaining ncRNA classes (miRNA/siRNA/piRNA/snoRNA/lncRNA/circRNA) or known lncRNA mechanisms (XIST, HOTAIR, NEAT1, MALAT1)
Version Compatibility
StringTie ≥2.2, gffcompare ≥0.12, CPC2 ≥1.0.1 (standalone), CPAT ≥3.0.4, CIRI2 v2.0.6, BWA ≥0.7.17, Python ≥3.10, pandas ≥2.0, numpy ≥1.24, R ≥4.3 with WGCNA ≥1.72.
Prerequisites
pip install pandas numpy matplotlib; StringTie/gffcompare/CIRI2/CPC2 or CPAT installed as CLI tools. Prior familiarity with GTF/GFF (bio-genome-intervals-gtf-gff-handling) and quantification (bio-rna-quantification-alignment-free-quant) is assumed.
ncRNA Landscape (for context)
| Class | Size | Function |
|---|
| miRNA | ~22 nt | mRNA silencing via RISC |
| piRNA | 26-31 nt | Transposon silencing (germline) |
| snoRNA | 60-300 nt | rRNA/tRNA modification |
| lncRNA | >200 nt | Chromatin, splicing, decoy, scaffold |
| circRNA | variable | Back-spliced loop; miRNA sponge |
lncRNA subtypes: lincRNA (intergenic), antisense (opposite strand overlap), intronic, eRNA (enhancer-derived), bidirectional. Known examples: XIST (X-inactivation via PRC2), HOTAIR (bridges PRC2/LSD1), NEAT1 (paraspeckle scaffold), MALAT1 (splicing regulation), H19 (imprinting).
Step 1 — Discover and Filter lncRNA Candidates
Goal: From a StringTie-assembled, gffcompare-classified transcript set, keep only plausible novel lncRNAs.
Approach: Run stringtie per sample, merge with stringtie --merge, then gffcompare -r reference.gtf merged.gtf to assign class codes. Keep codes u (intergenic), x (antisense), i (intronic); require length ≥200 nt, ≥2 exons, and TPM ≥0.1 in ≥2 samples to exclude assembly noise.
stringtie sample.bam -o sample.gtf -p 8
stringtie --merge -G gencode.v44.gtf sample1.gtf sample2.gtf -o merged.gtf
gffcompare -r gencode.v44.gtf -o gffcmp merged.gtf
import numpy as np
import pandas as pd
def filter_lncrna_candidates(transcripts: pd.DataFrame, min_length: int = 200,
min_exons: int = 2, min_tpm: float = 0.1) -> pd.DataFrame:
"""Filter a StringTie/gffcompare .tmap-derived DataFrame down to lncRNA candidates.
`transcripts` must have columns: class_code, length, exon_count, tpm_mean.
Candidate class codes: 'u' (intergenic/lincRNA), 'x' (antisense), 'i' (intronic).
"""
candidate_codes = {"u", "x", "i"}
mask = (
(transcripts["length"] >= min_length)
& (transcripts["exon_count"] >= min_exons)
& (transcripts["class_code"].isin(candidate_codes))
& (transcripts["tpm_mean"] >= min_tpm)
)
return transcripts.loc[mask].copy()
def demo():
"""Simulate a StringTie+gffcompare assembly and apply the lncRNA filters."""
rng = np.random.default_rng(7)
n = 8000
class_codes = rng.choice(
["=", "c", "j", "u", "x", "i", "o", "e", "p"], n,
p=[0.30, 0.15, 0.12, 0.20, 0.08, 0.06, 0.04, 0.03, 0.02],
)
transcripts = pd.DataFrame({
: [ i (n)],
: class_codes,
: rng.lognormal(np.log(), , n).astype(),
: rng.choice((, ), n,
p=[, , , , , , , , , , ]),
: rng.lognormal(np.log(), , n),
})
candidates = filter_lncrna_candidates(transcripts)
candidates[].() >=
candidates[].isin([, , ]).()
()
(candidates[].value_counts().to_string())
__name__ == :
demo()
Step 2 — Coding Potential Assessment
Goal: Confidently distinguish non-coding transcripts from unannotated/short protein-coding genes.
Approach: Extract candidate FASTA (gffread -w candidates.fa -g genome.fa merged.gtf), score with CPC2 or CPAT, then keep only transcripts labeled non-coding below the probability threshold (CPC2 default 0.5; CPAT human threshold 0.364). Cross-check borderline calls with PhyloCSF (negative score = non-coding evolution) or Ribo-seq (absence of 3-nt periodicity).
CPC2.py -i candidates.fa -o cpc2_out
cpat.py -x Human_Hexamer.tsv -d Human_logitModel.RData -g candidates.fa -o cpat_out
def classify_coding_potential(candidates: pd.DataFrame, cpc2_tsv: str,
prob_threshold: float = 0.5) -> pd.DataFrame:
"""Merge CPC2 output onto candidates and keep only confidently non-coding transcripts.
CPC2 output columns: #ID, transcript_length, peptide_length, Fickett_score,
pI, ORF_integrity, coding_probability, label.
"""
cpc2 = pd.read_csv(cpc2_tsv, sep="\t").rename(columns={"#ID": "transcript_id"})
merged = candidates.merge(
cpc2[["transcript_id", "coding_probability", "label"]],
on="transcript_id", how="inner",
)
return merged[(merged["label"] == "noncoding") & (merged["coding_probability"] < prob_threshold)]
Step 3 — circRNA Detection with CIRI2
Goal: Identify circular RNAs from back-splice junction (BSJ) reads and quantify them.
Approach: Align with a splice/chimeric-aware mapper, then run CIRI2 to detect reads spanning the BSJ (5' end of a downstream exon joined to 3' end of an upstream exon).
bwa mem -T 19 hg38.fa sample_R1.fq sample_R2.fq > sample.sam
perl CIRI2.pl -I sample.sam -O sample_ciri.txt -F hg38.fa -A gencode.v44.gtf -T 8
def load_ciri2_circrnas(ciri_txt: str, total_mapped_reads: int, min_junction_reads: int = 2) -> pd.DataFrame:
"""Load CIRI2 output, filter by junction-read support, and compute CPM.
CIRI2 columns include: circRNA_ID, chr, circRNA_start, circRNA_end,
junction_reads, SM_MS_SMS, gene_id, circRNA_type.
"""
ciri = pd.read_csv(ciri_txt, sep="\t")
ciri = ciri[ciri["junction_reads"] >= min_junction_reads].copy()
ciri["cpm"] = ciri["junction_reads"] / total_mapped_reads * 1e6
return ciri.sort_values("cpm", ascending=False)
Step 4 — Guilt-by-Association via Co-expression (R/WGCNA)
Goal: Infer putative lncRNA function from protein-coding genes co-expressed in the same module.
Approach: Build a weighted co-expression network across samples (lncRNAs + coding genes together), cut into modules, then inspect which coding genes share a module with each lncRNA.
library(WGCNA)
guilt_by_association <- function(expr_mat, power = 6, min_module_size = 30) {
adjacency_mat <- adjacency(expr_mat, power = power)
tom <- TOMsimilarity(adjacency_mat)
diss_tom <- 1 - tom
gene_tree <- hclust(as.dist(diss_tom), method = "average")
modules <- cutreeDynamic(dendro = gene_tree, distM = diss_tom,
minClusterSize = min_module_size, deepSplit = 2)
module_colors <- labels2colors(modules)
me_list <- moduleEigengenes(expr_mat, colors module_colors
modules module_colors eigengenes me_listeigengenes
Pitfalls
- Coordinate systems: BED is 0-based half-open; GTF/GFF/VCF are 1-based inclusive — mixing them causes off-by-one errors when intersecting candidate loci.
- Strand ambiguity: gffcompare's
u/x/i codes depend on correct strand annotation in the reference GTF; unstranded RNA-seq libraries make antisense calls unreliable.
- Coding-potential false negatives: CPC2/CPAT miss short-ORF micropeptides and non-canonical translation; treat calls near the threshold as ambiguous and confirm with PhyloCSF or Ribo-seq periodicity, not a single tool.
- circRNA BSJ reads need chimeric-aware alignment: plain
bwa mem/STAR default settings under-detect split reads; use BWA with soft-clip-tolerant parameters or STAR --chimSegmentMin, and require ≥2 independent junction reads before calling a circRNA.
- Multiple testing: apply FDR correction (Benjamini-Hochberg) when scoring thousands of assembled transcripts simultaneously.
See Also
bio-genome-intervals-gtf-gff-handlingbio-rna-quantification-alignment-free-quantbio-alternative-splicing-isoform-switchingbio-small-rna-seq-mirdeep2-analysis
