Biopython: MMTF writer

[harryjubb]

Hi,

Are there any plans to add an MMTF writer (MMTFIO?) to BioPython? Is this something that could be contributed?

My use case is for converting PDB files to MMTF. I had a go at taking a BioPython Structure and looping through the MCRA model into the MMTFEncoder, which for the most part worked apart from for bonding information (viewing in the NGL web app shows the atoms but no bonds).

I'm not sure if CONECT or otherwise perceived bond info is parsed/stored for PDB files in BioPython, from what I can see, this would need to be added also?

[peterjc]

As you've noticed, the mmft Python library has the MMTFEncoder which would support writing MMTF files.

Probably @abradle would be the best person to ask about how realistic this is, since he wrote the code for loading MMTF files into Biopython. It may be that too much data gets dropped on conversion to the PDB centric object model?

[abradle]

What you did makes sense @harryjubb and yes I think by default bonds are not parsed. I can't fully remember if other options enable bond parsing - but I seem to remember what you say is correct - and it is not possible and would need implementing.

If you want to convert PDB -> MMTF the simplest way would be via BioJava(https://biojava.org/) since this implementation works well for reading PDB files and writing MMTF and would create bonds.

[peterjc]

So right now (without also improving the MMTF parsing code in Biopython to handle bonds), it should at least be possible to at least output MMTF records recording the atom positions?

[harryjubb]

Thanks both for your replies :slightly_smiling_face:

It should at least be possible to at least output MMTF records recording the atom positions?

That's right, I could get a parsable MMTF file with all the atom positions:

Using:

import mmtf
from Bio.PDB import PDBParser

from mmtf.api.mmtf_writer import MMTFEncoder

def write_mmtf(biopython_structure, filename):

    encoder = MMTFEncoder()
    encoder.init_structure(
        total_num_bonds=0,
        total_num_atoms=len(list(biopython_structure.get_atoms())),
        total_num_groups=len(list(biopython_structure.get_residues())),
        total_num_chains=len(list(biopython_structure.get_chains())),
        total_num_models=len(list(biopython_structure.get_models())),
        structure_id=biopython_structure.id
    )

    # TODO: XTAL INFO
    encoder.set_xtal_info(
        space_group='',
        unit_cell=[0,0,0,0,0,0]
    )

    # TODO: BIO ASSEMBLY INFO

    # TODO: HEADER INFO
    encoder.set_header_info(
        r_free=0.0,
        r_work=0.0,
        resolution=0.0,
        title='',
        deposition_date='',
        release_date='',
        experimental_methods=''
    )

    for e, model in enumerate(biopython_structure.get_models()):

        encoder.set_model_info(
            model_id=e,
            chain_count=len(list(model.get_chains()))
        )

        for chain in model.get_chains():
            encoder.set_chain_info(
                chain_id=chain.get_id(),
                chain_name=chain.get_id(),
                num_groups=0
            )

            for residue in chain.get_residues():
                encoder.set_group_info(
                    group_name=residue.resname,
                    group_number=residue.id[1],
                    insertion_code=residue.id[2],
                    group_type='',
                    atom_count=len(list(residue.get_atoms())),
                    bond_count=0,
                    single_letter_code='',
                    sequence_index=0,
                    secondary_structure_type=0
                )   

                for atom in residue.get_atoms():
                    encoder.set_atom_info(
                        atom_name=atom.name,
                        serial_number=atom.serial_number,
                        alternative_location_id=atom.altloc,
                        x=atom.coord[0],
                        y=atom.coord[1],
                        z=atom.coord[2],
                        occupancy=atom.occupancy,
                        temperature_factor=atom.bfactor,
                        element=atom.element,
                        charge=0
                    )

    encoder.finalize_structure()
    encoder.write_file(filename)

pdb_parser = PDBParser()
structure = pdb_parser.get_structure('1xkk', './1xkk.pdb')

write_mmtf(structure, '1xkk.mmtf')

I imagine that while it doesn't cover all the edge cases, the above plus bond information plus the other missing metadata for the encoder would be a good start.

Presumably for bonds, the BioPython PDB parser would need to be extended to parse CONECT records, and fall back on a distance-based approach?

Thanks also @abradle, I'll take a look at the BioJava solution too :slightly_smiling_face:

[peterjc]

As a I recall (and it has been many years since I worked on PDB files), the covalent bonds are implicit via distance (and for proteins the residue type would also give you a prior for how the expected side chain would be connected) only a few special cases like cystine bridges would be annotated.

I take it MMTF expects the bonds to be recorded explicitly?

[abradle]

Yes it does. And distance based bond perception can be troubling. Bio java
uses the wwpdb bond information for assigning residue bonds. Ligands unless
in that dictionary would be more troublesome.

Harry - this discussion may also be worth keeping an eye on. Specifically
if you see Steve roughley at all he might be able to help you too.

https://github.com/rcsb/mmtf-java/issues/22#issuecomment-293087298

On Thu, 14 Mar 2019 at 15:50, Peter Cock [e-mail] wrote:

As a I recall (and it has been many years since I worked on PDB files),
the covalent bonds are implicit via distance (and for proteins the residue
type would also give you a prior for how the expected side chain would be
connected) only a few special cases like cystine bridges would be annotated.

I take it MMTF expects the bonds to be recorded explicitly?

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
https://github.com/biopython/biopython/issues/1956#issuecomment-472929906,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AFuRIRBDYSiUnFvj8MS25EKsly5bPVmzks5vWm-qgaJpZM4brd34
.

[harryjubb]

I can see where that would be problematic. Unfortunately my usage is outside a wwPDB context, there's unlikely to be existing mmCIF or MMTF files available, nor any ChemComp dictionary. There's probably no nice solutions! BioJava might be the way at the moment.

I'll keep an eye here and on the mmtf-java issue, thanks for linking to that. I'll continue to have a play though probably at "spare time project" rate :slightly_smiling_face:

[jgreener64]

This would be a cool addition, and would complete the last of the 6 {PDB, mmCIF, MMTF}{read, write} combinations.

[jgreener64]

There is a proposal for this in #2185.

[peterjc]

With #2185 merged I forgot to close this.