Contact map between protein residues and non-protein residues

[Former user]

Dear bio3d users, I am able to make the contact map for protein residues and for glycan residues (with selection "notprotein" ) seperately. Is the anyway that I can have the contact map between theses two groups of residues? Regards, Maryam

[Barry Grant]

Yes Maryam. You can use the combine.select() function to get at what you want. The example code below walks through this for using the cmap() function together with a pdb class object input.

> library(bio3d)

## Read an example PDB file
> pdb <- read.pdb("4q21")
  Note: Accessing on-line PDB file
  HEADER    ONCOGENE PROTEIN                        25-SEP-91   4Q21

## Note there are two non-protein 'ligands' in this structure (GDP and MG, as well as 78 water molecules).
> pdb

 Call:  read.pdb(file = "4q21")

   Total Models#: 1
     Total Atoms#: 1447,  XYZs#: 4341  Chains#: 1  (values: A)

     Protein Atoms#: 1340  (residues/Calpha atoms#: 168)
     Nucleic acid Atoms#: 0  (residues/phosphate atoms#: 0)

     Non-protein/nucleic Atoms#: 107  (residues: 80)
     Non-protein/nucleic resid values: [ GDP (1), HOH (78), MG (1) ]

   Protein sequence:
      MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSYRKQVVIDGETCLLDILDTAG
      QEEYSAMRDQYMRTGEGFLCVFAINNTKSFEDIHQYREQIKRVKDSDDVPMVLVGNKCDL
      AARTVESRQAQDLARSYGIPYIETSAKTRQGVEDAFYTLVREIRQHKL

+ attr: atom, helix, sheet, seqres, xyz,
        calpha, remark, call

# Note also that there are 168 Calpha atoms (see above) but actually 169 protein residues as the last ARG does not have a Calpha.
> head( tail(pdb$atom, 110) )
       type eleno elety  alt resid chain resno insert      x      y      z o
1338   ATOM  1338   CD1 <NA>   LEU     A   168   <NA> 49.504 55.704 27.166 1
1339   ATOM  1339   CD2 <NA>   LEU     A   168   <NA> 51.598 55.879 25.836 1
1340   ATOM  1340     N <NA>   ARG     A   169   <NA> 49.999 60.160 22.815 1
1341 HETATM  1342    MG <NA>    MG     A   273   <NA> 65.614 76.977 46.715 1
1342 HETATM  1343    PB <NA>   GDP     A   274   <NA> 62.667 77.781 47.505 1
1343 HETATM  1344   O1B <NA>   GDP     A   274   <NA> 61.587 77.413 46.626 1
         b segid elesy charge
1338 49.58  <NA>     C   <NA>
1339 50.30  <NA>     C   <NA>
1340 53.67  <NA>     N   <NA>
1341 20.79  <NA>    MG   <NA>
1342 31.08  <NA>     P   <NA>
1343 30.69  <NA>     O   <NA>


# A regular contact map function call will include these two ligands entries and all protein residues (i.e. 169+2=171)
> cm <- cmap(pdb)
> dim(cm)
[1] 171 171  


# Select protein atoms only
> prot.ind <- atom.select(pdb, "protein")
> cm.prot <- cmap(pdb, inds=prot.ind)
> dim(cm.prot)
[1] 169 169  


# Lets add the MG ion ligand to our previous protein selection
> mg.ind <- atom.select(pdb, resid="MG")
> inds <- combine.select(prot.ind, mg.ind, operator="OR")
 Union of selects
 *  Selected a total of: 1341 atoms  *

# Contact map
> cm.mg <- cmap(pdb, inds=inds)
> dim(cm.mg)
[1] 170 170

Hope this helps, Barry