The add_hydrogen_bond_interactions function assigns the 'atom_type' of all constructed hydrogen bonds as 'CA'

[wowaffle]

Describe the bug
The add_hydrogen_bond_interactions function assigns the atom_type of all constructed hydrogen bonds as CA, which is unexpected.

To Reproduce

1. Model a lab-synthesized peptide fragment. I'll paste the pdb raw file at the end of this issue.
2. Use Graphein to compute hydrogen bonds. The number and type of non-covalent interactions (e.g., hydrogen bonds, aromatic interactions) differ from what is shown in the PDB Mol* Viewer.
   According to Mol* Viewer, the peptide fragment should have 7 hydrogen bonds (including one S---H hydrogen bond).
3. I referred to issue #370 and tried using the modified code in that fork:
   Modified add_hydrogen_bond_interactions → hydrogen bond count = 0.
   Commenting out the final for loop in the modified function → hydrogen bonds reappear, but the count is higher than expected.

I suspected the issue might be related to donor/acceptor definitions. To debug, I printed atom_type for hydrogen bond edges in both the master branch and the modified version, codes are as below:

for r1, r2 in get_edges_by_bond_type(G, "hbond"):
    print(G.nodes[r1]["atom_type"], G.nodes[r2]["atom_type"])

Actual behavior
All hydrogen bond edges have atom_type = 'CA'.

Expected behavior
The atom_type should correspond to the correct atoms involved in hydrogen bonding (e.g., N, O, SG), not all 'CA'.

Additional info
Peptide fragment pdb file as below:

REMARK 1 CREATED WITH OPENMM 8.3, 2023-02-21
ATOM 1 N PRO A 1 200.011 221.987 299.505 1.00 0.00 N
ATOM 2 H2 PRO A 1 199.070 221.962 299.881 1.00 0.00 H
ATOM 3 H3 PRO A 1 200.035 222.706 298.793 1.00 0.00 H
ATOM 4 CA PRO A 1 200.972 222.339 300.565 1.00 0.00 C
ATOM 5 C PRO A 1 201.662 223.674 300.347 1.00 0.00 C
ATOM 6 O PRO A 1 202.075 224.310 301.324 1.00 0.00 O
ATOM 7 CB PRO A 1 201.968 221.170 300.540 1.00 0.00 C
ATOM 8 CG PRO A 1 201.769 220.502 299.226 1.00 0.00 C
ATOM 9 CD PRO A 1 200.333 220.695 298.878 1.00 0.00 C
ATOM 10 HA PRO A 1 200.525 222.353 301.426 1.00 0.00 H
ATOM 11 HB2 PRO A 1 202.873 221.511 300.617 1.00 0.00 H
ATOM 12 HB3 PRO A 1 201.772 220.559 301.267 1.00 0.00 H
ATOM 13 HG2 PRO A 1 202.338 220.920 298.561 1.00 0.00 H
ATOM 14 HG3 PRO A 1 201.976 219.558 299.306 1.00 0.00 H
ATOM 15 HD2 PRO A 1 200.225 220.736 297.918 1.00 0.00 H
ATOM 16 HD3 PRO A 1 199.794 219.988 299.266 1.00 0.00 H
ATOM 17 N ASP A 2 201.799 224.124 299.098 1.00 0.00 N
ATOM 18 CA ASP A 2 202.421 225.415 298.831 1.00 0.00 C
ATOM 19 C ASP A 2 201.622 226.586 299.391 1.00 0.00 C
ATOM 20 O ASP A 2 202.178 227.681 299.520 1.00 0.00 O
ATOM 21 CB ASP A 2 202.612 225.603 297.324 1.00 0.00 C
ATOM 22 CG ASP A 2 203.912 225.002 296.818 1.00 0.00 C
ATOM 23 OD1 ASP A 2 204.719 224.529 297.645 1.00 0.00 O
ATOM 24 OD2 ASP A 2 204.130 225.006 295.588 1.00 0.00 O
ATOM 25 H ASP A 2 201.542 223.701 298.395 1.00 0.00 H
ATOM 26 HA ASP A 2 203.294 225.434 299.252 1.00 0.00 H
ATOM 27 HB2 ASP A 2 201.880 225.172 296.856 1.00 0.00 H
ATOM 28 HB3 ASP A 2 202.622 226.552 297.123 1.00 0.00 H
ATOM 29 N ALA A 3 200.346 226.390 299.723 1.00 0.00 N
ATOM 30 CA ALA A 3 199.504 227.447 300.268 1.00 0.00 C
ATOM 31 C ALA A 3 199.259 227.296 301.764 1.00 0.00 C
ATOM 32 O ALA A 3 198.400 227.993 302.312 1.00 0.00 O
ATOM 33 CB ALA A 3 198.166 227.489 299.529 1.00 0.00 C
ATOM 34 H ALA A 3 199.940 225.637 299.638 1.00 0.00 H
ATOM 35 HA ALA A 3 199.943 228.300 300.128 1.00 0.00 H
ATOM 36 HB1 ALA A 3 197.573 228.108 299.983 1.00 0.00 H
ATOM 37 HB2 ALA A 3 197.778 226.600 299.526 1.00 0.00 H
ATOM 38 HB3 ALA A 3 198.319 227.787 298.619 1.00 0.00 H
ATOM 39 N LEU A 4 199.990 226.407 302.439 1.00 0.00 N
ATOM 40 CA LEU A 4 199.773 226.138 303.857 1.00 0.00 C
ATOM 41 C LEU A 4 201.002 226.459 304.701 1.00 0.00 C
ATOM 42 O LEU A 4 201.153 225.917 305.800 1.00 0.00 O
ATOM 43 CB LEU A 4 199.354 224.681 304.058 1.00 0.00 C
ATOM 44 CG LEU A 4 198.071 224.268 303.330 1.00 0.00 C
ATOM 45 CD1 LEU A 4 197.906 222.758 303.329 1.00 0.00 C
ATOM 46 CD2 LEU A 4 196.859 224.937 303.960 1.00 0.00 C
ATOM 47 H LEU A 4 200.625 225.942 302.093 1.00 0.00 H
ATOM 48 HA LEU A 4 199.048 226.697 304.176 1.00 0.00 H
ATOM 49 HB2 LEU A 4 200.068 224.107 303.739 1.00 0.00 H
ATOM 50 HB3 LEU A 4 199.213 224.528 305.006 1.00 0.00 H
ATOM 51 HG LEU A 4 198.125 224.561 302.407 1.00 0.00 H
ATOM 52 HD11 LEU A 4 198.674 222.357 302.893 1.00 0.00 H
ATOM 53 HD12 LEU A 4 197.846 222.447 304.246 1.00 0.00 H
ATOM 54 HD13 LEU A 4 197.095 222.530 302.848 1.00 0.00 H
ATOM 55 HD21 LEU A 4 196.921 225.895 303.825 1.00 0.00 H
ATOM 56 HD22 LEU A 4 196.846 224.737 304.909 1.00 0.00 H
ATOM 57 HD23 LEU A 4 196.055 224.595 303.540 1.00 0.00 H
ATOM 58 N PHE A 5 201.879 227.333 304.216 1.00 0.00 N
ATOM 59 CA PHE A 5 203.059 227.712 304.979 1.00 0.00 C
ATOM 60 C PHE A 5 202.664 228.496 306.226 1.00 0.00 C
ATOM 61 O PHE A 5 201.645 229.189 306.257 1.00 0.00 O
ATOM 62 CB PHE A 5 204.006 228.546 304.116 1.00 0.00 C
ATOM 63 CG PHE A 5 204.638 227.777 302.992 1.00 0.00 C
ATOM 64 CD1 PHE A 5 205.484 226.714 303.251 1.00 0.00 C
ATOM 65 CD2 PHE A 5 204.389 228.121 301.674 1.00 0.00 C
ATOM 66 CE1 PHE A 5 206.066 226.006 302.218 1.00 0.00 C
ATOM 67 CE2 PHE A 5 204.969 227.417 300.638 1.00 0.00 C
ATOM 68 CZ PHE A 5 205.808 226.360 300.910 1.00 0.00 C
ATOM 69 H PHE A 5 201.814 227.718 303.450 1.00 0.00 H
ATOM 70 HA PHE A 5 203.529 226.911 305.260 1.00 0.00 H
ATOM 71 HB2 PHE A 5 203.508 229.282 303.728 1.00 0.00 H
ATOM 72 HB3 PHE A 5 204.718 228.890 304.677 1.00 0.00 H
ATOM 73 HD1 PHE A 5 205.661 226.471 304.131 1.00 0.00 H
ATOM 74 HD2 PHE A 5 203.823 228.834 301.485 1.00 0.00 H
ATOM 75 HE1 PHE A 5 206.632 225.293 302.404 1.00 0.00 H
ATOM 76 HE2 PHE A 5 204.794 227.657 299.757 1.00 0.00 H
ATOM 77 HZ PHE A 5 206.200 225.884 300.214 1.00 0.00 H
ATOM 78 N GLN A 6 203.489 228.380 307.262 1.00 0.00 N
ATOM 79 CA GLN A 6 203.253 229.037 308.539 1.00 0.00 C
ATOM 80 C GLN A 6 204.297 230.119 308.779 1.00 0.00 C
ATOM 81 O GLN A 6 205.432 230.028 308.304 1.00 0.00 O
ATOM 82 CB GLN A 6 203.290 228.027 309.689 1.00 0.00 C
ATOM 83 CG GLN A 6 202.002 227.943 310.490 1.00 0.00 C
ATOM 84 CD GLN A 6 201.689 229.219 311.243 1.00 0.00 C
ATOM 85 OE1 GLN A 6 202.241 230.278 310.950 1.00 0.00 O
ATOM 86 NE2 GLN A 6 200.799 229.124 312.223 1.00 0.00 N
ATOM 87 H GLN A 6 204.210 227.912 307.248 1.00 0.00 H
ATOM 88 HA GLN A 6 202.378 229.454 308.531 1.00 0.00 H
ATOM 89 HB2 GLN A 6 203.468 227.147 309.324 1.00 0.00 H
ATOM 90 HB3 GLN A 6 204.001 228.278 310.300 1.00 0.00 H
ATOM 91 HG2 GLN A 6 201.265 227.762 309.886 1.00 0.00 H
ATOM 92 HG3 GLN A 6 202.081 227.225 311.139 1.00 0.00 H
ATOM 93 HE21 GLN A 6 200.434 228.366 312.399 1.00 0.00 H
ATOM 94 HE22 GLN A 6 200.587 229.820 312.681 1.00 0.00 H
ATOM 95 N GLY A 7 203.899 231.149 309.527 1.00 0.00 N
ATOM 96 CA GLY A 7 204.782 232.235 309.888 1.00 0.00 C
ATOM 97 C GLY A 7 205.422 232.033 311.251 1.00 0.00 C
ATOM 98 O GLY A 7 205.339 230.968 311.866 1.00 0.00 O
ATOM 99 H GLY A 7 203.103 231.235 309.841 1.00 0.00 H
ATOM 100 HA2 GLY A 7 205.490 232.306 309.231 1.00 0.00 H
ATOM 101 HA3 GLY A 7 204.285 233.067 309.899 1.00 0.00 H
ATOM 102 N ASN A 8 206.065 233.095 311.730 1.00 0.00 N
ATOM 103 CA ASN A 8 206.809 233.064 312.979 1.00 0.00 C
ATOM 104 C ASN A 8 206.216 234.059 313.965 1.00 0.00 C
ATOM 105 O ASN A 8 205.405 234.916 313.608 1.00 0.00 O
ATOM 106 CB ASN A 8 208.293 233.393 312.763 1.00 0.00 C
ATOM 107 CG ASN A 8 208.848 232.797 311.485 1.00 0.00 C
ATOM 108 OD1 ASN A 8 209.018 231.584 311.370 1.00 0.00 O
ATOM 109 ND2 ASN A 8 209.113 233.649 310.507 1.00 0.00 N
ATOM 110 H ASN A 8 206.084 233.861 311.339 1.00 0.00 H
ATOM 111 HA ASN A 8 206.747 232.179 313.371 1.00 0.00 H
ATOM 112 HB2 ASN A 8 208.398 234.356 312.718 1.00 0.00 H
ATOM 113 HB3 ASN A 8 208.806 233.038 313.506 1.00 0.00 H
ATOM 114 HD21 ASN A 8 208.978 234.491 310.617 1.00 0.00 H
ATOM 115 HD22 ASN A 8 209.423 233.354 309.764 1.00 0.00 H
ATOM 116 N THR A 9 206.633 233.923 315.222 1.00 0.00 N
ATOM 117 CA THR A 9 206.291 234.873 316.272 1.00 0.00 C
ATOM 118 C THR A 9 207.481 235.700 316.733 1.00 0.00 C
ATOM 119 O THR A 9 207.288 236.741 317.369 1.00 0.00 O
ATOM 120 CB THR A 9 205.695 234.143 317.482 1.00 0.00 C
ATOM 121 OG1 THR A 9 206.704 233.346 318.112 1.00 0.00 O
ATOM 122 CG2 THR A 9 204.534 233.258 317.059 1.00 0.00 C
ATOM 123 H THR A 9 207.125 233.273 315.496 1.00 0.00 H
ATOM 124 HA THR A 9 205.616 235.481 315.938 1.00 0.00 H
ATOM 125 HB THR A 9 205.362 234.796 318.118 1.00 0.00 H
ATOM 126 HG1 THR A 9 206.408 233.020 318.827 1.00 0.00 H
ATOM 127 HG21 THR A 9 204.088 232.900 317.843 1.00 0.00 H
ATOM 128 HG22 THR A 9 203.896 233.774 316.542 1.00 0.00 H
ATOM 129 HG23 THR A 9 204.857 232.522 316.517 1.00 0.00 H
ATOM 130 N GLY A 10 208.702 235.261 316.436 1.00 0.00 N
ATOM 131 CA GLY A 10 209.889 235.975 316.835 1.00 0.00 C
ATOM 132 C GLY A 10 211.046 235.764 315.881 1.00 0.00 C
ATOM 133 O GLY A 10 210.915 235.117 314.837 1.00 0.00 O
ATOM 134 H GLY A 10 208.862 234.538 315.998 1.00 0.00 H
ATOM 135 HA2 GLY A 10 209.697 236.925 316.877 1.00 0.00 H
ATOM 136 HA3 GLY A 10 210.163 235.676 317.716 1.00 0.00 H
ATOM 137 N PRO A 11 212.205 236.315 316.223 1.00 0.00 N
ATOM 138 CA PRO A 11 213.393 236.146 315.382 1.00 0.00 C
ATOM 139 C PRO A 11 214.023 234.772 315.572 1.00 0.00 C
ATOM 140 O PRO A 11 213.673 234.009 316.472 1.00 0.00 O
ATOM 141 CB PRO A 11 214.325 237.256 315.872 1.00 0.00 C
ATOM 142 CG PRO A 11 213.945 237.445 317.304 1.00 0.00 C
ATOM 143 CD PRO A 11 212.468 237.178 317.388 1.00 0.00 C
ATOM 144 HA PRO A 11 213.178 236.286 314.447 1.00 0.00 H
ATOM 145 HB2 PRO A 11 215.248 236.970 315.792 1.00 0.00 H
ATOM 146 HB3 PRO A 11 214.169 238.067 315.364 1.00 0.00 H
ATOM 147 HG2 PRO A 11 214.430 236.812 317.855 1.00 0.00 H
ATOM 148 HG3 PRO A 11 214.144 238.355 317.575 1.00 0.00 H
ATOM 149 HD2 PRO A 11 212.258 236.708 318.210 1.00 0.00 H
ATOM 150 HD3 PRO A 11 211.970 238.007 317.315 1.00 0.00 H
ATOM 151 OXT PRO A 11 214.892 234.544 314.676 1.00 0.00 O
TER 152 PRO A 11
END

Looking forward to ur reply, thank u~

[a-r-j]

Hi @wowaffle could you share a code snippet I can check out.

Also there are no sulfur-containing residues in the PDB file you've attached.

In any case, you can try the following:

from pathlib import Path
import networkx as nx
from graphein.protein.graphs import construct_graph
from graphein.protein.config import ProteinGraphConfig, GetContactsConfig
from graphein.protein.edges.intramolecular import _create_interaction_function

pdb_path = Path("./test.pdb")

# Create all interaction functions
hydrogen_bond_backbone_backbone = _create_interaction_function(
    "hbbb", "backbone to backbone hydrogen bond interactions"
)
hydrogen_bond_sidechain_backbone = _create_interaction_function(
    "hbsb", "side-chain to backbone hydrogen bond interactions"
)
hydrogen_bond_sidechain_sidechain = _create_interaction_function(
    "hbss", "side-chain to side-chain hydrogen bond interactions"
)


def hydrogen_bond(G: nx.Graph) -> nx.Graph:
    """
    Adds all types of hydrogen bonds to protein structure graph.
    Combines backbone-backbone, sidechain-backbone, and sidechain-sidechain hydrogen bonds.

    :param G: nx.Graph to add hydrogen bonds to.
    :return: nx.Graph with all hydrogen bond types added.
    """
    G = hydrogen_bond_backbone_backbone(G)
    G = hydrogen_bond_sidechain_backbone(G)
    G = hydrogen_bond_sidechain_sidechain(G)
    return G


getc = GetContactsConfig(
    get_contacts_path=Path("/Users/arian/github/getcontacts/"),
    contacts_dir=pdb_path.parent,
    pdb_dir=pdb_path.parent,
    granularity="atom",
)

cfg = ProteinGraphConfig(
    granularity="atom",
    exclude_waters=True,
    keep_hets=[],
    deprotonate=False,
    edge_construction_functions=[
        hydrogen_bond,
    ],
    node_metadata_functions=[],
    get_contacts_config=getc,
)

G = construct_graph(cfg, path=str(pdb_path))

print(G.edges(data=True))

Which outputs:

[('A:PRO:1:N', 'A:ASP:2:N', {'kind': {'hbbb'}, 'distance': np.float64(2.815905587824153)}), ('A:PRO:1:O', 'A:ALA:3:N', {'kind': {'hbbb'}, 'distance': np.float64(3.143104411097348)}), ('A:PRO:1:O', 'A:ASP:2:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2507574187824173)}), ('A:PRO:1:O', 'A:LEU:4:N', {'kind': {'hbbb'}, 'distance': np.float64(3.1603484147263448)}), ('A:ASP:2:N', 'A:ALA:3:N', {'kind': {'hbbb'}, 'distance': np.float64(2.7634449215649473)}), ('A:ASP:2:O', 'A:ALA:3:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2503597087160907)}), ('A:ALA:3:N', 'A:LEU:4:N', {'kind': {'hbbb'}, 'distance': np.float64(2.7392865208364094)}), ('A:ALA:3:O', 'A:LEU:4:N', {'kind': {'hbbb'}, 'distance': np.float64(2.249367261381758)}), ('A:LEU:4:N', 'A:PHE:5:N', {'kind': {'hbbb'}, 'distance': np.float64(2.7538184551437976)}), ('A:LEU:4:O', 'A:PHE:5:N', {'kind': {'hbbb'}, 'distance': np.float64(2.245237918685989)}), ('A:PHE:5:O', 'A:GLN:6:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2505143410929955)}), ('A:GLN:6:O', 'A:GLY:7:N', {'kind': {'hbbb'}, 'distance': np.float64(2.258878806552089)}), ('A:GLN:6:OE1', 'A:GLY:7:N', {'kind': {'hbsb'}, 'distance': np.float64(2.352140605132748)}), ('A:GLY:7:O', 'A:ASN:8:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2515971208621792)}), ('A:ASN:8:O', 'A:THR:9:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2580875341925406)}), ('A:THR:9:N', 'A:GLY:10:N', {'kind': {'hbbb'}, 'distance': np.float64(2.746787790390215)}), ('A:THR:9:O', 'A:GLY:10:N', {'kind': {'hbbb'}, 'distance': np.float64(2.2494989481824916)}), ('A:THR:9:OG1', 'A:GLY:10:N', {'kind': {'hbsb'}, 'distance': np.float64(3.2354636347424974)})

[wowaffle]

Thank you very much for your prompt reply.

Also, my apologies for providing an incorrect PDB file earlier. However, the issue is not specific to that file, so we can simply proceed with the previously shared PDB as an example.

Here is the code I originally used to identify hydrogen bonds in the peptide fragment:

import os
import networkx as nx
from pathlib import Path
from graphein.protein.graphs import construct_graph
from graphein.protein.config import ProteinGraphConfig, GetContactsConfig
from graphein.protein.edges.intramolecular import  hydrogen_bond
from graphein.protein.edges.distance import add_hydrogen_bond_interactions

def make_config(pdb_dir: Path) -> ProteinGraphConfig:
    getc = GetContactsConfig(
        get_contacts_path=Path("/mypathto/getcontacts"),
        contacts_dir=pdb_dir,
        pdb_dir=pdb_dir,
        granularity="CA"
    )

    edge_fns = [add_hydrogen_bond_interactions]
    #or using hb searching method in edges.intramolecular
    #edge_fns = [hydrogen_bond]

    cfg = ProteinGraphConfig(
        granularity="CA",
        keep_hets=[],
        exclude_waters=True,
        deprotonate=False,
        edge_construction_functions=edge_fns,
        node_metadata_functions=[],
        get_contacts_config=getc
    )
    return cfg

def build_graph_for_pdb(pdb_path: Path) -> nx.Graph:
    cfg = make_config(pdb_path.parent)
    G = construct_graph(cfg, path=str(pdb_path))
    three2one = {
        'ALA':'A','ARG':'R','ASN':'N','ASP':'D','CYS':'C','GLN':'Q','GLU':'E','GLY':'G',
        'HIS':'H','ILE':'I','LEU':'L','LYS':'K','MET':'M','PHE':'F','PRO':'P','SER':'S',
        'THR':'T','TRP':'W','TYR':'Y','VAL':'V', 'SEC':'U', 'PYL':'O'
    }
    for n, d in G.nodes(data=True):
        res3 = d.get("residue_name", "").upper()
        d["residue_type"] = three2one.get(res3, "X")
    return G

def print_graph_summary(G: nx.Graph):
    from collections import Counter
    kinds = []
    for u, v, d in G.edges(data=True):
        k = d.get("kind")
        if isinstance(k, (set, list, tuple)):
            for kk in k: kinds.append(str(kk))
        elif k is not None:
            kinds.append(str(k))
    cnt = Counter(kinds)
    print(f"Nodes: {G.number_of_nodes()}, Edges: {G.number_of_edges()}")
    print("Edge kinds (counts):", dict(cnt))

#main
G = build_graph_for_pdb("./test.pdb")
print_graph_summary(G)

Which outputs (using pdb file provided as above):

Nodes: 11, Edges: 10
Edge kinds (counts): {'peptide_bond': 10}

I tried both edge.intramolecular and edge.distance, but as you can see, no hb were found in this protein fragment. From my understanding, the interactions found by intramolecular should be fewer than those from distance, since the rules for intramolecular are presumably stricter. But in this case, they should at least found hb existed.

Another point is that I successfully reproduced your code and obtained the same results as you did. However, in the PDB Mol* Viewer, this protein only shows two hydrogen bonds:

('A:PRO:1:O', 'A:LEU:4:N', {'kind': {'hbbb'}, 'distance': np.float64(3.1603484147263448)})
('A:GLN:6:OE1', 'A:GLY:7:N', {'kind': {'hbsb'}, 'distance': np.float64(2.352140605132748)})

You can also load the previously mentioned PDB file into the PDB Mol* Viewer. There, you will see that it only considers two hydrogen bonds in the fragment, exactly corresponding to the two dashed interactions shown in the figure below.

I would like to hear your perspective on which is more reliable for identifying non-covalent interactions in such cases: PDB Mol* Viewer or Graphein?

I apologize again that my initial question did not include sufficient information. I hope you might be able to take another look and kindly share your thoughts. Thank you very much for your time and help ;)