Bio Applied Lncrna Classification

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)

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   # produces gffcmp.merged.gtf.tmap with class_code column

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({
        "transcript_id": [f"TCONS_{i+1:06d}" for i in range(n)],
        "class_code": class_codes,
        "length": rng.lognormal(np.log(1200), 0.9, n).astype(int),
        "exon_count": rng.choice(range(1, 12), n,
                                  p=[0.3, 0.25, 0.15, 0.1, 0.07, 0.05, 0.03, 0.02, 0.01, 0.01, 0.01]),
        "tpm_mean": rng.lognormal(np.log(1.5), 1.2, n),
    })
    candidates = filter_lncrna_candidates(transcripts)
    assert candidates["length"].min() >= 200
    assert candidates["class_code"].isin(["u", "x", "i"]).all()
    print(f"{n} assembled -> {len(candidates)} lncRNA candidates")
    print(candidates["class_code"].value_counts().to_string())


if __name__ == "__main__":
    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
# alternative: CPAT (needs prebuilt Human_Hexamer.tsv / Human_logitModel.RData)
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)

#' Build co-expression modules and return module assignments + eigengenes
#' expr_mat: samples x genes matrix (lncRNAs + protein-coding), variance-stabilized counts
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)
  list(modules = module_colors, eigengenes = me_list$eigengenes)
}

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-handling
• bio-rna-quantification-alignment-free-quant
• bio-alternative-splicing-isoform-switching
• bio-small-rna-seq-mirdeep2-analysis