[API] revisiting in-memory representation of `pdb`/ `mmcif`

[fkiraly]

We need a 1:1 in-memory representation of pdb / mmcif in python, for a clean design.

@satvshr in #148 (comment) claimed the biopython Structure object is not 1:1 but lossy.

We need to investigate and find an alternative 1:1 representation if this proves to be the case.

[fkiraly]

To summarize our discussion earlier today: @satvshr thinks the Structure object is lossy, because the logic in current struct_to_aaseq is lossy - pp.get_sequence retrieves not the full AA sequence but only a partial sequence, according to @KubiczekD, @JaBirke.

This does not logically imply that Strucure is lossy, i.e., that the full AA sequence cannot be retrieved in a different way somehow, but it looks like a strong indication that it is.

We should verify, e.g., by internal inspection of the Structure object whether it is indeed lossy - and if yes, what it is missing.

[KubiczekD]

I would like to summarize the current state from our point of view in the hope of making it clearer for everyone.
As we have already discussed several times, the amino acid sequence derived from the ATOM section of the PDB file can, and in many cases will, be lossy. We talked about this, for example, in Issue #8.
However, the PDB file will always contain the SEQRES section, which includes the full sequence.
→ We “just” need a way to extract the sequence provided in the SEQRES section. I am a little confused that there is no possibility.

From a user’s perspective, we would expect the following workflow in AptaNET:
Training
A dataset (e.g., a CSV file) contains known aptamer–protein pairs and the information on whether they bind or not. As a next step, we are also implementing the possibility to train and predict affinity (Kd value).
Prediction
Once the training is complete, we want to predict the binding of aptamers to a new target. At this stage, we provide:

• The PDB file representing our target, or alternatively, the option to input the sequence directly if the PDB file is not available. However using PDB files makes sense, as we plan to leverage their 3D structural information for downstream tasks such as docking simulations.
• The corresponding FASTA file

The combination of sequences from our FASTA file and the PDB file is what we currently designate as X_new.

[fkiraly]

I see, @KubiczekD, I assume you are referring to this post? #8 (comment)
"the sequence @satvshr prepared with biopython can be used but it is incomplete"

I did not catch that - or I might have read it that it was just a reference to some processing @satvshr did.

I think it is still not 100% clear whether the biopython Structure object represents pdf completely, but if it does not, it would be troubling and cannot be used as an in-memory correlate of pdb (or mmcif). Having had a cursory look, it is probably lossy in its entirety (as a python object), not just the ATOM sequence from it.

I am a little confused that there is no possibility.

I think this is a false statement and due to incomplete literature research. The natural place to look would have been what biopython uses internally to load the pdb files.

Looking into biopython, it seems to use the package mmtf-python underneath which loads pdb files. It was updated last in 2022 but it seems to be maintained by rcsb who also maintain the protein databank.
https://github.com/rcsb/mmtf-python
This seems to parse the file completely, so I think we should rely on this for now?

We should schedule this as a priority item and close out the data container question with clear directives.

[fkiraly]

@siddharth7113 mentioned this loader from openmm:
https://docs.openmm.org/7.0.0/api-python/generated/simtk.openmm.app.pdbfile.PDBFile.html
upon cursory look, I think this is based on the ATOM record and hence also lossy
(to be checked)

[fkiraly]

@rpgv mentioned that pdb files often do not have the full sequence, missing the SEQRES record

[fkiraly]

This would mean we have to take into account both pdb that have SEQRES and pdb that do not.

[NennoMP]

Maybe these discussions are related to our problem with BioPython, if you want to look at them:

• https://stackoverflow.com/questions/45466408/biopython-resseq-doesnt-match-pdb-file
• https://bioinformatics.stackexchange.com/questions/14101/extract-residue-sequence-from-pdb-file-in-biopython-but-open-to-recommendation

[fkiraly]

from discussion:

• we need to replace the loaders in the notebooks by mmtf-python
• in future designs, we should use the mmtf-python objects as primary in-memory representation

[fkiraly]

@satvshr pointed out that SeqIO is doing what mmtf-python isdoing.

@satvshr, can you please check whether the result of SeqIO is lossy or not - and whether it is an alternative to mmtf-python, or uses it under the hood.

[NennoMP]

Documentation for SeqIO: https://biopython.org/docs/1.75/api/Bio.SeqIO.PdbIO.html

[rpgv]

Discussion regarding Sequence completeness: In case SeqRes information is not present in all PDB files, it could be possible to implement UniProt sequence information.
i.e: From PDB id obtain UniProt id, from UniProt id retrieve complete sequence;

[JaBirke]

@rpgv I had similar thoughts. If the PDB format often lacks the SEQRES information, but we want to keep it as a standard file (needed for MAWS), this would be a good option. I found UniProdMapper, a Python package that might be suited to this purpose.
https://pypi.org/project/uniprot-id-mapper/#mapping-ids

[satvshr]

The design I think we should follow in terms of in-memory representation of file formats:

1. Use SeqIO as a loader for all possible file formats.
2. Use SeqRecord as the in-memory format (SeqIO loads a SeqRecord object too).

All notebooks, loaders, and anything else using Structure or String should be updated to use SeqRecord. A few questions I have:

1. Is SeqRecord a good representation as it does not capture 3d information as compared to Structure?
2. Should we change existing algorithms to handle Seq objects (comes along with SeqRecord) instead of String or not?

Regarding the pdf loader, it should be edited to handle a pdb file in this way (will change #148 based on this):

1. Use SEQRES records if available.
2. Use ATOM records if not.
3. Make the loader have an argument asking if the user wants to use the UniProt id to get the full amino acid sequence given the pdb id, if yes, make 2 API calls:
   1. To get the uniprot id using the user provided pdb id,
   2. To get the fasta file using the uniprot id and load it in memory using the fasta loader.

At the time of writing this I could not access the mmtf website (http://mmtf.rcsb.org/) or use the fetch function mentioned in their README file as the function was throwing a URLError (Python couldn’t resolve the hostname).