Bio Core Nucleic Acid Structure

Nucleic Acid Structure

When to Use

• Comparing A-DNA, B-DNA, and Z-DNA helix parameters (groove width, rise, handedness) for a sequence or crystal structure
• Estimating oligo duplex stability / melting temperature with nearest-neighbor thermodynamics instead of naive GC counting
• Parsing or predicting RNA secondary structure in dot-bracket notation (stems, hairpins, bulges, pseudoknots)
• Measuring base-pair or backbone geometry (C1'-C1' distance, rise, sugar pucker) from a DNA/RNA PDB structure
• Screening a promoter/regulatory sequence for TF motifs or Z-DNA-forming alternating purine-pyrimidine tracts

Version Compatibility

• Biopython >= 1.81 (Bio.Seq, Bio.PDB)
• NumPy >= 1.24
• Python >= 3.10
• ViennaRNA (RNA Python module) >= 2.6, optional, for real MFE folding

Prerequisites

• pip install biopython numpy; conda install -c bioconda viennarna for MFE folding
• Familiarity with FASTA/PDB formats (see bio-sequence-io-read-sequences, bio-structural-biology-structure-io)
• Basic thermodynamics (Gibbs free energy) and RNA dot-bracket notation

DNA Helix Forms

B-DNA is the standard physiological form. A-form is adopted by RNA duplexes and RNA-DNA hybrids (2'-OH forces C3'-endo pucker). Z-DNA forms transiently in alternating purine-pyrimidine tracts (e.g. (CG)n) under high salt or negative supercoiling. All four base pairs are distinguishable in the major groove (unique H-bond donor/acceptor pattern); A-T vs T-A and G-C vs C-G are hard to distinguish in the minor groove, which is why most transcription factors read the major groove.

Goal: compute physical dimensions (contour length, number of turns) of a DNA segment in a given helix form. Approach: look up per-bp rise and pitch for the form and scale by base-pair count.

HELIX_PARAMS = {
    'A-DNA': {'bp_per_turn': 11,   'rise': 2.6, 'diameter': 23},
    'B-DNA': {'bp_per_turn': 10.5, 'rise': 3.4, 'diameter': 20},
    'Z-DNA': {'bp_per_turn': 12,   'rise': 3.7, 'diameter': 18},
}


def helix_dimensions(num_bp, form='B-DNA'):
    """Return contour length (nm), number of helical turns, and diameter (A)
    for num_bp base pairs folded into the given DNA form."""
    p = HELIX_PARAMS[form]
    length_A = num_bp * p['rise']
    return {'form': form, 'length_nm': length_A / 10.0, 'turns': num_bp / p['bp_per_turn']}


# Human genome (~3.2 Gbp, B-DNA): ~2 m of DNA packed into a ~6 um nucleus
d = helix_dimensions(3.2e9, 'B-DNA')
print(f"{d['length_nm']/1e9:.2f} m, {d['turns']/1e6:.0f}M turns")

Nearest-Neighbor Duplex Stability

Stacking interactions between adjacent base pairs contribute more to duplex stability than individual H-bonds, so free energy must be summed over dinucleotide steps, not counted per base.

Goal: estimate duplex free energy (ΔG, kcal/mol at 37°C, 1M NaCl) for primer/probe design or Tm estimation. Approach: sum SantaLucia (1998) unified nearest-neighbor parameters over every dinucleotide step, plus an initiation penalty.

# SantaLucia (1998) unified nearest-neighbor stacking free energies
NN_ENERGIES = {
    'AA': -1.00, 'AT': -0.88, 'AG': -1.28, 'AC': -1.44,
    'TA': -0.58, 'TT': -1.00, 'TG': -1.45, 'TC': -1.30,
    'GA': -1.30, 'GT': -1.44, 'GG': -1.84, 'GC': -2.24,
    'CA': -1.45, 'CT': -1.28, 'CG': -2.17, 'CC': -1.84,
}


def nearest_neighbor_dG(seq):
    """Estimate duplex free energy (kcal/mol) via the nearest-neighbor model.
    More accurate than base counting; initiation penalty is +1.96 kcal/mol."""
    seq = seq.upper()
    dG = 1.96  # initiation penalty
    for i in range(len(seq) - 1):
        dG += NN_ENERGIES.get(seq[i:i + 2], -1.0)
    return dG


for s in ['GCGCGCGC', 'AAAATTTT', 'GATCGATC']:
    print(f"{s}: dG = {nearest_neighbor_dG(s):.2f} kcal/mol")

RNA Secondary Structure: Dot-Bracket Parsing and Folding

Dot-bracket notation: ( pairs with the matching ), . = unpaired. Elements: stem (paired region), hairpin loop (stem capped by unpaired loop, most common), bulge (unpaired bases on one strand), internal loop (mismatches both strands), multi-branch junction (3+ stems meet), pseudoknot (crossing pairs, needs extended []{} alphabet — most tools, including the parser below, cannot handle them).

Goal: parse a dot-bracket structure into base-pair coordinates, and predict a structure from sequence alone when no experimental structure exists. Approach: a stack-based parser recovers pairs in O(n); a Nussinov dynamic-program maximizes pair count (illustrative — not thermodynamically accurate; use ViennaRNA's RNA.fold for real MFE structures).

def parse_dot_bracket(structure):
    """Parse dot-bracket notation into a sorted list of (i, j) base-pair
    indices (0-indexed, i < j). Raises ValueError on unmatched brackets."""
    pairs, stack = [], []
    for i, ch in enumerate(structure):
        if ch == '(':
            stack.append(i)
        elif ch == ')':
            if not stack:
                raise ValueError(f"Unmatched ')' at position {i}")
            pairs.append((stack.pop(), i))
    if stack:
        raise ValueError(f"Unmatched '(' at positions {stack}")
    return sorted(pairs)


def nussinov_fold(sequence, min_loop=3):
    """Simplified Nussinov algorithm: maximize base-pair count subject to a
    minimum hairpin loop size. Returns a dot-bracket string. Illustrates the
    DP approach; not thermodynamically accurate (use ViennaRNA for MFE)."""
    seq = sequence.upper()
    n = len(seq)
    can_pair = {('A', 'U'), ('U', 'A'), ('G', 'C'), ('C', 'G'), ('G', 'U'), ('U', 'G')}
    dp = [[0] * n for _ in range(n)]

    for length in range(min_loop + 2, n + 1):
        for i in range(n - length + 1):
            j = i + length - 1
            dp[i][j] = dp[i][j - 1]  # j unpaired
            for k in range(i, j - min_loop):
                if (seq[k], seq[j]) in can_pair:
                    left = dp[i][k - 1] if k > i else 0
                    inside = dp[k + 1][j - 1] if k + 1 <= j - 1 else 0
                    dp[i][j] = max(dp[i][j], left + inside + 1)

    structure = ['.'] * n

    def traceback(i, j):
        if i >= j or dp[i][j] == 0:
            return
        if dp[i][j] == dp[i][j - 1]:
            traceback(i, j - 1)
            return
        for k in range(i, j - min_loop):
            if (seq[k], seq[j]) in can_pair:
                left = dp[i][k - 1] if k > i else 0
                inside = dp[k + 1][j - 1] if k + 1 <= j - 1 else 0
                if left + inside + 1 == dp[i][j]:
                    structure[k], structure[j] = '(', ')'
                    traceback(i, k - 1)
                    traceback(k + 1, j - 1)
                    return

    traceback(0, n - 1)
    return ''.join(structure)


# GGGAAACCC should fold into a single 3-bp stem + 3-nt loop: (((...)))
predicted = nussinov_fold("GGGAAACCC")
print(predicted, parse_dot_bracket(predicted))

# Prefer real MFE folding when ViennaRNA is available:
try:
    import RNA
    structure, mfe = RNA.fold("GGGAAACCC")
    print(f"MFE structure: {structure} ({mfe:.2f} kcal/mol)")
except ImportError:
    pass  # fall back to nussinov_fold above

Measuring Base-Pair Geometry from a PDB Structure

Goal: verify that a crystal structure adopts the expected helix form by measuring real inter-strand distances. Approach: the Dickerson dodecamer 1BNA (CGCGAATTCGCG) is the classic canonical B-DNA structure; its two chains are antiparallel, so residue i in chain A pairs with residue N-1-i in chain B.

from Bio.PDB import PDBParser, PDBList
import numpy as np


def measure_base_pair_distance(residue1, residue2):
    """C1'-C1' distance (A) between two nucleotides; canonical WC pairs
    fall in ~9.5-11.5 A regardless of helix form."""
    if "C1'" in residue1 and "C1'" in residue2:
        c1_1 = residue1["C1'"].get_vector().get_array()
        c1_2 = residue2["C1'"].get_vector().get_array()
        return np.linalg.norm(c1_1 - c1_2)
    return None


pdbl = PDBList()
pdb_file = pdbl.retrieve_pdb_file('1BNA', pdir='pdb_files', file_format='pdb')
structure = PDBParser(QUIET=True).get_structure('dna', pdb_file)
model = structure[0]
chain_a, chain_b = list(model.get_chains())[:2]
res_a = [r for r in chain_a if r.id[0] == ' ']
res_b = [r for r in chain_b if r.id[0] == ' ']
n = min(len(res_a), len(res_b))
for i in range(n):
    d = measure_base_pair_distance(res_a[i], res_b[n - 1 - i])
    if d is not None:
        print(f"{res_a[i].get_resname():>3s} - {res_b[n-1-i].get_resname():>3s}: {d:.2f} A")

Pitfalls

• Strand polarity: DNA is synthesized 5'->3'; duplex strands are antiparallel. Always search both strands for motifs (Bio.Seq.Seq.reverse_complement).
• RNA MFE is only one possibility: the same RNA can fold multiple ways in vivo. Use ViennaRNA/RNA.fold for MFE; consider suboptimal structures (RNA.subopt) for riboswitches and regulatory RNAs.
• Pseudoknots: parse_dot_bracket above only handles (); pseudoknotted structures need extended []{} brackets and most standard tools cannot represent them.
• CpG islands: define with CpG observed/expected > 0.6 and GC > 55% over > 200 bp. Most CpGs in mammalian genomes are methylated and depleted by 5mC->T deamination, so naive CpG counts underestimate island boundaries.
• Coordinate systems: BED is 0-based half-open; VCF/GFF/PDB residue numbering is 1-based inclusive — off-by-one errors are common when mixing these with structure files.
• Z-DNA scoring: alternating purine-pyrimidine tracts (CG, CA) are necessary but not sufficient; actual B-to-Z transition also depends on salt and supercoiling, so sequence-only scores are propensity, not prediction.

See Also

• bio-structural-biology-structure-io — loading/parsing PDB/mmCIF structures with Bio.PDB
• bio-structural-biology-geometric-analysis — distances, angles, and torsions on parsed structures
• bio-rna-structure-secondary-structure-prediction — production-grade RNA folding (ViennaRNA, MFE, suboptimal ensembles)
• bio-sequence-manipulation-motif-search — IUPAC motif search and TF binding site scanning