search

MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.32k stars 652 forks source link

PDB not importing all bonds #4100

Closed jamesdalg closed 1 year ago

jamesdalg commented 1 year ago

Expected behavior

All of the bonds in molecule are imported into the PDB object. At the moment, I'm importing the human polymerase, but I'm getting very few of the bonds to actually come through. I would imagine that, based on the fact alone that it is made of proteins... there should be more than 35 or 36 bonds in the entirety of a molecule with 32892 or 32712 atoms. It's a bit odd... is this typical?

Actual behavior

PDB files for 7OL0 and 5FLM import with less than 40 bonds in total.

Code to reproduce the behavior

(base) jd@oban:/code/python/reprex-3-29-2023$ pip install --upgrade MDAnalysis
Requirement already satisfied: MDAnalysis in /home/jd/.local/lib/python3.10/site-packages (2.3.0)
Collecting MDAnalysis
  Using cached MDAnalysis-2.4.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.0 MB)
Requirement already satisfied: packaging in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (21.3)
Requirement already satisfied: biopython>=1.80 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.80)
Requirement already satisfied: joblib>=0.12 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.2.0)
Requirement already satisfied: tqdm>=4.43.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (4.64.1)
Requirement already satisfied: gsd>=1.9.3 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (2.6.1)
Requirement already satisfied: matplotlib>=1.5.1 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (3.6.2)
Requirement already satisfied: GridDataFormats>=0.4.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.0.1)
Requirement already satisfied: fasteners in /home/jd/miniconda3/lib/python3.10/site-packages (from MDAnalysis) (0.18)
Requirement already satisfied: numpy>=1.20.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.23.4)
Requirement already satisfied: networkx>=2.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (2.8.8)
Requirement already satisfied: mmtf-python>=1.0.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.1.3)
Requirement already satisfied: scipy>=1.5.0 in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (1.9.3)
Requirement already satisfied: threadpoolctl in /home/jd/.local/lib/python3.10/site-packages (from MDAnalysis) (3.1.0)
Requirement already satisfied: mrcfile in /home/jd/.local/lib/python3.10/site-packages (from GridDataFormats>=0.4.0->MDAnalysis) (1.4.3)
Requirement already satisfied: fonttools>=4.22.0 in /home/jd/.local/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (4.38.0)
Requirement already satisfied: pyparsing>=2.2.1 in /home/jd/miniconda3/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (3.0.9)
Requirement already satisfied: contourpy>=1.0.1 in /home/jd/.local/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (1.0.6)
Requirement already satisfied: cycler>=0.10 in /home/jd/.local/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (0.11.0)
Requirement already satisfied: pillow>=6.2.0 in /home/jd/miniconda3/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (9.4.0)
Requirement already satisfied: python-dateutil>=2.7 in /home/jd/miniconda3/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (2.8.2)
Requirement already satisfied: kiwisolver>=1.0.1 in /home/jd/.local/lib/python3.10/site-packages (from matplotlib>=1.5.1->MDAnalysis) (1.4.4)
Requirement already satisfied: msgpack>=1.0.0 in /home/jd/.local/lib/python3.10/site-packages (from mmtf-python>=1.0.0->MDAnalysis) (1.0.4)
Requirement already satisfied: six>=1.5 in /home/jd/miniconda3/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib>=1.5.1->MDAnalysis) (1.16.0)
Installing collected packages: MDAnalysis
  Attempting uninstall: MDAnalysis
    Found existing installation: MDAnalysis 2.3.0
    Uninstalling MDAnalysis-2.3.0:
      Successfully uninstalled MDAnalysis-2.3.0
Successfully installed MDAnalysis-2.4.2
(base) jd@oban:/code/python/reprex-3-29-2023$ python
Python 3.10.8 (main, Nov 24 2022, 14:13:03) [GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import MDAnalysis as mda
>>> pdb_u_7Ol0_nb=mda.Universe('7OL0.pdb')
/home/jd/miniconda3/lib/python3.10/site-packages/MDAnalysis/coordinates/PDB.py:432: UserWarning: 1 A^3 CRYST1 record, this is usually a placeholder. Unit cell dimensions will be set to None.
  warnings.warn("1 A^3 CRYST1 record,"
>>> pdb_u_5flm_nb=mda.Universe('5FLM.pdb')
>>> pdb_u_7Ol0_nb.bonds
<TopologyGroup containing 36 bonds>
>>> pdb_u_7Ol0_nb.atoms
<AtomGroup with 32892 atoms>
>>> pdb_u_5flm_nb.bonds
<TopologyGroup containing 35 bonds>
>>> pdb_u_5flm_nb.atoms
<AtomGroup with 32712 atoms>
>>> pdb_u_5flm_nb.bonds.atom1.resnames
array(['CYS', 'CYS', 'CYS', 'HIS', 'CYS', 'CYS', 'CYS', 'ASP', 'ASP',
       'ASP', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS',
       'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS',
       'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'C'], dtype=object)
>>> pdb_u_7Ol0_nb.bonds.atom1.resnames
array(['CYS', 'CYS', 'CYS', 'HIS', 'CYS', 'CYS', 'CYS', 'CYS', 'ASP',
       'ASP', 'ASP', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS',
       'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS',
       'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'CYS', 'A'],
      dtype=object)
>>> pdb_u_5flm_nb.bonds.atom1.resnums
array([  71,   74,   81,   84,  111,  114,  184,  495,  499,  499, 1119,
       1122, 1137, 1140,   88,   90,   94,   97,   17,   20,   39,   42,
         86,   89,  114,  119,    7,   10,   44,   45,   19,   22,   36,
         39,   20])
>>> pdb_u_7Ol0_nb.bonds.atom1.resnums
array([  71,   74,   81,   84,  111,  114,  154,  184,  495,  497,  499,
       1196, 1199, 1214, 1217,   88,   90,   94,   97,   17,   20,   39,
         42,   86,   89,  114,  119,    7,   10,   44,   45,   19,   22,
         36,   39,   46])
>>> 

....

Current version of MDAnalysis

jamesdalg commented 1 year ago

The sequence is known for all chains. Is there some way to fix this with the fasta file (which gives the amino acid sequence)?

>5FLM_2|Chain B|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB2|BOS TAURUS (9913)
MYDADEDMQYDEDDDEITPDLWQEACWIVISSYFDEKGLVRQQLDSFDEFIQMSVQRIVEDAPPIDLQAEAQHASGEVEEPPRYLLKFEQIYLSKPTHWERDGAPSPMMPNEARLRNLTYSAPLYVDITKTVIKEGEEQLQTQHQKTFIGKIPIMLRSTYCLLNGLTDRDLCELNECPLDPGGYFIINGSEKVLIAQEKMATNTVYVFAKKDSKYAYTGECRSCLENSSRPTSTIWVSMLARGGQGAKKSAIGQRIVATLPYIKQEVPIIIVFRALGFVSDRDILEHIIYDFEDPEMMEMVKPSLDEAFVIQEQNVALNFIGSRGAKPGVTKEKRIKYAKEVLQKEMLPHVGVSDFCETKKAYFLGYMVHRLLLAALGRRELDDRDHYGNKRLDLAGPLLAFLFRGMFKNLLKEVRIYAQKFIDRGKDFNLELAIKTRIISDGLKYSLATGNWGDQKKAHQARAGVSQVLNRLTFASTLSHLRRLNSPIGRDGKLAKPRQLHNTLWGMVCPAETPEGHAVGLVKNLALMAYISVGSQPSPILEFLEEWSMENLEEISPAAIADATKIFVNGCWVGIHKDPEQLMNTLRKLRRQMDIIVSEVSMIRDIREREIRIYTDAGRICRPLLIVEKQKLLLKKRHIDQLKEREYNNYSWQDLVASGVVEYIDTLEEETVMLAMTPDDLQEKEVAYCSTYTHCEIHPSMILGVCASIIPFPDHNQSPRNTYQSAMGKQAMGVYITNFHVRMDTLAHVLYYPQKPLVTTRSMEYLRFRELPAGINSIVAIASYTGYNQEDSVIMNRSAVDRGFFRSVFYRSYKEQESKKGFDQEEVFEKPTRETCQGMRHAIYDKLDDDGLIAPGVRVSGDDVIIGKTVTLPENEDELEGTNRRYTKRDCSTFLRTSETGIVDQVMVTLNQEGYKFCKIRVRSVRIPQIGDKFASRHGQKGTCGIQYRQEDMPFTCEGITPDIIINPHAIPSRMTIGHLIECLQGKVSANKGEIGDATPFNDAVNVQKISNLLSDYGYHLRGNEVLYNGFTGRKITSQIFIGPTYYQRLKHMVDDKIHSRARGPIQILNRQPMEGRSRDGGLRFGEMERDCQIAHGAAQFLRERLFEASDPYQVHVCNLCGIMAIANTRTHTYECRGCRNKTQISLVRMPYACKLLFQELMSMSIAPRMMSV
>5FLM_14|Chain N[auth P]|RNA, DNA-RNA ELONGATION SCAFFOLD|SYNTHETIC CONSTRUCT (32630)
UAUAUGCAUAAAGACCAGGC
>5FLM_1|Chain A|DNA-DIRECTED RNA POLYMERASE|BOS TAURUS (9913)
MHGGGPPSGDSACPLRTIKRVQFGVLSPDELKRMSVTEGGIKYPETTEGGRPKLGGLMDPRQGVIERTGRCQTCAGNMTECPGHFGHIELAKPVFHVGFLVKTMKVLRCVCFFCSKLLVDSNNPKIKDILAKSKGQPKKRLTHVYDLCKGKNICEGGEEMDNKFGVEQPEGDEDLTKEKGHGGCGRYQPRIRRSGLELYAEWKHVNEDSQEKKILLSPERVHEIFKRISDEECFVLGMEPRYARPEWMIVTVLPVPPLSVRPAVVMQGSARNQDDLTHKLADIVKINNQLRRNEQNGAAAHVIAEDVKLLQFHVATMVDNELPGLPRAMQKSGRPLKSLKQRLKGKEGRVRGNLMGKRVDFSARTVITPDPNLSIDQVGVPRSIAANMTFAEIVTPFNIDRLQELVRRGNSQYPGAKYIIRDNGDRIDLRFHPKPSDLHLQTGYKVERHMCDGDIVIFNRQPTLHKMSMMGHRVRILPWSTFRLNLSVTTPYNADFDGDEMNLHLPQSLETRAEIQELAMVPRMIVTPQSNRPVMGIVQDTLTAVRKFTKRDVFLERGEVMNLLMFLSTWDGKVPQPAILKPRPLWTGKQIFSLIIPGHINCIRTHSTHPDDEDSGPYKHISPGDTKVVVENGELIMGILCKKSLGTSAGSLVHISYLEMGHDITRLFYSNIQTVINNWLLIEGHTIGIGDSIADSKTYQDIQNTIKKAKQDVIEVIEKAHNNELEPTPGNTLRQTFENQVNRILNDARDKTGSSAQKSLSEYNNFKSMVVSGAKGSKINISQVIAVVGQQNVEGKRIPFGFKHRTLPHFIKDDYGPESRGFVENSYLAGLTPTEFFFHAMGGREGLIDTAVKTAETGYIQRRLIKSMESVMVKYDATVRNSINQVVQLRYGEDGLAGESVEFQNLATLKPSNKAFEKKFRFDYTNERALRRTLQEDLVKDVLSNAHIQNELEREFERMREDREVLRVIFPTGDSKVVLPCNLLRMIWNAQKIFHINPRLPSDLHPIKVVEGVKELSKKLVIVNGDDPLSRQAQENATLLFNIHLRSTLCSRRMAEEFRLSGEAFDWLLGEIESKFNQAIAHPGEMVGALAAQSLGEPATQMTLNTFHYAGVSAKNVTLGVPRLKELINISKKPKTPSLTVFLLGQSARDAERAKDILCRLEHTTLRKVTANTAIYYDPNPQSTVVAEDQEWVNVYYEMPDFDVARISPWLLRVELDRKHMTDRKLTMEQIAEKINAGFGDDLNCIFNDDNAEKLVLRIRIMNSDENKMQEEEEVVDKMDDDVFLRCIESNMLTDMTLQGIEQISKVYMHLPQTDNKKKIIITEDGEFKALQEWILETDGVSLMRVLSEKDVDPVRTTSNDIVEIFTVLGIEAVRKALERELYHVISFDGSYVNYRHLALLCDTMTCRGHLMAITRHGVNRQDTGPLMKCSFEETVDVLMEAAAHGESDPMKGVSENIMLGQLAPAGTGCFDLLLDAEKCKYGMEIPTNIPGLGAAGPTGMFFGSAPSPMGGISPAMTPWNQGATPAYGAWSPSVGSGMTPGAAGFSPSAASDASGFSPGYSPAWSPTPGSPGSPGPSSPYIPSPGGAMSPSYSPTSPAYEPRSPGGYTPQSPSYSPTSPSYSPTSPSYSPTSPNYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPSYSPTSPNYSPTSPNYTPTSPSYSPTSPSYSPTSPNYTPTSPNYSPTSPSYSPTSPSYSPTSPSYSPSSPRYTPQSPTYTPSSPSYSPSSPSYSPTSPKYTPTSPSYSPSSPEYTPTSPKYSPTSPKYSPTSPKYSPTSPTYSPTTPKYSPTSPTYSPTSPVYTPTSPKYSPTSPTYSPTSPKYSPTSPTYSPTSPKGSTYSPTSPGYSPTSPTYSLTSPAISPDDSDDEN
>5FLM_10|Chain J|DNA-DIRECTED RNA POLYMERASES I, II, AND III SUBUNIT RPABC5|BOS TAURUS (9913)
MIIPVRCFTCGKIVGNKWEAYLGLLQAEYTEGDALDALGLKRYCCRRMLLAHVDLIEKLLNYAPLEK
>5FLM_11|Chain K|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB11|BOS TAURUS (9913)
MNAPPAFESFLLFEGEKKITINKDTKVPNACLFTINKEDHTLGNIIKSQLLKDPQVLFAGYKVPHPLEHKIIIRVQTTPDYSPQEAFTNAITDLISELSLLEERFRVAIKDKQEGIE
>5FLM_12|Chain L|DNA-DIRECTED RNA POLYMERASES I, II, AND III SUBUNIT RPABC4|BOS TAURUS (9913)
MDTQKDVQPPKQQPMIYICGECHTENEIKSRDPIRCRECGYRIMYKKRTKRLVVFDAR
>5FLM_13|Chain M[auth N]|DNA, DNA-RNA ELONGATION SCAFFOLD|SYNTHETIC CONSTRUCT (32630)
GGCAGTACTAGTAAACTAGTATTGAAAGTACTTGAGCTT
>5FLM_15|Chain O[auth T]|DNA, DNA-RNA ELONGATION SCAFFOLD|SYNTHETIC CONSTRUCT (32630)
AAGCTCAAGTACTTAAGCCTGGTCATTACTAGTACTGCC
>5FLM_3|Chain C|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB3|BOS TAURUS (9913)
MPYANQPTVRITELTDENVKFIIENTDLAVANSIRRVFIAEVPIIAIDWVQIDANSSVLHDEFIAHRLGLIPLTSDDIVDKLQYSRDCTCEEFCPECSVEFTLDVRCNEDQTRHVTSRDLISNSPRVIPVTSRNRDNDPNDYVEQDDILIVKLRKGQELRLRAYAKKGFGKEHAKWNPTAGVAFEYDPDNALRHTVYPKPEEWPKSEYSELDEDESQAPYDPNGKPERFYYNVESCGSLRPETIVLSALSGLKKKLSDLQTQLSHEIQSDVLTIN
>5FLM_4|Chain D|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB4|BOS TAURUS (9913)
MAAGGSDPRSGDVEEDASQLIFPKEFETAETLLNSEVHMLLEHRKQQNESAEDEQELSEVFMKTLNYTARFSRFKNRETIASVRSLLLQKKLHKFELACLANLCPETAEESKALIPSLEGRFEDEELQQILDDIQTKRSFQY
>5FLM_5|Chain E|DNA-DIRECTED RNA POLYMERASES I, II, AND III SUBUNIT RPABC1|BOS TAURUS (9913)
MDDEEETYRLWKIRKTIMQLCHDRGYLVTQDGLDQTLEEFKAQFGGKPSEGRPRRTDLTVLVAHNDDPTDQMFVFFPEEPKVGIKTIKVYCQRMQEENITRALIVVQQGMTPSAKQSLVDMAPKYILEQFLQQELLINITEHELVPEHVVMTKEEVTELLARYKLRENQLPRIQAGDPVARYFGIKRGQVVKIIRPSETAGRYITYRLVQ
>5FLM_6|Chain F|DNA-DIRECTED RNA POLYMERASES I, II, AND III SUBUNIT RPABC2|BOS TAURUS (9913)
MSDNEDNFDGDDFDDVEEDEGLDDLENAEEEGQENVEILPSGERPQANQKRITTPYMTKYERARVLGTRALQIAMCAPVMVELEGETDPLLIAMKELKARKIPIIIRRYLPDGSYEDWGVDELIITD
>5FLM_7|Chain G|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB7|BOS TAURUS (9913)
MFYHISLEHEILLHPRYFGPNLLNTVKQKLFTEVEGTCTGKYGFVIAVTTIDNIGAGVIQPGRGFVLYPVKYKAIVFRPFKGEVVDAVVTQVNKVGLFTEIGPMSCFISRHSIPSEMEFDPNSNPPCYKTMDEDIVIQQDDEIRLKIVGTRVDKNDIFAIGSLMDDYLGLVS
>5FLM_8|Chain H|DNA-DIRECTED RNA POLYMERASES I, II, AND III SUBUNIT RPABC3|BOS TAURUS (9913)
MAGILFEDIFDVKDIDPEGKKFDRVSRLHCESESFKMDLILDVNIQIYPVDLGDKFRLVIASTLYEDGTLDDGEYNPTDDRPSRADQFEYVMYGKVYRIEGDETSTEAATRLSAYVSYGGLLMRLQGDANNLHGFEVDSRVYLLMKKLAF
>5FLM_9|Chain I|DNA-DIRECTED RNA POLYMERASE II SUBUNIT RPB9|BOS TAURUS (9913)
MEPDGTYEPGIVGIRFCQECNNMLYPKEDKENRILLYACRNCDYQQEADNSCIYVNKITHEVDELTQIIADVSQDPTLPRTEDHPCQKCGHKEAVFFQSHSARAEDAMRLYYVCTAPHCGHRWTE
richardjgowers commented 1 year ago

hi @jamesdalg I suspect it's only reading the CONECT records, which it looks like you've got for your sulphurs. If you do mda.Universe(..., guess_bonds=True) it will guess bonds based on coordinates in addition to the CONECT records, this might get what you expect.

IAlibay commented 1 year ago

Just to add here, the current bond guessing is a little bit dumb (looks at things like the distance between atoms), indeed one could improve this by using the expected bonds from known residues. There is currently work in progress to do this which should hopefully make it to a release of MDAnalysis soon (see: https://github.com/MDAnalysis/mdanalysis/pull/3866).

I am closing this issue as resolved for now, please do re-open if you think this needs further discussions / there are other underlying issues not currently addressed.

jamesdalg commented 1 year ago

Thanks! This was insightful and helpful.

jamesdalg commented 1 year ago

@IAlibay If you get this fixed, let me know. right now, it draws connections between the DNA and the tail of the polymerase. If you know of a decent workaround that would be helplful to know as well. Is there some format that will have the topology correct from the get-go that I could convert the PDB to? What do you do in your simulations? I've tried XPDB, GSD (which has the topology correct, but drops the residues), FHIAIMS, and several others. Is there a conversion tool you'd recommend that could work to convert it to something that I could have the correct topology and residues in MDAnalysis?

IAlibay commented 1 year ago

@jamesdalg - am I correct in understanding that you are getting these weird bonds after guessing them?

One quick option might be to guess on atomgroups for parts of the system that you know are contiguous, i.e. protein = u.select_atoms('protein') protein.guess_bonds() dna = u.select_atoms('resnames DA DC DG DI DT DU') dna.guess_bonds()

I think that should avoid guessing bonds between the contiguous polymers.

Alternatively re: bond containing files, you'd have to use "fully informative" molecular dynamics topologies (e.g. Gromacs TPR, and AMBER PARM7). What is your current entry point / simulation method for your files?

jamesdalg commented 1 year ago

@IAlibay I think I'll try the guess_bonds on individual segments. That's the most sensible way of going about it. Here's my approach: I'm using MDAnalysis to pull in a PDB from https://www.rcsb.org/structure/5flm, from there, I'm attempting to coarse grain it, taking the average x,y,z coordinates and plugging them into the positions in a universe object. From there, I use the existing, unique bonds between residues (using pandas with merges on a table of atoms and a table of bonds) and then add them to a new universe where each residue has one particle per residue, then output the file to LAMMPS. I could easily try it in GROMACS or CHARMM, which is what's nice about MDAnalysis-- once all the pieces are put in correctly, you can export it anywhere. I think I'll try to convert the PDB to gromacs first and then import it from there. I'll try AMBER too.

Great suggestions! Thanks!