# Bio Applied Lncrna Classification

> 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.

- Skill: `pavel-kravchenko/bio-applied-lncrna-classification` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-lncrna-classification`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-lncrna-classification/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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)
| 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.

```bash
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
```

```python
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).

```bash
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
```

```python
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).

```bash
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
```

```python
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.

```r
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`