Re: [ccp4bb]: surface residue average B

[James Irving]

***  For details on how to be removed from this list visit the  ***
***          CCP4 home page http://www.ccp4.ac.uk         ***


Hi Bernhard,

The following perlscript will take the output from the DSSPcmbi
(http://swift.cmbi.ru.nl/gv/dssp/) command-line program, picks
residues according to a minimum absolute (rather than relative)
accessibility, and extracts B values and occupancy from a PDB file.

Cheers,
James

#!/usr/bin/perl

#usage:
# avg_surf.pl dsspin structure.dssp pdbin structure.pdb acc 20

#mainchain atoms; all others are defined as sidechain
$main="C N CA O OXT OXT1 OXT2";

#process command line arguments
@[EMAIL PROTECTED];
for ($i=0; $i<@tmp; $i+=2) { $tmp[$i]=~tr/A-Z/a-z/; }
print join(";",$tmp);
[EMAIL PROTECTED];

#read dsspcmbi file
open(I,$cmd{dsspin}) ||
 die "DSSP file open failed; use: dsspin file.dssp\n";

#minimum accessibility must be input
if ($cmd{acc} !~/\d/) {
 die "Minimum accessibility not defined; use: acc 20\n";
}

#find column header line
$acpos=-1; $respos=-1;
while (($_=<I>)&&($acpos==-1)) {
 if ($_=~/RESIDUE.+AA.+STRUCTURE/) {
   $acpos=index($_,"ACC");
   $respos=index($_,"RESIDUE");
 }
}

#identify residues >= accepted surface accessibility
if (($acpos != -1)&&($respos != -1)) {
 while ($_=<I>) {
   if (substr($_,$acpos,3) >= $cmd{acc}) {
     $tmp=substr($_,$respos,7);
     $tmp=~s/\s+//g;
     $resok.="$tmp ";
   }
 }
}

close(I);

#open and read pdb file
open(I,$cmd{pdbin}) ||
 die "PDB file open failed; use: pdbin file.pdb\n";

while ($_=<I>) {
 if ($_=~/^ATOM/) {
   $res=substr($_,22,4).substr($_,21,1);
   $res=~s/\s+//g;
   $occ=substr($_,54,6);
   $bf=substr($_,60,6);
   #ensure atom is from a shortlisted residue
   if ($resok=~/\b$res\b/) {
     $atm=substr($_,12,4);
     $atm=~s/\s+//g;
     #if defined as mainchain...
     if ($main=~/\b$atm\b/) {
       $mainav+=$bf*$occ;
       $mainn+=$occ;
     } else {
     #else assumed to be sidechain
       $sideav+=$bf*$occ;
       $siden+=$occ;
     }
   }
 }
}

close(I);

$mainout = ($mainn>0) ? "$mainn atoms with average B of ".
 sprintf("%0.2f",$mainav/$mainn) : "none";
$sideout = ($siden>0) ? "$siden atoms with average B of ".
 sprintf("%0.2f",$sideav/$siden) : "none";

print STDERR "Files: $cmd{dsspin} $cmd{pdbin}
Minimum accessibility: $cmd{acc}
Mainchain: $mainout
Sidechain: $sideout\n";