OneScience-Group/flex_ddG_tutorial

1#!/usr/bin/env python32"""3Generate a copy of ddG-backrub.xml that additionally reports per-chain energies.4 5Why6---7The `unbound_wt` / `unbound_mut` totals in ddG.db3 are whole-pose scores. InterfaceDdGMover8builds the unbound state by rigid-body translating the moving chain(s) 1000 A apart, without9repacking or re-minimizing, and every chain stays in the pose10(protocols/features/InterfaceDdGMover.cc::unbind). So11 12    unbound_X_total = sum over chains of intra(chain)        for X in {wt, mut}13 14and `unbound_mut - unbound_wt` is the summed intramolecular ddG of *every* chain, including15the ones that were never mutated. Those only differ between the wild type and mutant branches16because the whole pose is re-minimized independently in each, so they contribute artifact.17 18What this adds19--------------20A `TotalEnergyMetric` per chain, reported into ddG.db3 by the same ReportToDB mover that21InterfaceDdGMover already applies to all four states. This is *reporting only* -- it adds no22sampling, consumes no random numbers, and leaves ddG_bind bit-identical (verified against the23unmodified protocol with -constant_seed).24 25Reading the result26------------------27On the **unbound** states the chains are 1000 A apart, so there are no cross-chain pair28energies and the reported value is exactly that chain's intramolecular energy. That is the29number to use:30 31    per-chain ddG(chain X) = intra_mut(X) - intra_wt(X)32 33On the **bound** states the value is intra(chain) + roughly half the interface energy, because34Rosetta's residue_total_energies splits each two-body term between its two residues. Those rows35are useful as a cross-check (they sum to the pose total) but should not be read as per-chain36stability numbers.37 38Caveat this does NOT fix39------------------------40The unbound state is still never relaxed. This gives you the intramolecular strain difference41*in the bound backbone conformation*, not a folding ddG of the free monomer. For that you want42a dedicated monomer protocol (e.g. cartesian_ddg) run on the isolated chain.43"""44 45import os46 47 48def chains_in_pdb(pdb_path):49    """Chain IDs in the order they first appear in the PDB."""50    chains = []51    with open(pdb_path) as f:52        for line in f:53            if line.startswith(('ATOM', 'HETATM')):54                chain = line[21]55                if chain not in chains:56                    chains.append(chain)57    return chains58 59 60def write_per_chain_protocol(base_xml_path, chains, out_path, scorefxn='fa_talaris2014'):61    """Write a copy of base_xml_path with a per-chain TotalEnergyMetric for each chain.62 63    Each chain gets its own database table (chain_<X>_energy) with a single `total_energy`64    column. Separate tables rather than one table with prefixed columns is deliberate: as of65    Rosetta 2022.45, SimpleMetricFeatures builds schema column names as `custom_type + name`66    but builds the INSERT as `custom_type + "_" + name`, so any non-empty custom_type produces67    "table simple_metrics has no column named ..." at report time.68    """69    if not chains:70        raise ValueError('No chains given')71 72    # Rosetta lowercases table names, so chains differing only in case would collide.73    lowered = [c.lower() for c in chains]74    if len(set(lowered)) != len(lowered):75        raise ValueError('Chain IDs differ only by case, which collides in table names: %s' % chains)76 77    xml = open(base_xml_path).read()78 79    def substitute(text, anchor, addition):80        if anchor not in text:81            raise ValueError('Could not find anchor in %s:\n  %s' % (base_xml_path, anchor))82        return text.replace(anchor, anchor + addition, 1)83 84    selectors = ''.join(85        '\n    <Chain name="chain_%s" chains="%s"/>' % (c, c) for c in chains)86    xml = substitute(87        xml, '    <StoredResidueSubset name="restore_neighbor_shell" subset_name="neighbor_shell"/>',88        selectors)89 90    metrics = '\n'.join(91        '    <TotalEnergyMetric name="chain_%s_energy" residue_selector="chain_%s" scorefxn="%s"/>'92        % (c, c, scorefxn) for c in chains)93    anchor = '  <FILTERS>'94    if anchor not in xml:95        raise ValueError('Could not find FILTERS block in %s' % base_xml_path)96    xml = xml.replace(anchor, '  <SIMPLE_METRICS>\n%s\n  </SIMPLE_METRICS>\n\n%s' % (metrics, anchor), 1)97 98    reporters = ''.join(99        '\n      <SimpleMetricFeatures metrics="chain_%s_energy" table_name="chain_%s_energy"/>'100        % (c, c) for c in chains)101    xml = substitute(102        xml, '      <StructureScoresFeatures scorefxn="%s"/>' % scorefxn, reporters)103 104    with open(out_path, 'w') as f:105        f.write(xml)106    return out_path107 108 109if __name__ == '__main__':110    import sys111    if len(sys.argv) != 4:112        sys.exit('usage: per_chain_protocol.py <base.xml> <input.pdb> <out.xml>')113    base_xml, pdb, out = sys.argv[1:4]114    found = chains_in_pdb(pdb)115    write_per_chain_protocol(base_xml, found, out)116    print('Wrote %s with per-chain metrics for chains: %s' % (out, ', '.join(found)))117