Residues added to N-terminus of a chain do not appear in Model.polymers

[Tristan Croll]

If I use the (new today, still a bit messy) code below (add_amino_acid_residue) to do a quick-and-dirty addition of a new residue to the N-terminus of an existing chain, then *most* things work - bonds and cartoons draw correctly, and simulations run - but the new residue doesn't appear in Model.polymers(missing_structure_treatment=Model.PMS_NEVER_CONNECTS).

Code below. The call is something like add_amino_acid_residue(model, 'PHE', next_res=selected_residues(session)[0]).

While playing around with attempts to debug this, I also noticed that Model.reorder_residues() segfaults. Will raise a separate ticket.

def set_new_atom_style(session, atoms):
    from chimerax.atomic import selected_atoms, Atom
    from chimerax.core.commands import run
    atoms.draw_modes = Atom.STICK_STYLE
    current_sel = selected_atoms(session)
    session.selection.clear()
    atoms.selected = True
    run(session, "color sel bychain; color sel byhetero")
    session.selection.clear()
    current_sel.selected = True

def fix_amino_acid_protonation_state(residue):
    if residue.name not in ('GLU', 'ASP'):
        raise TypeError('This method is only applicable to GLU and ASP residues!')
    if residue.name == 'GLU':
        a = residue.find_atom('HE2')
        if a is not None:
            a.delete()
    elif residue.name == 'ASP':
        a = residue.find_atom('HD2')
        if a is not None:
            a.delete()

def new_residue_from_template(model, template, chain_id, center, residue_number=None, insert_code = ' ', b_factor = 50):
    '''
    Create a new residue based on a template, and add it to the model.
    '''
    if residue_number is None:
        if chain_id in m.residues.chain_ids:
            residue_number = suggest_new_residue_number_for_ligand(model, chain_id)
        else:
            residue_number = 0
    import numpy
    from chimerax.atomic import Atom
    t_coords = numpy.array([a.coord for a in template.atoms])
    t_center = t_coords.mean(axis=0)
    t_coords += numpy.array(center) - t_center
    tatom_to_atom = {}
    r = model.new_residue(template.name, chain_id, residue_number, insert=insert_code)
    for i, ta in enumerate(template.atoms):
        a = tatom_to_atom[ta] = model.new_atom(ta.name, ta.element)
        a.coord = t_coords[i]
        a.bfactor = b_factor
        r.add_atom(a)
        for tn in ta.neighbors:
            n = tatom_to_atom.get(tn, None)
            if n is not None:
                model.new_bond(a, n)
    set_new_atom_style(model.session, r.atoms)
    return r

def add_amino_acid_residue(model, resname, prev_res=None, next_res=None,
        chain_id=None, number=None, center=None, insertion_code=' ', b_factor=50,
        occupancy=1, phi=-135, psi=135):
    session = model.session
    if (not chain_id or not number or not center) and (not prev_res and not next_res):
        raise TypeError('If no anchor residues are specified, chain ID, '
            'number and center must be provided!')
    if prev_res and next_res:
        raise TypeError('Cannot specify both previous and next residues!')
    if prev_res:
        for n in prev_res.find_atom('C').neighbors:
            if n.residue != prev_res:
                raise TypeError('This residue already has another bonded to its '
                    'C terminus!')
        chain_id = prev_res.chain_id
        oxt = prev_res.find_atom('OXT')
        if oxt is not None:
            oxt.delete()
    elif next_res:
        for n in next_res.find_atom('N').neighbors:
            if n.residue != next_res:
                raise TypeError('This residue already has another bonded to its '
                    'N terminus!')
        chain_id = next_res.chain_id
        for hname in ('H2', 'H3'):
            h = next_res.find_atom(hname)
            if h is not None:
                h.delete()
            h = next_res.find_atom('H1')
            if h is not None:
                h.name='H'
            if next_res.name == 'PRO':
                h = next_res.find_atom('H')
                if h:
                    h.delete()
    if number is None:
        if prev_res:
            number = prev_res.number + 1
        elif next_res:
            number = next_res.number - 1

    from chimerax.atomic import mmcif
    tmpl = mmcif.find_template_residue(session, resname)
    from .place_ligand import new_residue_from_template
    import numpy
    # delete extraneous atoms
    r = new_residue_from_template(model, tmpl, chain_id, [0,0,0], number,
            insert_code=insertion_code, b_factor=b_factor)
    r.atoms[numpy.in1d(r.atoms.names, ['OXT', 'HXT', 'H2'])].delete()

    # Translate and rotate residue to (roughly) match the desired position
    if not next_res and not prev_res:
        r.atoms.coords += numpy.array(center) - r.atoms.coords.mean(axis=0)
    else:
        from chimerax.core.geometry import align_points
        if prev_res:
            model.new_bond(r.find_atom('N'), prev_res.find_atom('C'))
            n_pos = _find_next_N_position(prev_res)
            ca_pos = _find_next_CA_position(n_pos, prev_res)
            c_pos = _find_next_C_position(ca_pos, n_pos, prev_res, phi)
            target_coords = numpy.array([n_pos, ca_pos, c_pos])
            align_coords = numpy.array([r.find_atom(a).coord for a in ['N', 'CA', 'C']])
        elif next_res:
            model.new_bond(r.find_atom('C'), next_res.find_atom('N'))
            c_pos = _find_prev_C_position(next_res, psi)
            ca_pos = _find_prev_CA_position(c_pos, next_res)
            o_pos = _find_prev_O_position(c_pos, next_res)
            target_coords = numpy.array([c_pos, ca_pos, o_pos])
            align_coords = numpy.array([r.find_atom(a).coord for a in ['C', 'CA', 'O']])

        tf = align_points(align_coords, target_coords)[0]
        r.atoms.coords = tf*r.atoms.coords
    if r.name in ('GLU', 'ASP'):
        fix_amino_acid_protonation_state(r)
    if r.name == 'PRO':
        r.atoms[r.atoms.names=='H'].delete()
    r.atoms.bfactors = b_factor
    r.atoms.occupancies = occupancy

    set_new_atom_style(session, r.atoms)
    return r



def _find_next_N_position(prev_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.34
    angle = 120
    dihedral = 180
    a1 = prev_res.find_atom('C')
    a2 = prev_res.find_atom('CA')
    a3 = prev_res.find_atom('O')
    return find_pt(*[a.coord for a in [a1, a2, a3]], bond_length, angle, dihedral)

def _find_next_CA_position(n_pos, prev_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.48
    angle = 124
    omega = 180
    c = prev_res.find_atom('C')
    ca = prev_res.find_atom('CA')
    return find_pt(n_pos, *[a.coord for a in [c, ca]], bond_length, angle, omega)

def _find_next_C_position(ca_pos, n_pos, prev_res, phi):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.53
    angle = 120
    c = prev_res.find_atom('C')
    return find_pt(ca_pos, n_pos, c.coord, bond_length, angle, phi)


def _find_prev_C_position(next_res, psi):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.34
    angle = 120
    a1 = next_res.find_atom('N')
    a2 = next_res.find_atom('CA')
    a3 = next_res.find_atom('C')
    return find_pt(*[a.coord for a in [a1, a2, a3]], bond_length, angle, psi)

def _find_prev_CA_position(c_pos, next_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.53
    angle = 120
    omega = 180
    n = next_res.find_atom('N')
    ca = next_res.find_atom('CA')
    return find_pt(c_pos, *[a.coord for a in [n, ca]], bond_length, angle, omega)

def _find_prev_O_position(c_pos, next_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.22
    angle = 120
    dihedral = 0
    n = next_res.find_atom('N')
    ca = next_res.find_atom('CA')
    return find_pt(c_pos, *[a.coord for a in (n, ca)], bond_length, angle, dihedral)

[pett]

I am having trouble reproducing this. Certainly, if I take a structure with two polymer chains and form a polymeric bond between them, I get one polymer as the result. I could not use the code you provided because it is not complete. Namely:

from .place_ligand import new_residue_from_template

Perhaps you could provide a self-contained example that reproduces the error? (Preferably as an attachment rather than inline code!)

--Eric

[Tristan Croll]

Ah - sorry about the unreproducability! I tried to make it self-contained by bundling in new_residue_from_template(), but forgot to remove the import statement. Anyway, it looks like this was fixed by your addition of the precedes keyword to new_residue() (an aside: for the sake of code simplicity it might be nice to also have a follows keyword...?). Anyway, this code is all in the ISOLDE builds I released last night:

from chimerax.core.commands import open as cxopen

m = cxopen.open(session, '3io0')[0]

from chimerax.isolde.atomic.building.build_utils import add_amino_acid_residue

r = add_amino_acid_residue(m, 'ALA', next_res=m.residues[0])

pols = m.polymers(missing_structure_treatment=m.PMS_NEVER_CONNECTS)

pols
Out[6]: [(<chimerax.atomic.molarray.Residues at 0x7fa01872ef90>, 1)]

r in pols[0][0]
Out[7]: True

[pett]

I typically don't like to add mutually exclusive keywords to a command. The only time you need 'follows' is at the end of all residues, which is the same as omitting the 'precedes' keyword. I guess I do understand the convenience in some usage scenarios.