scripts/xplor2gmx.pl

   1 #!/usr/bin/perl -w
   2
   3 #
   4 # This script reads an XPLOR input file with distance restraint data
   5 # as sometimes is found in the pdb database (http://www.pdb.org).
   6 # From this input file dihedral restrints should be removed, such that
   7 # only distance restraints are left. The script can handle ambiguous restraints.
   8 # It converts the distance restraints to GROMACS format.
   9 # A typical command line would be
  10 # ./xplor2gmx.pl 33 conf.pdb < restraints.dat > disre.itp
  11 # You can turn off debugging info, but please do verify your output is correct!
  12 #
  13 # David van der Spoel (spoel@gromacs.org), July 2002
  14 #
  15
  16 # Turn debugging off (0), or on ( > 0).
  17 $debug = 1;
  18 # Turn atom name translation on and off
  19 $trans = 1;
  20
  21 $res0  = shift;# || die "I need the residue offset\n";
  22 $pdb   = shift || die "I need the name of the pdb file with correct atom numbers\n";
  23 $core  = shift || "core.ndx";
  24 $tbl   = "/home/spoel/gmxdev/share/top/atom_nom.tbl";
  25
  26 printf "[ distance_restraints ]\n";
  27 printf "; Read an xplor distance restraint file, and output GROMACS distance restraints.\n";
  28 printf "; This also needs a pdb file with correct GROMACS atom numbering.\n";
  29 printf "; I used $pdb for the atom numbers\n";
  30 printf "; This also needs the offset in residues.\n";
  31 printf "; I used $res0 for the residue offset\n";
  32
  33 # Read the pdb file
  34 # If things go wrong, check whether your pdb file is read correctly.
  35 $natom = 0;
  36 $nres  = 0;
  37 open(PDB,$pdb) || die "Can not open file $pdb\n";
  38 while ($line = <PDB>) {
  39     if (index($line,"ATOM") >= 0) {
  40         @tmp = split(' ',$line);
  41         $aname[$natom] = $tmp[2];
  42         $resnr[$natom] = $tmp[4];
  43         if ($resnr[$natom] > $nres) {
  44             $rname[$nres] = $tmp[3];
  45             $nres++;
  46         }
  47         $natom++;
  48     }
  49 }
  50 close PDB;
  51 printf "; I found $natom atoms in the pdb file $pdb\n";
  52 printf "; I found $nres residues in the pdb file $pdb\n";
  53 if ($debug > 1) {
  54     for ($i = 0; ($i < $natom); $i ++) {
  55         printf("; %5d  %5s  %5s  %5d\n",$i+1,$aname[$i],
  56                $rname[$resnr[$i]-1],$resnr[$i]);
  57     }
  58 }
  59
  60 #
  61 # Read the name translation table.
  62 # Source for this comes from: http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl
  63 # It was adapted slightly for GROMACS use, but only as regards formatting,
  64 # not for content
  65 #
  66 open(TBL,$tbl) || die "Can not open atom-name translation table $tbl\n";
  67 $ntbl=0;
  68 while ($line = <TBL>) {
  69     @ttt = split('#',$line);
  70     @tmp = split(' ',$ttt[0]);
  71     if ($#tmp == 2) {
  72         # New table entry
  73         $tblres[$ntbl] = $tmp[0];
  74         $tblxplor[$ntbl] = $tmp[1];
  75         $tblgmx[$ntbl] = $tmp[2];
  76         $ntbl++;
  77     }
  78 }
  79 close TBL;
  80 printf "; Read $ntbl entries from $tbl\n";
  81
  82 @templates = (
  83  [ "HA#", "HA1", "HA2" ],
  84  [ "HA*", "HA1", "HA2" ],
  85  [ "HB#",  "HB",         "HB1", "HB2"   ],
  86  [ "HB*",  "HB",         "HB1", "HB2"   ],
  87  [ "HG#",  "HG",         "HG1", "HG2",  "HG11", "HG12", "HG13", "HG21", "HG22", "HG23"  ],
  88  [ "HG*",  "HG",         "HG1", "HG2",  "HG11", "HG12", "HG13", "HG21", "HG22", "HG23"  ],
  89  [ "HG1#", "HG11",      "HG12", "HG13"  ],
  90  [ "HG1*", "HG11",      "HG12", "HG13"  ],
  91  [ "HG2#", "HG21",      "HG22", "HG23"  ],
  92  [ "HG2*", "HG21",      "HG22", "HG23"  ],
  93  [ "HD#",  "HD1",       "HD2",  "HD11", "HD12", "HD13", "HD21", "HD22", "HD23" ],
  94  [ "HD*",  "HD1",       "HD2",  "HD11", "HD12", "HD13", "HD21", "HD22", "HD23" ],
  95  [ "HD1#", "HD11",      "HD12"  ],
  96  [ "HD1*", "HD11",      "HD12"  ],
  97  [ "HD2#", "HD21",      "HD22"  ],
  98  [ "HD2*", "HD21",      "HD22"  ],
  99  [ "HE#",  "HE",        "HE1",  "HE2"   ],
 100  [ "HE*",  "HE",        "HE1",  "HE2"   ],
 101  [ "HH#",  "HH11",      "HH12", "HH21", "HH22" ],
 102  [ "HH*",  "HH11",      "HH12", "HH21", "HH22" ],
 103  [ "HZ#",  "HZ",         "HZ1", "HZ2",  "HZ3"   ],
 104  [ "HZ*",  "HZ",         "HZ1", "HZ2",  "HZ3"   ],
 105  [ "HN",   "H" ]
 106 );
 107
 108 $ntranslated = 0;
 109 $nuntransltd = 0;
 110 sub transl_aname {
 111     my $atom   = shift(@_);
 112     my $resnr  = shift(@_);
 113
 114     if ( $trans == 1 ) {
 115         for(my $i = 0; ($i < $ntbl); $i ++) {
 116             if ($tblres[$i] eq $rname[$resnr]) {
 117                 if ($tblxplor[$i] eq $atom) {
 118                     $ntranslated++;
 119                     return $tblgmx[$i];
 120                 }
 121             }
 122         }
 123     }
 124     $nuntransltd++;
 125     if ($debug > 1) {
 126         printf "; No name change for $rname[$resnr] $atom\n";
 127     }
 128     return $atom;
 129 }
 130
 131 sub expand_template {
 132     my $atom  = shift(@_);
 133     my $rnum  = shift(@_)-1;
 134     my $bdone = 0;
 135     my @atoms;
 136     my $jj = 0;
 137
 138     for (my $tt=0; (($tt <= $#templates) && ($bdone == 0)); $tt++) {
 139         $templ = $templates[$tt];
 140         if ($atom eq $templ->[0]) {
 141             for ($jj = 1; ($jj <= $#{$templ}); $jj++) {
 142                 push @atoms, transl_aname($templ->[$jj],$rnum);
 143             }
 144             $bdone  = 1;
 145         }
 146     }
 147     if ($bdone == 0) {
 148         push @atoms, transl_aname($atom,$rnum);
 149     }
 150     if ($debug > 0) {
 151         my $natom = $#atoms+1;
 152         printf("; Found $natom elements for atom $rname[$rnum] %d $atom:",
 153                $rnum+1);
 154         for $aa ( @atoms ) {
 155             printf " $aa";
 156         }
 157         printf "\n";
 158     }
 159     return @atoms;
 160 }
 161
 162 if ($debug > 1) {
 163     printf "; There are $#templates template entries\n";
 164     for ($tt=0; ($tt <= $#templates); $tt++) {
 165         $templ  = $templates[$tt];
 166         $ntempl = $#{$templ};
 167         printf "; Item $tt ($templates[$tt][0]) has $ntempl entries\n";
 168     }
 169 }
 170
 171 # This index file holds numbers of restraints involving core atoms
 172 @protcore = ( "H", "HN", "HA", "HA1", "HA2", "HB", "HB1", "HB2", "HB3", "HG", "HG1", "HG2", "HG3", "N", "O"  );
 173 open(CORE,">$core") || die "Can not open $core\n";
 174 print CORE "[ core-restraints ]\n";
 175 $ncore = 0;
 176
 177 $myindex     = 0;
 178 $linec       = 0;
 179 $npair       = 0;
 180 $nunresolved = 0;
 181 while ($line = <STDIN>) {
 182     @ttt = split('!',$line);
 183     if ((index($ttt[0],"assign") >= 0) && (index($ttt[0],"!assign") < 0)) {
 184         @tmp = split('\(',$ttt[0]);
 185         # Find first argument
 186         if (($rhaak  = index($tmp[1],')')) < 0) {
 187             printf "No ) in '$tmp[1]'\n";
 188         }
 189         $atres1 = substr($tmp[1],0,$rhaak);
 190         @at1 = split('or',$atres1);
 191
 192         # Find second argument
 193         if (($rhaak  = index($tmp[2],')')) < 0) {
 194             printf "No ) in '$tmp[2]'\n";
 195         }
 196         $atres2 = substr($tmp[2],0,$rhaak);
 197         @at2 = split('or',$atres2);
 198
 199         @expdata = split('\)',$tmp[2]);
 200         @dist    = split(' ',$expdata[1]);
 201
 202         $bOK   = 0;
 203         $bCore = 1;
 204         foreach $a1 ( @at1 ) {
 205             @info1  = split(' ',$a1);
 206             $r1     = $info1[1] - $res0;
 207             @atoms1 = expand_template($info1[4],$r1);
 208
 209             foreach $a2 ( @at2 ) {
 210                 @info2  = split(' ',$a2);
 211                 $r2     = $info2[1] - $res0;
 212                 @atoms2 = expand_template($info2[4],$r2);
 213
 214                 for ($i = 0; ($i < $natom) && ($resnr[$i] < $r1); $i++) { ; }
 215                 for ( ; ($i < $natom) && ($resnr[$i] == $r1); $i++) {
 216                     foreach $ii ( @atoms1 ) {
 217                         if ($ii eq $aname[$i]) {
 218                             $bCoreI = 0;
 219                             for $pp ( @protcore ) {
 220                                 if ($ii eq $pp) {
 221                                     $bCoreI = 1;
 222                                 }
 223                             }
 224
 225                             for ($j = 0; ($j < $natom) && ($resnr[$j] < $r2); $j++) { ; }
 226                             for ( ; ($j < $natom) && ($resnr[$j] == $r2); $j++) {
 227                                 foreach $jj ( @atoms2 ) {
 228                                     if ($jj eq $aname[$j]) {
 229                                         $dd     = 0.1*$dist[0];
 230                                         $dminus = 0.1*$dist[1];
 231                                         $dplus  = 0.1*$dist[2];
 232                                         $low    = $dd-$dminus;
 233                                         $up1    = $dd+$dplus;
 234                                         $up2    = $up1+1;
 235                                         printf("%5d %5d %5d %5d %5d %7.3f %7.3f %7.3f 1.0; res $r1 $ii - res $r2 $jj\n",$i+1,$j+1,1,$myindex,1,$low,$up1,$up2);
 236                                         # Do some checks
 237                                         $bOK    = 1;
 238                                         $bCoreJ = 0;
 239                                         for $pp ( @protcore ) {
 240                                             if ($jj eq $pp) {
 241                                                 $bCoreJ = 1;
 242                                             }
 243                                         }
 244                                         if (($bCoreI == 0) || ($bCoreJ == 0)) {
 245                                             if ($debug > 0) {
 246                                                 printf "; No core $ii ($bCoreI) $jj ($bCoreJ)\n";
 247                                             }
 248                                             $bCore = 0;
 249                                         }
 250                                         $npair++;
 251                                     }
 252                                 }
 253                             }
 254                         }
 255                     }
 256                 }
 257             }
 258         }
 259         if ($bCore == 1) {
 260             printf CORE "$myindex\n";
 261             $ncore++;
 262         }
 263         if ($bOK == 0) {
 264             printf "; UNRESOLVED: $ttt[0]\n";
 265             $nunresolved++;
 266         }
 267         else {
 268             $myindex++;
 269         }
 270     }
 271     $linec++;
 272 }
 273 close CORE;
 274
 275 printf "; A total of $myindex lines with distance restraints were read from $linec input lines\n";
 276 printf "; From this, $npair actual restraints were generated.\n";
 277 printf "; A total of $nunresolved lines could not be interpreted\n";
 278 printf "; In the process $ntranslated atoms had a name change\n";
 279 printf "; However for $nuntransltd no name change could be found\n";
 280 printf "; usually these are either HN->H renamings or not-existing atoms\n";
 281 printf "; generated by the template expansion (e.g. HG#)\n";
 282 printf "; A total of $ncore restraints were classified as core restraints\n";