Annotator

add_common_output_arguments(parser)

Add common output and formatting arguments to an ArgumentParser.

These options control where BPSEQ/CSV/JSON/DOT/PyMOL outputs are written and enable extra CSV exports for stem parameters.

Source code in src/rnapolis/annotator.py

def add_common_output_arguments(parser: argparse.ArgumentParser):
    """
    Add common output and formatting arguments to an ArgumentParser.

    These options control where BPSEQ/CSV/JSON/DOT/PyMOL outputs are written
    and enable extra CSV exports for stem parameters.
    """
    parser.add_argument("-b", "--bpseq", help="(optional) path to output BPSEQ file")
    parser.add_argument("-c", "--csv", help="(optional) path to output CSV file")
    parser.add_argument(
        "-j",
        "--json",
        help="(optional) path to output JSON file",
    )
    parser.add_argument(
        "-e",
        "--extended",
        action="store_true",
        help=(
            "(optional) if set, the program will print extended secondary "
            "structure to the standard output"
        ),
    )
    parser.add_argument("-d", "--dot", help="(optional) path to output DOT file")
    parser.add_argument(
        "-p", "--pml", help="(optional) path to output PyMOL PML script for stems"
    )
    parser.add_argument(
        "--inter-stem-csv",
        help="(optional) path to output CSV file for inter-stem parameters",
    )
    parser.add_argument(
        "--stems-csv",
        help="(optional) path to output CSV file for stem details",
    )
    parser.add_argument(
        "--decompose-pseudoknot-free",
        action="store_true",
        help=(
            "(optional) if set, structural elements (stems, hairpins, loops, "
            "single strands) are decomposed from the pseudoknot-free structure, "
            "but dot-bracket strings retain pseudoknot characters; pseudoknotted "
            "stems are reported separately in the JSON output"
        ),
    )

angle_between_vectors(v1, v2)

Compute the angle in radians between two vectors.

The angle is calculated using the dot product and vector norms.

Parameters:

Name	Type	Description	Default
v1	NDArray[floating]	First vector.	required
v2	NDArray[floating]	Second vector.	required

Returns:

Name	Type	Description
float	float	Angle between the two vectors in radians.

Source code in src/rnapolis/annotator.py

def angle_between_vectors(
    v1: numpy.typing.NDArray[numpy.floating],
    v2: numpy.typing.NDArray[numpy.floating],
) -> float:
    """
    Compute the angle in radians between two vectors.

    The angle is calculated using the dot product and vector norms.

    Args:
        v1: First vector.
        v2: Second vector.

    Returns:
        float: Angle between the two vectors in radians.
    """
    return math.acos(numpy.dot(v1, v2) / numpy.linalg.norm(v1) / numpy.linalg.norm(v2))

detect_bph_br_classification(donor_residue, donor, acceptor)

Classify base–phosphate or base–ribose interactions.

The classification is based on donor and acceptor atom names, residue type and a torsion angle around selected atoms.

Parameters:

Name	Type	Description	Default
donor_residue	Residue3D	Residue that contains the donor atom.	required
donor	Atom	Donor atom in the interaction.	required
acceptor	Atom	Acceptor atom in the interaction.	required

Returns:

Type	Description
Optional[int]	Integer code describing the interaction type, or None if the interaction does not match any known class.

Source code in src/rnapolis/annotator.py

def detect_bph_br_classification(
    donor_residue: Residue3D, donor: Atom, acceptor: Atom
) -> Optional[int]:
    """
    Classify base–phosphate or base–ribose interactions.

    The classification is based on donor and acceptor atom names, residue
    type and a torsion angle around selected atoms.

    Args:
        donor_residue: Residue that contains the donor atom.
        donor: Donor atom in the interaction.
        acceptor: Acceptor atom in the interaction.

    Returns:
        Integer code describing the interaction type, or ``None`` if the interaction does not match any known class.
    """
    # source: Classification and energetics of the base-phosphate interactions
    # in RNA. Craig L. Zirbel, Judit E. Sponer, Jiri Sponer, Jesse Stombaugh
    # and Neocles B. Leontis
    if donor_residue.one_letter_name == "A":
        if donor.name == "C2":
            return 2
        if donor.name == "N6":
            n1 = donor_residue.find_atom("N1")
            c6 = donor_residue.find_atom("C6")
            if n1 is not None and c6 is not None:
                torsion = math.degrees(torsion_angle(n1, c6, donor, acceptor))
                return 6 if -90.0 < torsion < 90.0 else 7
        if donor.name == "C8":
            return 0

    if donor_residue.one_letter_name == "G":
        if donor.name == "N1":
            return 5
        if donor.name == "N2":
            n3 = donor_residue.find_atom("N3")
            c2 = donor_residue.find_atom("C2")
            if n3 is not None and c2 is not None:
                torsion = math.degrees(torsion_angle(n3, c2, donor, acceptor))
                return 1 if -90.0 < torsion < 90.0 else 3
        if donor.name == "C8":
            return 0

    if donor_residue.one_letter_name == "C":
        if donor.name == "N4":
            n3 = donor_residue.find_atom("N3")
            c4 = donor_residue.find_atom("C4")
            if n3 is not None and c4 is not None:
                torsion = math.degrees(torsion_angle(n3, c4, donor, acceptor))
                return 6 if -90.0 < torsion < 90.0 else 7
        if donor.name == "C5":
            return 9
        if donor.name == "C6":
            return 0

    if donor_residue.one_letter_name == "U":
        if donor.name == "N3":
            return 5
        if donor.name == "C5":
            return 9
        if donor.name == "C6":
            return 0

    if donor_residue.one_letter_name == "T":
        if donor.name == "N3":
            return 5
        if donor.name == "C6":
            return 0
        if donor.name == "C7":
            return 9

    return None

detect_cis_trans(residue_i, residue_j)

Classify a base pair as cis or trans based on a torsion angle.

The function uses C1'–N1/N9 atoms of both residues to compute a torsion angle and classifies the pair as cis or trans.

Parameters:

Name	Type	Description	Default
residue_i	Residue3D	First residue in the pair.	required
residue_j	Residue3D	Second residue in the pair.	required

Returns:

Type	Description
Optional[str]	"c" for cis, "t" for trans, or None if the required atoms are missing.

Source code in src/rnapolis/annotator.py

def detect_cis_trans(residue_i: Residue3D, residue_j: Residue3D) -> Optional[str]:
    """
    Classify a base pair as cis or trans based on a torsion angle.

    The function uses C1'–N1/N9 atoms of both residues to compute a torsion
    angle and classifies the pair as cis or trans.

    Args:
        residue_i: First residue in the pair.
        residue_j: Second residue in the pair.

    Returns:
        ``"c"`` for cis, ``"t"`` for trans, or ``None`` if the required atoms are missing.
    """
    c1p_i = residue_i.find_atom("C1'")
    c1p_j = residue_j.find_atom("C1'")

    if residue_i.one_letter_name in "AG":
        n9n1_i = residue_i.find_atom("N9")
    else:
        n9n1_i = residue_i.find_atom("N1")

    if residue_j.one_letter_name in "AG":
        n9n1_j = residue_j.find_atom("N9")
    else:
        n9n1_j = residue_j.find_atom("N1")

    if c1p_i is None or c1p_j is None or n9n1_i is None or n9n1_j is None:
        return None

    torsion = math.degrees(torsion_angle(c1p_i, n9n1_i, n9n1_j, c1p_j))
    return "c" if -90.0 < torsion < 90.0 else "t"

extract_base_interactions(tertiary_structure, model=None)

Extract all base-level interactions from a 3D structure.

This includes base pairs, base–phosphate, base–ribose and stacking interactions and returns them wrapped in a BaseInteractions object.

Parameters:

Name	Type	Description	Default
tertiary_structure	Structure3D	Input 3D structure.	required
model	Optional[int]	Optional model number to analyse.	None

Returns:

Name	Type	Description
BaseInteractions	BaseInteractions	Container with all detected interactions.

Source code in src/rnapolis/annotator.py

def extract_base_interactions(
    tertiary_structure: Structure3D, model: Optional[int] = None
) -> BaseInteractions:
    """
    Extract all base-level interactions from a 3D structure.

    This includes base pairs, base–phosphate, base–ribose and stacking
    interactions and returns them wrapped in a BaseInteractions object.

    Args:
        tertiary_structure: Input 3D structure.
        model: Optional model number to analyse.

    Returns:
        BaseInteractions: Container with all detected interactions.
    """
    base_pairs, base_phosphate, base_ribose = find_pairs(tertiary_structure, model)
    stackings = find_stackings(tertiary_structure, model)
    return BaseInteractions.from_structure3d(
        tertiary_structure, base_pairs, stackings, base_ribose, base_phosphate, []
    )

find_pairs(structure, model=None)

Find base–base, base–phosphate and base–ribose interactions.

The function scans all residues in the structure, builds a KDTree of donors and acceptors and applies geometric filters to detect hydrogen bonds and classify them into interaction types.

Parameters:

Name	Type	Description	Default
structure	Structure3D	3D structure to analyse.	required
model	Optional[int]	Optional model number to restrict the search to a single model.	None

Returns:

Type	Description
Tuple[List[BasePair], List[BasePhosphate], List[BaseRibose]]	Lists of base pairs, base–phosphate and base–ribose interactions.

Source code in src/rnapolis/annotator.py

def find_pairs(
    structure: Structure3D, model: Optional[int] = None
) -> Tuple[List[BasePair], List[BasePhosphate], List[BaseRibose]]:
    """
    Find base–base, base–phosphate and base–ribose interactions.

    The function scans all residues in the structure, builds a KDTree of
    donors and acceptors and applies geometric filters to detect hydrogen
    bonds and classify them into interaction types.

    Args:
        structure: 3D structure to analyse.
        model: Optional model number to restrict the search to a single model.

    Returns:
        Lists of base pairs, base–phosphate and base–ribose interactions.
    """
    # put all donors and acceptors into a KDTree
    coordinates = []
    coordinates_atom_map: Dict[Tuple[float, float, float], Atom] = {}
    coordinates_type_map: Dict[Tuple[float, float, float], str] = {}
    coordinates_residue_map: Dict[Tuple[float, float, float], Residue3D] = {}
    for residue in structure.residues:
        if model is not None and residue.model != model:
            continue
        acceptors = (
            BASE_ACCEPTORS.get(residue.one_letter_name, [])
            + RIBOSE_ACCEPTORS
            + PHOSPHATE_ACCEPTORS
        )
        donors = BASE_DONORS.get(residue.one_letter_name, [])
        for atom_name in acceptors + donors:
            atom = residue.find_atom(atom_name)
            if atom:
                xyz = (atom.x, atom.y, atom.z)
                coordinates.append(xyz)
                coordinates_atom_map[xyz] = atom
                coordinates_type_map[xyz] = (
                    "acceptor" if atom_name in acceptors else "donor"
                )
                coordinates_residue_map[xyz] = residue

    if len(coordinates) < 2:
        return [], [], []

    kdtree = KDTree(coordinates)

    # find all hydrogen bonds
    hydrogen_bonds = []
    base_phosphate_pairs = []
    base_ribose_pairs = []
    used_atoms: Set[Atom] = set()
    for i, j in kdtree.query_pairs(HYDROGEN_BOND_MAX_DISTANCE):
        type_i = coordinates_type_map[coordinates[i]]
        type_j = coordinates_type_map[coordinates[j]]

        # process only acceptor/donor pairs, not acceptor/acceptor or donor/donor
        if type_i == type_j:
            continue

        atom_i = coordinates_atom_map[coordinates[i]]
        atom_j = coordinates_atom_map[coordinates[j]]

        # skip spurious hydrogen bonds in the same residue
        if (
            atom_i.label is not None
            and atom_i.label is not None
            and atom_i.label == atom_j.label
        ):
            continue
        if (
            atom_i.auth is not None
            and atom_i.auth is not None
            and atom_i.auth == atom_j.auth
        ):
            continue

        residue_i = coordinates_residue_map[coordinates[i]]
        residue_j = coordinates_residue_map[coordinates[j]]
        logging.debug(
            f"Checking pair {residue_i.full_name} {atom_i.name} - "
            f"{residue_j.full_name} {atom_j.name}"
        )

        # check for base-phosphate contacts
        if (
            (atom_i.name in PHOSPHATE_ACCEPTORS or atom_j.name in PHOSPHATE_ACCEPTORS)
            and atom_i not in used_atoms
            and atom_j not in used_atoms
        ):
            logging.debug("Checking base-phosphate interaction")
            if type_i == "donor":
                donor_residue, acceptor_residue = residue_i, residue_j
                donor_atom, acceptor_atom = atom_i, atom_j
            else:
                donor_residue, acceptor_residue = residue_j, residue_i
                donor_atom, acceptor_atom = atom_j, atom_i
            bph = detect_bph_br_classification(donor_residue, donor_atom, acceptor_atom)
            if bph is not None:
                used_atoms.add(atom_i)
                used_atoms.add(atom_j)
                base_phosphate_pairs.append((donor_residue, acceptor_residue, bph))
            continue

        # check for base-ribose contacts
        if (
            (atom_i.name in RIBOSE_ACCEPTORS or atom_j.name in RIBOSE_ACCEPTORS)
            and atom_i not in used_atoms
            and atom_j not in used_atoms
        ):
            logging.debug("Checking base-ribose interaction")
            if type_i == "donor":
                donor_residue, acceptor_residue = residue_i, residue_j
                donor_atom, acceptor_atom = atom_i, atom_j
            else:
                donor_residue, acceptor_residue = residue_j, residue_i
                donor_atom, acceptor_atom = atom_j, atom_i
            br = detect_bph_br_classification(donor_residue, donor_atom, acceptor_atom)
            if br is not None:
                used_atoms.add(atom_i)
                used_atoms.add(atom_j)
                base_ribose_pairs.append((donor_residue, acceptor_residue, br))
            continue

        # check for base-base contacts
        if residue_i.base_normal_vector is None or residue_j.base_normal_vector is None:
            continue

        logging.debug("Checking base-base interaction")
        vector = atom_i.coordinates - atom_j.coordinates
        angle1 = math.degrees(
            angle_between_vectors(residue_i.base_normal_vector, vector)
        )
        angle2 = math.degrees(
            angle_between_vectors(residue_j.base_normal_vector, vector)
        )
        logging.debug(
            f"Angles between normals and hydrogen bond: {angle1:.2f} and {angle2:.2f}"
        )
        if (
            HYDROGEN_BOND_ANGLE_RANGE[0] < angle1 < HYDROGEN_BOND_ANGLE_RANGE[1]
            and HYDROGEN_BOND_ANGLE_RANGE[0] < angle2 < HYDROGEN_BOND_ANGLE_RANGE[1]
        ):
            hydrogen_bonds.append((atom_i, atom_j, residue_i, residue_j))

    # match hydrogen bonds with base edges
    labels = []
    for atom_i, atom_j, residue_i, residue_j in hydrogen_bonds:
        edges_i = BASE_EDGES.get(residue_i.one_letter_name, dict()).get(
            atom_i.name, None
        )
        edges_j = BASE_EDGES.get(residue_j.one_letter_name, dict()).get(
            atom_j.name, None
        )
        if edges_i is None or edges_j is None:
            continue

        # detect cis/trans
        cis_trans = detect_cis_trans(residue_i, residue_j)
        if cis_trans is None:
            continue

        logging.debug(
            f"Matched {residue_i.full_name} with {residue_j.full_name} as "
            f"{cis_trans} {edges_i} {edges_j}"
        )

        if residue_i < residue_j:
            for edge_i in edges_i:
                for edge_j in edges_j:
                    labels.append((residue_i, residue_j, cis_trans, edge_i, edge_j))
        else:
            for edge_i in edges_i:
                for edge_j in edges_j:
                    labels.append((residue_j, residue_i, cis_trans, edge_j, edge_i))

    # create a list of base pairs
    base_base_pairs = []
    occupied = set()

    counter = Counter(labels)
    for interaction, hydrogen_bond_count in counter.most_common():
        if hydrogen_bond_count < 2:
            continue

        residue_i, residue_j, cis_trans, edge_i, edge_j = interaction

        if (residue_i, edge_i) in occupied:
            continue
        if (residue_j, edge_j) in occupied:
            continue

        occupied.add((residue_i, edge_i))
        occupied.add((residue_j, edge_j))

        lw = LeontisWesthof[f"{cis_trans}{edge_i}{edge_j}"]
        base_base_pairs.append((residue_i, residue_j, lw))

    base_pairs = []
    for residue_i, residue_j, lw in sorted(base_base_pairs):
        base_pairs.append(
            BasePair(
                Residue(residue_i.label, residue_i.auth),
                Residue(residue_j.label, residue_j.auth),
                lw,
                Saenger.from_leontis_westhof(
                    residue_i.one_letter_name, residue_j.one_letter_name, lw
                ),
            )
        )

    bph_map = merge_and_clean_bph_br(sorted(base_phosphate_pairs))
    base_phosphates = []
    for pair, bphs in bph_map.items():
        residue_i, residue_j = pair
        for bph in bphs:
            base_phosphates.append(
                BasePhosphate(
                    Residue(residue_i.label, residue_i.auth),
                    Residue(residue_j.label, residue_j.auth),
                    BPh[f"_{bph}"],
                )
            )

    br_map = merge_and_clean_bph_br(sorted(base_ribose_pairs))
    base_riboses = []
    for pair, brs in br_map.items():
        residue_i, residue_j = pair
        for br in brs:
            base_riboses.append(
                BaseRibose(
                    Residue(residue_i.label, residue_i.auth),
                    Residue(residue_j.label, residue_j.auth),
                    BR[f"_{br}"],
                )
            )

    return base_pairs, base_phosphates, base_riboses

find_stackings(structure, model=None)

Find base stacking interactions in a 3D structure.

Stacking is detected using geometric centres of bases, distance thresholds and angles between base normals and connecting vectors.

Parameters:

Name	Type	Description	Default
structure	Structure3D	3D structure to analyse.	required
model	Optional[int]	Optional model number to restrict the search.	None

Returns:

Type	Description
List[Stacking]	List of detected stacking interactions.

Source code in src/rnapolis/annotator.py

def find_stackings(
    structure: Structure3D, model: Optional[int] = None
) -> List[Stacking]:
    """
    Find base stacking interactions in a 3D structure.

    Stacking is detected using geometric centres of bases, distance
    thresholds and angles between base normals and connecting vectors.

    Args:
        structure: 3D structure to analyse.
        model: Optional model number to restrict the search.

    Returns:
        List of detected stacking interactions.
    """
    # put all nitrogen ring centers into a KDTree
    coordinates = []
    coordinates_residue_map: Dict[Tuple[float, float, float], Residue3D] = {}
    for residue in structure.residues:
        if model is not None and residue.model != model:
            continue
        base_atoms = BASE_ATOMS.get(residue.one_letter_name, [])
        xs, ys, zs = [], [], []
        for atom_name in base_atoms:
            atom = residue.find_atom(atom_name)
            if atom is not None:
                xs.append(atom.x)
                ys.append(atom.y)
                zs.append(atom.z)
        if len(xs) > 0:
            geometric_center = (sum(xs) / len(xs), sum(ys) / len(ys), sum(zs) / len(zs))
            coordinates.append(geometric_center)
            coordinates_residue_map[geometric_center] = residue

    if len(coordinates) < 2:
        return []

    kdtree = KDTree(coordinates)

    # find all stacking interaction
    pairs = []
    for i, j in kdtree.query_pairs(STACKING_MAX_DISTANCE):
        residue_i = coordinates_residue_map[coordinates[i]]
        residue_j = coordinates_residue_map[coordinates[j]]

        # check angle between normals
        normal_i = residue_i.base_normal_vector
        normal_j = residue_j.base_normal_vector
        if normal_i is None or normal_j is None:
            continue

        angle = min(
            [
                angle_between_vectors(normal_i, normal_j),
                angle_between_vectors(-normal_i, normal_j),
            ]
        )
        if math.degrees(angle) > STACKING_MAX_ANGLE_BETWEEN_NORMALS:
            continue

        vector = numpy.array([coordinates[i][k] - coordinates[j][k] for k in (0, 1, 2)])
        angle = min(
            angle_between_vectors(vector, normal_i),
            angle_between_vectors(vector, normal_j),
        )
        if math.degrees(angle) > STACKING_MAX_ANGLE_BETWEEN_VECTOR_AND_NORMAL:
            continue

        same_direction = True if numpy.dot(normal_i, normal_j) > 0.0 else False

        if residue_i < residue_j:
            if same_direction:
                pairs.append((residue_i, residue_j, "upward"))
            else:
                pairs.append((residue_i, residue_j, "inward"))
        else:
            if same_direction:
                pairs.append((residue_j, residue_i, "downward"))
            else:
                pairs.append((residue_j, residue_i, "outward"))

    stackings = []
    for residue_i, residue_j, topology in sorted(pairs):
        nt1 = Residue(residue_i.label, residue_i.auth)
        nt2 = Residue(residue_j.label, residue_j.auth)
        stackings.append(Stacking(nt1, nt2, StackingTopology[topology]))

    return stackings

generate_pymol_script(mapping, stems)

Generate a PyMOL script to visualise stems as cylinders and dihedrals.

The script creates selections, pseudoatoms on stem centroids and CGO cylinders between them, plus optional dihedral measurements.

Parameters:

Name	Type	Description	Default
mapping	Mapping2D3D	Mapping between 2D indices and 3D residues.	required
stems	List[Stem]	List of stems to visualise.	required

Returns:

Name	Type	Description
str	str	Multi-line PyMOL script.

Source code in src/rnapolis/annotator.py

def generate_pymol_script(mapping: Mapping2D3D, stems: List[Stem]) -> str:
    """
    Generate a PyMOL script to visualise stems as cylinders and dihedrals.

    The script creates selections, pseudoatoms on stem centroids and CGO
    cylinders between them, plus optional dihedral measurements.

    Args:
        mapping: Mapping between 2D indices and 3D residues.
        stems: List of stems to visualise.

    Returns:
        str: Multi-line PyMOL script.
    """
    pymol_commands = []
    radius = 0.5
    r, g, b = 1.0, 0.0, 0.0  # Red color

    for stem_idx, stem in enumerate(stems):
        # Get residues for selection string
        try:
            res5p_first = mapping.bpseq_index_to_residue_map[stem.strand5p.first]
            res5p_last = mapping.bpseq_index_to_residue_map[stem.strand5p.last]
            res3p_first = mapping.bpseq_index_to_residue_map[stem.strand3p.first]
            res3p_last = mapping.bpseq_index_to_residue_map[stem.strand3p.last]

            # Prefer auth chain/number if available
            chain5p = (
                res5p_first.auth.chain if res5p_first.auth else res5p_first.label.chain
            )
            num5p_first = (
                res5p_first.auth.number
                if res5p_first.auth
                else res5p_first.label.number
            )
            num5p_last = (
                res5p_last.auth.number if res5p_last.auth else res5p_last.label.number
            )

            chain3p = (
                res3p_first.auth.chain if res3p_first.auth else res3p_first.label.chain
            )
            num3p_first = (
                res3p_first.auth.number
                if res3p_first.auth
                else res3p_first.label.number
            )
            num3p_last = (
                res3p_last.auth.number if res3p_last.auth else res3p_last.label.number
            )

            # Format selection string: select stem0, A/1-5/ or A/10-15/
            selection_str = (
                f"{chain5p}/{num5p_first}-{num5p_last}/ or "
                f"{chain3p}/{num3p_first}-{num3p_last}/"
            )
            pymol_commands.append(f"select stem{stem_idx}, {selection_str}")

        except (KeyError, AttributeError) as e:
            logging.warning(
                f"Could not generate selection string for stem {stem_idx}: "
                f"Missing residue data ({e})"
            )

        centroids = mapping.get_stem_coordinates(stem)

        # Need at least 2 centroids to draw a segment
        if len(centroids) < 2:
            # Removed warning log for stems with < 2 base pairs
            continue

        # Create pseudoatoms for each centroid
        for centroid_idx, centroid in enumerate(centroids):
            x, y, z = centroid
            pseudoatom_name = f"stem{stem_idx}_centroid{centroid_idx}"
            pymol_commands.append(
                f"pseudoatom {pseudoatom_name}, pos=[{x:.3f}, {y:.3f}, {z:.3f}]"
            )

        # Draw cylinders between consecutive centroids
        for seg_idx in range(len(centroids) - 1):
            p1 = centroids[seg_idx]
            p2 = centroids[seg_idx + 1]
            x1, y1, z1 = p1
            x2, y2, z2 = p2
            # Format: [CYLINDER, x1, y1, z1, x2, y2, z2, radius, r1, g1, b1, r2, g2, b2]
            # Use 9.0 for CYLINDER code
            # Use same color for both ends
            cgo_object = (
                f"[ 9.0, {x1:.3f}, {y1:.3f}, {z1:.3f}, "
                f"{x2:.3f}, {y2:.3f}, {z2:.3f}, {radius}, "
                f"{r}, {g}, {b}, {r}, {g}, {b} ]"
            )
            pymol_commands.append(
                f'cmd.load_cgo({cgo_object}, "stem_{stem_idx}_seg_{seg_idx}")'
            )

        # Calculate and display dihedral angles between consecutive centroids
        if len(centroids) >= 4:
            for i in range(len(centroids) - 3):
                pa1 = f"stem{stem_idx}_centroid{i}"
                pa2 = f"stem{stem_idx}_centroid{i + 1}"
                pa3 = f"stem{stem_idx}_centroid{i + 2}"
                pa4 = f"stem{stem_idx}_centroid{i + 3}"
                dihedral_name = f"stem{stem_idx}_dihedral{i}"
                pymol_commands.append(
                    f"dihedral {dihedral_name}, {pa1}, {pa2}, {pa3}, {pa4}"
                )

    return "\n".join(pymol_commands)

handle_output_arguments(args, structure2d, mapping, input_filename)

Handle all output options based on parsed CLI arguments.

Depending on the flags, this function writes BPSEQ, CSV, JSON, DOT, PyMOL scripts and optional CSV summaries for stems and inter-stem parameters.

Parameters:

Name	Type	Description	Default
args	Namespace	Parsed command-line arguments.	required
structure2d	Structure2D	Secondary structure produced by the annotator.	required
mapping	Mapping2D3D	Mapping between 2D indices and 3D residues.	required
input_filename	str	Name of the input PDB/mmCIF file.	required

Source code in src/rnapolis/annotator.py

def handle_output_arguments(
    args: argparse.Namespace,
    structure2d: Structure2D,
    mapping: Mapping2D3D,
    input_filename: str,
):
    """
    Handle all output options based on parsed CLI arguments.

    Depending on the flags, this function writes BPSEQ, CSV, JSON, DOT,
    PyMOL scripts and optional CSV summaries for stems and inter-stem
    parameters.

    Args:
        args: Parsed command-line arguments.
        structure2d: Secondary structure produced by the annotator.
        mapping: Mapping between 2D indices and 3D residues.
        input_filename: Name of the input PDB/mmCIF file.
    """
    input_basename = os.path.basename(input_filename)
    if args.csv:
        write_csv(args.csv, structure2d)

    if args.json:
        write_json(args.json, structure2d)

    if args.bpseq:
        write_bpseq(args.bpseq, structure2d.bpseq)

    if args.extended:
        print(structure2d.extended_dot_bracket)
    else:
        print(structure2d.dot_bracket)

    if args.dot:
        print(structure2d.bpseq.graphviz)

    if args.pml:
        pml_script = generate_pymol_script(mapping, structure2d.stems)
        with open(args.pml, "w") as f:
            f.write(pml_script)

    if args.inter_stem_csv:
        if structure2d.inter_stem_parameters:
            # Convert list of dataclasses to list of dicts
            params_list = [
                {
                    "stem1_idx": p.stem1_idx,
                    "stem2_idx": p.stem2_idx,
                    "type": p.type,
                    "torsion": p.torsion,
                    "min_endpoint_distance": p.min_endpoint_distance,
                    "torsion_angle_pdf": p.torsion_angle_pdf,
                    "min_endpoint_distance_pdf": p.min_endpoint_distance_pdf,
                    "coaxial_probability": p.coaxial_probability,
                }
                for p in structure2d.interStemParameters
            ]
            df = pd.DataFrame(params_list)
            df["input_basename"] = input_basename
            # Reorder columns to put input_basename first
            cols = ["input_basename"] + [
                col for col in df.columns if col != "input_basename"
            ]
            df = df[cols]
            df.to_csv(args.inter_stem_csv, index=False)
        else:
            logging.warning(
                f"No inter-stem parameters calculated for {input_basename}, "
                f"CSV file '{args.inter_stem_csv}' will be empty or not created."
            )

    if args.stems_csv:
        if structure2d.stems:
            stems_data = []
            for i, stem in enumerate(structure2d.stems):
                try:
                    res5p_first = mapping.bpseq_index_to_residue_map.get(
                        stem.strand5p.first
                    )
                    res5p_last = mapping.bpseq_index_to_residue_map.get(
                        stem.strand5p.last
                    )
                    res3p_first = mapping.bpseq_index_to_residue_map.get(
                        stem.strand3p.first
                    )
                    res3p_last = mapping.bpseq_index_to_residue_map.get(
                        stem.strand3p.last
                    )

                    stems_data.append(
                        {
                            "stem_idx": i,
                            "strand5p_first_nt_id": res5p_first.full_name
                            if res5p_first
                            else None,
                            "strand5p_last_nt_id": res5p_last.full_name
                            if res5p_last
                            else None,
                            "strand3p_first_nt_id": res3p_first.full_name
                            if res3p_first
                            else None,
                            "strand3p_last_nt_id": res3p_last.full_name
                            if res3p_last
                            else None,
                            "strand5p_sequence": stem.strand5p.sequence,
                            "strand3p_sequence": stem.strand3p.sequence,
                        }
                    )
                except KeyError as e:
                    logging.warning(
                        f"Could not find residue for stem {i} (index {e}), "
                        "skipping stem details."
                    )
                    continue

            if stems_data:
                df_stems = pd.DataFrame(stems_data)
                df_stems["input_basename"] = input_basename
                # Reorder columns
                stem_cols = ["input_basename", "stem_idx"] + [
                    col
                    for col in df_stems.columns
                    if col not in ["input_basename", "stem_idx"]
                ]
                df_stems = df_stems[stem_cols]
                df_stems.to_csv(args.stems_csv, index=False)
            else:
                logging.warning(
                    f"No valid stem data generated for {input_basename}, "
                    f"CSV file '{args.stems_csv}' will be empty or not created."
                )
        else:
            logging.warning(
                f"No stems found for {input_basename}, CSV file "
                f"'{args.stems_csv}' will be empty or not created."
            )

main()

Command-line entry point for the annotator tool.

It reads a 3D structure, extracts all base-level interactions, builds a secondary structure and writes outputs requested by the user.

Source code in src/rnapolis/annotator.py

def main():
    """
    Command-line entry point for the annotator tool.

    It reads a 3D structure, extracts all base-level interactions,
    builds a secondary structure and writes outputs requested by the user.
    """
    parser = argparse.ArgumentParser()
    parser.add_argument("input", help="Path to PDB or mmCIF file")
    parser.add_argument(
        "-f",
        "--find-gaps",
        action="store_true",
        help=(
            "(optional) if set, the program will detect gaps and break the PDB "
            "chain into two or more strands; the gap is defined as O3'-P "
            f"distance greater then {1.5 * AVERAGE_OXYGEN_PHOSPHORUS_DISTANCE_COVALENT}"
        ),
    )
    add_common_output_arguments(parser)
    args = parser.parse_args()

    file = handle_input_file(args.input)
    structure3d = read_3d_structure(file, None)
    base_interactions = extract_base_interactions(structure3d)
    structure2d, mapping = structure3d.extract_secondary_structure(
        base_interactions, args.find_gaps, args.decompose_pseudoknot_free
    )

    handle_output_arguments(args, structure2d, mapping, args.input)

merge_and_clean_bph_br(pairs)

Merge and simplify base–phosphate and base–ribose classifications.

Multiple classifications for the same residue pair are merged into a set and then cleaned according to simple rules (for example 3+5 → 4).

Parameters:

Name	Type	Description	Default
pairs	List[Tuple[Residue3D, Residue3D, int]]	List of residue pairs with their classification codes.	required

Returns:

Type	Description
Dict[Tuple[Residue3D, Residue3D], OrderedSet[int]]	Mapping from residue pairs to a cleaned set of classification codes.

Source code in src/rnapolis/annotator.py

def merge_and_clean_bph_br(
    pairs: List[Tuple[Residue3D, Residue3D, int]],
) -> Dict[Tuple[Residue3D, Residue3D], OrderedSet[int]]:
    """
    Merge and simplify base–phosphate and base–ribose classifications.

    Multiple classifications for the same residue pair are merged into a set
    and then cleaned according to simple rules (for example 3+5 → 4).

    Args:
        pairs: List of residue pairs with their classification codes.

    Returns:
        Mapping from residue pairs to a cleaned set of classification codes.
    """
    bph_br_map: Dict[Tuple[Residue3D, Residue3D], OrderedSet[int]] = defaultdict(
        OrderedSet
    )
    for residue_i, residue_j, classification in pairs:
        bph_br_map[(residue_i, residue_j)].add(classification)
    for bphs_brs in bph_br_map.values():
        # 3BPh and 5BPh simultanously means that it is actually 4BPh
        if 3 in bphs_brs and 5 in bphs_brs:
            bphs_brs.remove(3)
            bphs_brs.remove(5)
            bphs_brs.add(4)
        # 7BPh and 9BPh simultanously means that it is actually 8BPh
        if 7 in bphs_brs and 9 in bphs_brs:
            bphs_brs.remove(7)
            bphs_brs.remove(9)
            bphs_brs.add(8)
    for key, bphs_brs in bph_br_map.items():
        if len(bphs_brs) > 1:
            bph_br_map[key] = OrderedSet([bphs_brs[0]])
    return bph_br_map

write_bpseq(path, bpseq)

Write a BpSeq object to a plain text file.

The BpSeq object is converted to a string using its __str__ method.

Source code in src/rnapolis/annotator.py

def write_bpseq(path: str, bpseq: BpSeq):
    """
    Write a BpSeq object to a plain text file.

    The BpSeq object is converted to a string using its ``__str__`` method.
    """
    with open(path, "w") as f:
        f.write(str(bpseq))

write_csv(path, structure2d)

Write base interactions from Structure2D to a CSV file.

The CSV contains base pairs, stacking, base–phosphate and base–ribose interactions in a flat table.

Source code in src/rnapolis/annotator.py

def write_csv(path: str, structure2d: Structure2D):
    """
    Write base interactions from Structure2D to a CSV file.

    The CSV contains base pairs, stacking, base–phosphate and base–ribose
    interactions in a flat table.
    """
    with open(path, "w") as f:
        writer = csv.writer(f)
        writer.writerow(["nt1", "nt2", "type", "classification-1", "classification-2"])
        for base_pair in structure2d.base_pairs:
            writer.writerow(
                [
                    base_pair.nt1.full_name,
                    base_pair.nt2.full_name,
                    "base pair",
                    base_pair.lw.value,
                    (
                        base_pair.saenger.value or ""
                        if base_pair.saenger is not None
                        else ""
                    ),
                ]
            )
        for stacking in structure2d.stackings:
            writer.writerow(
                [
                    stacking.nt1.full_name,
                    stacking.nt2.full_name,
                    "stacking",
                    stacking.topology.value if stacking.topology is not None else "",
                    "",
                ]
            )
        for base_phosphate in structure2d.base_phosphate_interactions:
            writer.writerow(
                [
                    base_phosphate.nt1.full_name,
                    base_phosphate.nt2.full_name,
                    "base-phosphate interaction",
                    base_phosphate.bph.value if base_phosphate.bph is not None else "",
                    "",
                ]
            )
        for base_ribose in structure2d.base_ribose_interactions:
            writer.writerow(
                [
                    base_ribose.nt1.full_name,
                    base_ribose.nt2.full_name,
                    "base-ribose interaction",
                    base_ribose.br.value if base_ribose.br is not None else "",
                    "",
                ]
            )
        for other in structure2d.other_interactions:
            writer.writerow(
                [
                    other.nt1.full_name,
                    other.nt2.full_name,
                    "other interaction",
                    "",
                    "",
                ]
            )

write_json(path, structure2d)

Write a Structure2D object to a JSON file using orjson.

The function copies the object, replaces residues with simpler Residue objects and serialises everything with options that handle numpy types and non-string dictionary keys.

Source code in src/rnapolis/annotator.py

def write_json(path: str, structure2d: Structure2D):
    """
    Write a Structure2D object to a JSON file using orjson.

    The function copies the object, replaces residues with simpler Residue
    objects and serialises everything with options that handle numpy types
    and non-string dictionary keys.
    """
    processed = copy.deepcopy(structure2d)
    processed.bpseq_index = {
        k: Residue(v.label, v.auth) for k, v in structure2d.bpseq_index.items()
    }

    with open(path, "wb") as f:
        # Add OPT_SERIALIZE_NUMPY to handle numpy types like float64
        # Add OPT_NON_STR_KEYS to preserve integer keys in dictionaries
        f.write(
            orjson.dumps(
                processed, option=orjson.OPT_SERIALIZE_NUMPY | orjson.OPT_NON_STR_KEYS
            )
        )