Implementing gemmi-based mmcif reader (with easy extension to PDB/PDBx and mmJSON)
#4712
Conversation
|
Hello @marinegor! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2024-10-25 11:17:29 UTC |
Linter Bot Results:Hi @marinegor! Thanks for making this PR. We linted your code and found the following: Some issues were found with the formatting of your code.
Please have a look at the Please note: The |
| pass | ||
|
|
||
|
|
||
| class MMCIFWriter(base.WriterBase): |
There was a problem hiding this comment.
I wouldn't include this at this stage, Writer is optional
| from .base import TopologyReaderBase | ||
|
|
||
|
|
||
| def _into_idx(arr: list[int]) -> list[int]: |
There was a problem hiding this comment.
Document what this does, ideally with an example
| record_types, # res.het_flag | ||
| tempfactors, # at.b_iso | ||
| residx, # _into_idx(res.seqid.num) | ||
| ) = map( # this is probably not pretty, but it's efficient -- one loop over the mmcif |
There was a problem hiding this comment.
Are all the fields here guaranteed in a valid pdbx? One benefit to working column by column is that you can do optional columns
There was a problem hiding this comment.
Do you have an example of a PDBx in mind, or like a test set for them? I've never actually worked with the format, since in RCSB afaik we have only pdb or mmcif
There was a problem hiding this comment.
PDBx is mmcif. The download links here will give you an example file: https://www.rcsb.org/structure/4ake we use 4ake elsewhere in the testsuite. In my experience, sometimes the PDB / mmcif versions of the same entry aren't completely identical, so I wouldn't worry about trying to align the PDB & PDBx tests.
| np.array, | ||
| list( | ||
| zip( | ||
| *[ |
There was a problem hiding this comment.
I'm struggling to follow the logic here, a comment breaking down what this double nested loop iteration into a zip is doing would be nice
package/pyproject.toml
Outdated
| @@ -78,6 +78,7 @@ extra_formats = [ | |||
| "pytng>=0.2.3", | |||
| "gsd>3.0.0", | |||
| "rdkit>=2020.03.1", | |||
| "gemmi", # for mmcif format | |||
There was a problem hiding this comment.
This will probably be optional, so other imports will have to respect that too
|
and my bad for the messed up commit history, I certainly have to not try switching git tools in the middle of a PR🥲 |
|
Thanks to @chmnk's contribution expressed as a test file, I've added multimodel structure warnings, and now I think codecov must be happy as well. |
|
all checks done, tests added -- pinging @yuxuanzhuang and @IAlibay (and also @richardjgowers could you please unassign yourself if you don't have the bandwidth?) |
|
I ran into a strange bug and don’t have time to investigate it right now. Could you take a look and see if you can fix it? import MDAnalysis as mda
from MDAnalysisTests.datafiles import MMCIF as MMCIF_FOLDER
u = mda.Universe(f"{MMCIF_FOLDER}/1YJP.cif")
print(u.atoms.resids)
> [ 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4
4 4 4 4 4 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 7 7
7 7 7 7 7 7 7 7 7 7 7 8 9 10 11 12 13 14]
# weird weird
print(u.select_atoms('resid 1'))
> <AtomGroup []>
# this should print both resid 1 and 2.
print(u.select_atoms('resid 1-2'))
<AtomGroup [<Atom 1: N of type N of resname GLY, resid 1 and segid A and altLoc A>, <Atom 2: CA of type CA of resname GLY, resid 1 and segid A and altLoc A>, <Atom 3: C of type C of resname GLY, resid 1 and segid A and altLoc A>, <Atom 4: O of type O of resname GLY, resid 1 and segid A and altLoc A>]> |
|
The issue can be traced back to how selections handle insertion codes (icodes) In mmCIF, icodes are read and stored as a list of A potential fix would be to strip all icodes during MMCIFParsing:
|
| ] | ||
| ) | ||
| if len(structure) > 1: | ||
| warnings.warn( # FIXME: add tests for this |
There was a problem hiding this comment.
Wait, you do have multimodel_warning.cif tested---I guess codecov is outdated?
|
I found a strange edge case that makes the parser do weird things even when using @yuxuanzhuang 's fix above: The This leads to cases where multiple different residues in the Gemmi struct have the same For example, from 1BD2: These two residues have a different As a result, for this PDB entry's cif file, when this line is called in residx = np.array(_into_idx(residx))len(np.unique(residx))
> 853However, the corresponding arrays in len(resids)
> 854This creates an indexing issue in the atom->residue mapping for all atoms after the atoms in the two residues shown above. The same issue does not exist in the equivalent legacy PDB file likely because the residues are uniquely identified as 54 and 54A (although I haven't read the PDB parser so this is a bit of a guess) To demonstrate the issue: u = mda.Universe(
"1bd2.pdb"
)
tmp = u.select_atoms("segid D and name CA")
tmp
> <AtomGroup with 190 atoms>But with cif: u_broken = mda.Universe(
"1bd2.cif"
)
tmp_brk = u_broken.select_atoms("segid D and name CA")
tmp_brk
> <AtomGroup with 191 atoms>
tmp_brk[-1]
# this is nonsense, doesn't exist in the topology
> <Atom 4518: CA of type CA of resname HOH, resid 223 and segid D and altLoc A>The actual atom 4518 from the topology: I found a few more such examples of PDB entries with cif files that I'm fairly certain are producing nonsense selections for the same reasion: Gemmi has some limited notes on auth vs label but otherwise I'm not quite sure if this is a PDB error or a parsing error that needs to be fixed- if it is a parsing error, either |
|
More concise examples: u = mda.Universe(
"1bd2.pdb"
)
tmp = u.select_atoms(f"resid 54 and segid D")
tmp.select_atoms("name CA")
> <AtomGroup with 1 atom>
tmp = u.select_atoms(f"resid 54A and segid D")
tmp.select_atoms("name CA")
> <AtomGroup with 1 atom>u_broken = mda.Universe(
"1bd2.cif"
)
tmp_brk = u_broken.select_atoms(f"resid 54 and segid D")
tmp_brk.select_atoms("name CA")
> <AtomGroup with 2 atoms>Before the broken res (54), the atoms we get from selections actually exist in the CIF file (although the repr of the atom seems to be broken, using atm = u_broken.select_atoms(f"resid 53 and segid D")[0]
print(atm)
print(atm.id)
> <Atom 3490: N of type N of resname LYS, resid 53 and segid D and altLoc A>
> 3476From CIF: After the broken res in the CIF, we get nonsense that doesn't exist in the cif: atm_brk = u_broken.select_atoms(f"resid 55 and segid D")[0]
print(atm_brk)
print(atm_brk.id)
> <Atom 3516: N of type N of resname ASN, resid 55 and segid D and altLoc A>
> 3502From CIF: |
Fixes #2367 and also extends #4303
Changes made in this Pull Request:
gemmilibrary (link) to parse mmcif filesclass MMCIFReader(base.SingleFrameReaderBase)andclass MMCIFParser(TopologyReaderBase)classes for thatAs a bonus, this implementation would potentially allow to read any of the gemmi-supported formats (source):
Also, this (with slight modifications) also would allow reading mmcif with multiple models sharing the same topology, as well as more feature-rich parsing of PDBs (the same code without changes can be used for parsing altlocs, charges, etc, from all of these formats).
However, I'm slightly lost on what's to be done next for this PR to be merged, so I'm asking if someone could help me navigate here (tagging @richardjgowers here as author of original PDBx implementation 4303).
PR Checklist
Developers certificate of origin
📚 Documentation preview 📚: https://mdanalysis--4712.org.readthedocs.build/en/4712/