Bio Rna Structure Ncrna Search

Version Compatibility

Reference examples tested with: Infernal 1.1.4+, BioPython 1.83+, 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.

ncRNA Search

"Search my sequences for known non-coding RNA families" -> Score candidates against Rfam covariance models that capture both sequence and consensus secondary structure, or build a custom model for a novel family.

• CLI: cmscan for querying sequences against the Rfam CM database
• CLI: cmsearch for one CM against a sequence database
• CLI: cmbuild + cmcalibrate + cmpress for a custom covariance model

The governing principle: a covariance model is only worth it for STRUCTURED RNA

A covariance model (CM) is a profile stochastic context-free grammar that models a family's consensus secondary structure AND primary sequence jointly. Its power source is covariation: a base pair is scored by a joint distribution over pair states, so a compensatory double mutation (C-G to G-C, or a G-U wobble) that PRESERVES the pair scores well even though both positions changed -- a BLAST or profile-HMM sees two mismatches and loses the signal. This is why Infernal detects remote homologs of structured RNAs (tRNA, rRNA, riboswitches, SRP, RNase P, ribozymes, snoRNA) past the sequence twilight zone.

The decisive corollary: a CM offers NO advantage when there is little conserved secondary structure to exploit. Many lncRNAs (R-scape finds no significant covariation in HOTAIR/Xist/SRA), mature miRNAs (~22 nt, no base-pairing in the mature strand), and primary-sequence-only motifs gain nothing from a CM -- it is slow and calibration-heavy, and a profile-HMM (nhmmer) or BLASTN is the correct, faster tool. A CM built from a structure-FREE alignment (no real #=GC SS_cons pairs) collapses to an HMM and buys nothing but runtime.

Infernal toolchain

• cmbuild model.cm aln.sto -- build a CM from a structure-annotated Stockholm alignment; REQUIRES a #=GC SS_cons line.
• cmcalibrate model.cm -- fit E-value statistics; SLOW (minutes to hours) but MANDATORY before any E-value is meaningful. Without it, search still runs but only bit-score thresholding is valid.
• cmpress model.cm -- build the binary index cmscan requires (cmsearch does not need it).
• cmsearch CM seqdb -- one CM vs a sequence database (e.g. one family vs a genome).
• cmscan CMdb seqs -- query sequences vs a CM database (e.g. all of Rfam.cm); this is the genome/transcript ncRNA-annotation layout.
• cmalign CM seqs -- align/fold hits back to the CM consensus to recover each hit's secondary structure.

Rfam.cm ships PRE-CALIBRATED: run cmpress on it, but do NOT cmcalibrate it. Calibration is only for locally built custom models, and must be re-run after every rebuild.

E-values depend on database size; gathering thresholds do not

A CM E-value scales linearly with the searched database size (Z): the SAME hit gets a different E-value depending on what was searched. cmsearch Z defaults to the target sequence-DB size counted on both strands (total residues x2); cmscan Z defaults to (query length x2 x number of models). So a cmscan E-value and a cmsearch E-value for the same locus are NOT comparable, and a flat -E 1e-5 is exactly what curated thresholds replace. Fix Z explicitly with -Z <Mb> for reproducible cross-run comparison.

Bit scores are DB-size-INDEPENDENT, which is precisely why Rfam stores its curated cutoffs as bit scores:

--cut_ga beats a flat E-value for Rfam because each family has a different signal-to-noise profile (a 70 nt tRNA vs a 2900 nt rRNA vs a short riboswitch); one flat cutoff over- or under-calls per family, and reintroduces the DB-size dependence GA was designed to avoid.

The canonical Rfam annotation command

# One-time setup
wget https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz && gunzip Rfam.cm.gz
wget https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.clanin
cmpress Rfam.cm   # NOT cmcalibrate -- Rfam.cm is pre-calibrated

# Annotate a genome. -Z = 2 x genome size in Mb keeps E-values reproducible; --rfam is the
# large-DB strict filter; --nohmmonly forces CM mode so GA cutoffs stay valid for every model.
cmscan -Z 100 --cut_ga --rfam --nohmmonly --fmt 2 --clanin Rfam.clanin \
    --tblout genome.tblout Rfam.cm genome.fa > genome.cmscan

# Clan deoverlapping: --fmt 2 adds the 'olp' column; the documented Rfam filter drops hits
# marked '=' (dominated by a higher-scoring clanmate), keeping '^' (best of an overlap) and '*' (no overlap).
grep -v ' = ' genome.tblout > genome.deoverlapped.tblout

--clanin with --fmt 2 plus the grep -v ' = ' post-filter is the documented deoverlap path; --oclan is a valid in-tool alternative but not what the modern Rfam pipeline uses, so do not treat it as mandatory.

Reading --fmt 2 output (the column shift that breaks fmt-1 parsers)

--fmt 2 prepends an idx column and inserts a clan name column versus the default --fmt 1, so every downstream field index shifts -- a parser written for fmt 1 silently reads the wrong columns. Verified 0-based fmt-2 cmscan indices: idx 0, target/family name 1, accession 2, query/seq name 3, query accession 4, clan 5, mdl type 6, mdl_from 7, mdl_to 8, seq_from 9, seq_to 10, strand 11, trunc 12, pass 13, gc 14, bias 15, score 16, E-value 17, inc 18, olp 19. In cmscan the model name is column 1 and the sequence name is column 3; in cmsearch they are reversed -- a parser must branch on which tool produced the file.

Interpretive columns: trunc is no, 5', 3', or 5'&3' for hits running off a contig end (often REAL incomplete genes worth keeping, not a discard flag; --anytrunc allows truncation at any internal position (E-values become less accurate), --notrunc disables it entirely); bias is the composition correction already subtracted (a large bias flags a low-complexity hit); mdl_from/mdl_to are consensus coordinates, so a small model span is a partial match even when trunc is no.