Bio Rna Structure Structure Probing

Version Compatibility

Reference examples tested with: ShapeMapper2 2.1.5+, ViennaRNA 2.6+, SEISMIC-RNA 0.20+, matplotlib 3.8+, pandas 2.2+, numpy 1.26+

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.

Structure Probing

"Process my SHAPE-MaP experiment to get RNA reactivity profiles" -> Convert per-nucleotide chemical-modification signal into reactivities, normalize them, then feed them as SOFT restraints into thermodynamic folding.

• CLI: shapemapper (ShapeMapper2) for end-to-end SHAPE-MaP / DMS-MaP processing
• CLI: RNAfold --shape (ViennaRNA) for reactivity-restrained folding
• CLI: seismic (SEISMIC-RNA) for DMS-MaPseq and multi-conformation clustering

The governing principle: reactivity is flexibility, not pairing, and it is a constraint, not a structure

SHAPE reagents acylate the ribose 2'-OH at a rate set by local nucleotide FLEXIBILITY / conformational dynamics; DMS methylates the Watson-Crick face of A (N1) and C (N3). High reactivity means flexible/accessible, low means constrained. The routine over-interpretation is "low reactivity = base-paired": a nucleotide can be unreactive because it is base-paired OR because it is tertiary-contacted, protein-bound, ligand-occluded, or stacked. Reactivity probes CONSTRAINT, and it cannot by itself distinguish pairing from other protection. Before interpreting any low-reactivity region, rule out the TECHNICAL causes first: a long flat run can be low read depth / no-data or global undermodification, not structure -- check effective depth and confirm low-depth positions are carried as -999, not 0.

Two consequences govern everything below:

• Reactivities are a RESTRAINT on folding, not a structure. They enter as pseudo-free-energy terms that bias the thermodynamic fold and raise accuracy substantially, but they do not yield a structure on their own and cannot disambiguate pairing from other protection.
• A per-position profile is a POPULATION AVERAGE. If the RNA samples more than one structure (riboswitches, dynamic mRNAs), the averaged profile can match no real structure. Detecting and deconvolving multiple conformations requires per-read data (mutational profiling) and clustering (SEISMIC-RNA / DREEM), not a single profile.

Reagent and readout choices

SHAPE acts on the backbone 2'-OH, so it reports all four bases; DMS reads only A/C by default (G/U carry no Watson-Crick-face signal and must be masked to no-data before folding). Newer DMS-MaPseq protocols recover G(N1)/U(N3) signal for a four-base readout (via mutation-signature filtering plus optimized buffer, not the buffer alone) -- only treat DMS as four-base if the protocol and analysis explicitly enable it; otherwise mask G/U.

Mutational profiling (MaP) vs RT-stop is a fundamental analysis fork, not a detail: in MaP the reverse transcriptase reads THROUGH the adduct (Mn2+/TGIRT/Marathon-RT) and misincorporates, encoding each modification as a point mutation; in RT-stop the adduct truncates the cDNA and the read 5'-end is counted. They need different scoring, and a pipeline tuned for one is wrong for the other.

icSHAPE / Structure-seq / DMS-seq are RT-stop methods; do not push them through a MaP mutation-rate pipeline (ShapeMapper2/SEISMIC).

ShapeMapper2: the three samples and verified flags

The three sample roles are distinct: MODIFIED is the signal; UNTREATED subtracts background (SNPs, RT errors, intrinsic damage); DENATURED normalizes sequence-dependent reactivity bias (divide modified by denatured when no good no-data reference exists). The untreated control is mandatory; the denatured control improves normalization.

# ShapeMapper2 is Linux-only; on macOS use Docker/Singularity (see usage-guide).
shapemapper \
    --target target_rna.fa \
    --name my_rna \
    --modified --R1 mod_R1.fastq.gz --R2 mod_R2.fastq.gz \
    --untreated --R1 unmod_R1.fastq.gz --R2 unmod_R2.fastq.gz \
    --out results/ \
    --nproc 8 \
    --min-depth 5000

ShapeMapper2 writes a results/ tree. The reactivity table is <name>_<RNA>_profile.txt; the folder-ready files are SEPARATE: <name>_<RNA>.shape (2 columns: position, normalized reactivity; excluded = -999) and <name>_<RNA>.map (4 columns: position, normalized reactivity, stderr, base). There is no combined _map.shape file.

Key profile.txt columns: Nucleotide is the POSITION integer (1-based); Sequence is the base character; Reactivity_profile is raw; Norm_profile is after normalization. Fold with Norm_profile (or the .shape/.map file), never the raw Reactivity_profile -- the pseudo-energy parameters assume normalized input.

Normalization: per-transcript, and why raw reactivities are not comparable

Raw reactivity (background-subtracted modified mutation rate) sits on an arbitrary, experiment-specific scale set by reagent dose, RT efficiency, and depth. The standard 2-8% / box-plot normalization excludes outliers (top ~2% as a whisker cap), then scales by the mean of the next most-reactive ~8-10% of nucleotides, so normalized values mostly fall ~0-2 with ~1.0 = average reactivity. This scale factor is PER TRANSCRIPT: raw reactivities from different transcripts or experiments are NOT comparable, so never pool or compare raw reactivities across them. To compare two conditions (e.g. +/- ligand), use delta-SHAPE at matched positions with the per-nt standard errors, not raw subtraction. Low-depth nucleotides must become no-data (-999), not zero.

Reactivity-restrained folding and the pseudo-energy model

The Deigan model adds a soft pseudo-energy to every nucleotide in a stacked pair: deltaG = m * ln(1 + reactivity) + b. It is a restraint, not a hard constraint -- a nucleotide can still pair against the data if the global fold demands it.

The m=1.8, b=-0.6 pair is the Hajdin et al. 2013 standard and the ViennaRNA "Deigan" DEFAULT -- it is NOT Deigan et al. 2009's own values (m=2.6, b=-0.8); cite it correctly. For DMS, apply the restraint ONLY to A/C and set G/U to -999, or the model invents constraints at bases that carry no signal. The folded energy a tool reports after a SHAPE/DMS restraint INCLUDES the pseudo-energy bonus, so it is not comparable to an unrestrained MFE -- judge the result by structure agreement, not by a more-negative energy.

# Fold directly from the ShapeMapper2 .shape file (already normalized)
RNAfold --shape=results/my_rna_my_rna.shape --shapeMethod="Dm1.8b-0.6" --noPS < target_rna.fa

In-cell vs in-vitro: the occupancy trap

In-vitro (refolded, deproteinized) RNA reports pure thermodynamics; in-cell reports the RNA as it exists, with bound proteins, ligands, and chaperone-remodeled states all altering reactivity. An in-cell PROTECTED nucleotide may be protein-bound or ligand-occluded, not base-paired -- the single biggest in-cell misinterpretation. Cells also actively unfold mRNA: genome-wide in-vivo DMS showed mRNAs are MORE unfolded in vivo than in vitro (Rouskin 2014). The in-cell-minus-in-vitro difference is itself the signal for protein/ligand footprints (Spitale 2015).