Csi2Aria.py
Csi2Aria.py — Python Source, 74 KB (76071 bytes)
File contents
#!/usr/local/bin/perl # Set tab spacing to 4 when viewing this file! ###################### # AUTHOR INFORMATION # ###################### # Problems, questions or suggestions? Contact the author: # # Wim Vranken (wim@bri.nrc.ca) # Biomolecular NMR # Biotechnology Research Institute # Montreal # Canada ########## # USAGE: # ########## # aria.overview.pl <options> # # Start this script in your <ARIAproject>/run?? directory. # # Use aria.overview.pl -h for list of <options> ########## # OUTPUT # ########## # Description of output files (default for <output> is 'overview') # # <output>.info Info on unambiguous restraints long-range restraints that initially # had an ambiguous short-range possibility # Info on number of unambiguous long-range restraints between residues # Info on not well supported unambiguous long-range restraints. # Info on number of intra/seq/medium/long range contacts per spectrum # # <output).newunambig New unambiguous assignments made by ARIA # # <output).newambig New ambiguous assignments made by ARIA # # <output>.info.detail Detailed info on long-range possibilities in ambiguous restraints. # # <output>.aria Overview of aria output (pdb, violations, ...) # # <output>.ps Postscript file displaying condensed restraint information # # <output>.restr (optional) Info on long/med/seq/intra contacts per spectrum & per residue # # <output>.pdb.xmgr XMGR (now GRACE) graph file displaying information on PDB files # <output>.restraints.xmgr XMGR (now GRACE) graph file displaying information on restraints ################# # MODIFICATIONS # ################# # Kevin Gardner 26/10/2001 - two modifications for more flexible .tbl reading. # Version 1.1 (06/02/2001) # Includes modifications to 1.0 by A. Bonvin (18/01/2001) # (can be set with -nw (print new unamb/amb assignments) option) # Xmgr graphs more universal format # Prints number of intra/seq/medium/long range contacts per spectrum to .info # # Version 1.2 (13/06/2001) # Added short violations list sorted by ppm value to .info file # Added -pr option to print file with number of restraints per residue # Fixed 'Odd number of elements in hash list' warning on SGI machines. # Patched count number of intra/seq/med/long based on unambig/ambig, # program now also distinguishes between these types. # Contact map now generated (as default) from unambig.tbl and ambig.tbl. # Can use old way (using spectrum data) with -sp on flag. # Fixed additional ',[,]' printing problem in .info file ############################################################# # FOLLOWING ARE OPTIONS FOR PROGRAM # ############################################################# $helptext .=" -h print help -f <name> output file name (DEFAULT 'overview') -n xx number of iterations to check (DEFAULT 8) -ns xx start number for iterations to check (DEFAULT 0) -d on/off print details on long-range ambiguous assignments (DEFAULT OFF) -c on/off print calibration info -nw on/off print list new unambiguous/unambiguous assignments -vl on/off print consistent violations for LAST iteration -va on/off print consistent violations for ALL iterations (DEFAULT OFF & overrides -vl) -vc on/off use 'compressed' strong violations format -r on/off print restraint overviews -rc on/off print spectra restraint overviews -p on/off print pdb file energy overview -pr on/off print number of restraints per residues to .restr file (DEFAULT OFF) -np xx number of lowest overall energy pdb files for overview -sp on/off use spectra data instead of ambig/unambig.tbl for generation of contact map (DEFAULT OFF) Default is 'on' unless specified. "; ############################################################# # Main routines # ------------- # Initialize program &initialize; # Get command flag data &getflags; # Open output files &openoutput; # Get directory names, specify number of iterations. &getiterationdirs($it_total,$it_start); # Get spectrum names. &getspecinfo; # Get info for restraint set in data directory &getorigrestr; # Get information for each iteration &getiterationinfo; # Print violation and calibration info per spectrum. &printspectruminfo; # Print restraint overview (conditional) &printrestraintoverview; # Print pdb file information overview (conditional) &printpdboverview; # Print xmgr graphs (conditional) &printxmgrgraphs; # Detailed reports &printdetails; # Calculate matrices &calculatematrices; # Report new unambiguous assignments (conditional) &reportunambig; # Report new ambiguous assignments (conditional) &reportambig; # Report 'suspect' contacts based on contact map &reportsuspects; # Print ps graph. &printps; # The end print "\n\nDone!\n"; ################################ # Main routines in above order # ################################ sub initialize { # Initialize some program variables # Elements found ARIA output pdb file headers (in order of appearance)... @PDBinfolist=('E(overall)','E(bonds)','E(angles)','E(improper)','E(vdw)','E(noe)', 'E(cdih)','E(coup)','E(sani)','','rms(bonds)','rms(angles)', 'rms(improper)','rms(noe)','rms(cdih)','rms(coup)','rms(sani)', '','viol(noe)','viol(cdih)','viol(coup)','viol(sani)'); # XMGR (now GRACE) graph elements that are displayed @Egraphlist=('E(overall)','E(noe)','E(vdw)','E(cdih)'); @rmsgraphlist=('rms(angles)','rms(improper)','rms(noe)','rms(cdih)'); @violgraphlist=('viol(noe)','viol(cdih)'); # Methyl names for atom name compression. $methyls="ALA HB|ILE HG2|ILE HD1|LEU HD1|LEU HD2|MET HE|THR HG2|VAL HG1|VAL HG2"; # For title printing @itname=('Initial iteration','First iteration','Second iteration','Third iteration','Fourth iteration','Fifth iteration','Sixth iteration','Seventh iteration','Eight iteration','Ninth iteration','Tenth iteration'); # Defaults for program run $outputfile="overview"; $it_start=0; $it_total=8; $pdbbest=10; $violprintall = 'off'; $detailprint = "off"; $newasnprint = 'on'; $printrestperres = 'off'; $graph='on'; $specmatrix='off'; $matrixtype='restraints'; # Identify script name, run date and user ($scriptname,$scriptname_short)=&getscriptname; ($date_time,$date)=&gettime; ($user)=&iduser; # Get project name and run number ($projectname,$run)=&getprojectinfo; } sub getscriptname { my($sn,$sns); # Get scriptname with and without path name $sn=$sns=$0; $sn=~ s/\/\//\//g; $sns=~ s/\/.+\///g; return $sn,$sns; } sub getflags { for ($a=0;$a<=$#ARGV;$a+=2) { if ($ARGV[$a] eq '-h') { &print_help; exit; } elsif ($ARGV[$a] eq '-f') { $outputfile=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-d' && $ARGV[$a+1]=~ /on|off/) { $detailprint=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-ns' && $ARGV[$a+1]=~ /\d+/) { $it_start=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-n' && $ARGV[$a+1]=~ /\d+/) { $it_total=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-c' && $ARGV[$a+1]=~ /on|off/) { $calibprint=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-nw' && $ARGV[$a+1]=~ /on|off/) { $newasnprint=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-vl' && $ARGV[$a+1]=~ /on|off/) { $violprintlast=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-va' && $ARGV[$a+1]=~ /on|off/) { $violprintall=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-vc' && $ARGV[$a+1]=~ /on|off/) { $violcompress=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-r' && $ARGV[$a+1]=~ /on|off/) { $restrover=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-rc' && $ARGV[$a+1]=~ /on|off/) { $restrover_spec=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-p' && $ARGV[$a+1]=~ /on|off/) { $pdbover=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-np' && $ARGV[$a+1]=~ /\d+/) { $pdbbest=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-pr' && $ARGV[$a+1]=~ /on|off/) { $printrestperres=$ARGV[$a+1]; } elsif ($ARGV[$a] eq '-sp' && $ARGV[$a+1]=~ /on|off/) { $specmatrix=$ARGV[$a+1]; if ($specmatrix eq 'on') { $matrixtype="spectra"; } } else { &print_help; print "Unrecognized option: $ARGV[$a] with $ARGV[$a+1]. Please check above options...\n"; exit; } } } sub openoutput { # Open output file(s) &openoutputfile("$outputfile.aria",OUTPUT,"$scriptname called by $user on $date.\n\n"); &openoutputfile("$outputfile.info",REPORT,"$scriptname called by $user on $date.\n\n"); if ($newasnprint eq 'on') { &openoutputfile("$outputfile.newunambig",NEWUNAMBIG,"$scriptname called by $user on $date.\n\n"); &openoutputfile("$outputfile.newambig",NEWAMBIG,"$scriptname called by $user on $date.\n\n"); } if ($detailprint eq 'on') { &openoutputfile("$outputfile.info.detail",DETAIL,"$scriptname called by $user on $date.\n\n"); } if ($printrestperres eq 'on') { &openoutputfile("$outputfile.restr",RESTPERRES,"$scriptname called by $user on $date.\n\n"); } } sub getiterationdirs { # Gets the existing iteration directories # so that program still works if some already deleted. my ($it_nums,$it_start) = @_; my ($it,$it_hit) = 0; for ($it=$it_start;$it<=($it_start+$it_nums);$it++) { $status=opendir (DIR, "structures/it$it"); if ($status == 0) {&send_warning("Not using iteration $it: $!.\n",STDOUT);next;} @dirlist = readdir DIR ; closedir DIR; if ($#dirlist <= 2) { &send_warning("Directory for iteration $it empty.\n",STDOUT); next; } else { # Check if ambig and unambig.tbl files generated - if not, skip directory @ambig = grep(/ambig\.tbl/, @dirlist); if ($#ambig < 1) { &send_warning("No ambig.tbl and unambig.tbl file for iteration $it.\n",STDOUT); next; } } $it_dir[$it_hit]=$it;$it_hit++; } # Check if final_data directory in last iteration $status=opendir (DIR, "structures/it".$it_dir[$#it_dir]."/finaldata"); if ($status == 0) { &send_warning("No finaldata directory found for iteration ".$it_dir[$#it_dir].".\n",STDOUT); $it_final=$it_dir[$#it_dir]; } else { @dirlist = readdir DIR ; closedir DIR; if ($#dirlist < 5) { &send_warning("Finaldata directory iteration ".$it_dir[$#it_dir]." empty.\n\n",STDOUT); $it_final=$it_dir[$#it_dir]; } else { $it_final=$it_dir[$#it_dir]+1; $it_dir[$#it_dir+1]=$it_final; $itname[$it_final]="Final data"; } } if (!defined @it_dir) {&send_fatal("No iterations found. Make sure you are in an aria 'run' directory!\n");} } sub getspecinfo { my ($spectot,$i,$j,$it)=0; my (@dirlist,@spec_list_temp) =""; print ("Getting spectrum names...\n\n"); $status=open (RUNCNS,"run.cns"); if ($status == 1) { while (defined ($_=<RUNCNS>)) { if (/^\{===>\} aspectrum_/) { s/^\{===>\} aspectrum_\d//; s/_pro//; s/=|\"|\;\s*//g;chomp; if ($_ ne "") { push(@spec_list_temp,$_); } } elsif (/{===>} assignstruct_/) { s/{===>} assignstruct_//;s/_pro//;s/ //g;chomp; $it=$_;$it =~ s/=.*//; s/\d+=//g;chop; $strongviol_structurestotal[$it]=$_; $strongviol_structurestotal[$it+1]=$_; } } } else { &send_warning("run.cns file not found: $!.\n",STDOUT); # What follows is not a very quick routine, but does make sure no # spectrum is missed (if already deleted somewhere). for ($i=0;$i<=$#it_dir;$i++) { $it=$it_dir[$i]; if ($itname[$it] !~ /Final/) { $filetext="it$it"; } else { $filetext="it".($it-1)."/finaldata"; } opendir (DIR, "structures/$filetext"); @dirlist = grep /^calib_\w+_it\d+.out$|\.tbl$/, readdir DIR ; closedir DIR; while ($j <= $#dirlist) { if ($dirlist[$j] =~ /ambig.tbl$/) { splice(@dirlist,$j,1);next; } elsif ($dirlist[$j] =~ /.tbl$/) { $dirlist[$j] =~ s/.tbl$//; } elsif ($dirlist[$j] =~ /^calib_\w+_it\d+.out$/) { $dirlist[$j] =~ s/^calib\_//; $dirlist[$j] =~ s/_it$it.out$//; } if ((grep(/^$dirlist[$j]$/,@spec_list_temp))==0) { push(@spec_list_temp,$dirlist[$j]); } $j++; } } } close RUNCNS; @spec_list = sort @spec_list_temp; } sub getorigrestr { # Gets the information of the original restraint list(s) from the data directory # Note that %assignment is not used in rest of program for now. my ($status,$origres_num) =0; my (@peakinfo) = ""; print ("\n"); foreach $spectrum (@spec_list) { $origres_num = 0; $status=open (ORIGRES,"data/$spectrum/$spectrum.tbl"); if ($status == 0) {&send_warning("No original .tbl file for spectrum $spectrum: $!.\n",STDOUT);next;} print ("Reading data/$spectrum/$spectrum.tbl...\n"); while (defined ($_=<ORIGRES>)) { # Skip comment lines next if (/\s*\!/); chomp; # Modification by Kevin Gardner 26/10/2001. if (/^\s*assi/ || /^\s*ASSI/) { $origres_num++;$peakbox3=0;$ppm3=0; $assignment{$spectrum}[$origres_num]=$_."\n"; # read original info if (/resid/ || /\{/) { $peakbox1=0; next; } s/assi \( attr store1 \<//; s/and attr store1 \>//; ($high_ppm,$low_ppm,@void)=split(' '); $peakbox1=($high_ppm-$low_ppm); next; } elsif (/^\s*OR/) { $assignment{$spectrum}[$origres_num].=$_."\n"; } elsif (/^\s*\(/) { $assignment{$spectrum}[$origres_num].=$_."\n"; if (/resid/ || /\{/) { $peakbox2=0; next; } s/^\s+\( attr store1 \<//; s/and attr store1 \>//; ($high_ppm,$low_ppm,@void)=split(' '); $peakbox2=($high_ppm-$low_ppm); next; } elsif (/^\s+and/) { $assignment{$spectrum}[$origres_num].=$_."\n"; s/^\s+and bondedto\s+\( attr store1 \<//; s/and attr store1 \>//; ($high_ppm,$low_ppm,@void)=split(' '); $peakbox3=($high_ppm-$low_ppm); $ppm3=($high_ppm+$low_ppm)/2; next; } elsif (/ peak /) { $peakinfo{$spectrum}[$origres_num]=$_; if (/weight/) { ( /peak\s+(\S+)/ ) && ( $peaknum=$1 ); ( /volume\s+(\S+)/ ) && ( $vol=$1 ); ( /ppm1\s+(\S+)/ ) && ( $ppm1=$1 ); ( /ppm2\s+(\S+)/ ) && ( $ppm2=$1 ); #This block Copyright 2002 Christian Gunning and Ryan Richter #Does the following line still work? Above handles change better. ($v0,$v1,$v2,$v3,$peaknum,$v5,$v6,$v7,$vol,$v9,$ppm1,$v11,$ppm2,@void)=split(' '); } else { ( /peak\s+(\S+)/ ) && ( $peaknum=$1 ); ( /volume\s+(\S+)/ ) && ( $vol=$1 ); ( /ppm1\s+(\S+)/ ) && ( $ppm1=$1 ); ( /ppm2\s+(\S+)/ ) && ( $ppm2=$1 ); #This block Copyright 2002 Christian Gunning and Ryan Richter #The following line doesn't handle an nmrview setup: #($v0,$v1,$v2,$v3,$peaknum,$v5,$vol,$v7,$ppm1,$v9,$ppm2)=split(' '); } $origres_ppm1{$spectrum}[$peaknum]=$ppm1; $origres_ppm2{$spectrum}[$peaknum]=$ppm2; $origres_ppm3{$spectrum}[$peaknum]=$ppm3; $origlistid{$spectrum}[$origres_num]=$peaknum; $origres_vol{$spectrum}[$peaknum]=$vol; $origres_pb1{$spectrum}[$peaknum]=$peakbox1; $origres_pb2{$spectrum}[$peaknum]=$peakbox2; $origres_pb3{$spectrum}[$peaknum]=$peakbox3; $origres_listid{$spectrum}[$peaknum]=$origres_num; #$peakinspec{$peaknum.','.$ppm1.','.$ppm2}=$spectrum; next; } } } } sub getiterationinfo { # Go through the 'normal' iterations. Print errors if files not found. for ($it=0;$it<=$#it_dir;$it++) { $it_cur=$it_dir[$it]; if ($itname[$it_cur] !~ /Final/ && $pdbover ne 'off') { $pdbread[$it_cur]=&getpdbinfo($it_cur); } &readitfiles($it_cur); } } sub printspectruminfo { print "\nWriting spectrum overviews...\n"; foreach $spectrum (@spec_list) { &printcalibinfo($spectrum) if ($calibprint ne 'off'); &printviolationinfo($spectrum) if ($violprintall ne 'off' || $violprintlast ne 'off'); } } sub printrestraintoverview { if ($restrover ne 'off') { print "\nWriting restraint overviews...\n"; &print_title("Restraint information overview",OUTPUT); print OUTPUT ("\n * n/u = restraints not used compared to iteration 0"); print OUTPUT ("\n * ratio = ambiguous restraints/unambiguous restraints"); print OUTPUT ("\n * avg_amb = ambiguous possibilities/ambiguous restraints\n"); for ($it=0;$it<=$#it_dir;$it++) { $it_cur=$it_dir[$it]; if ($itname[$it_cur] =~ /Final/) { $it_text="Final data"; } else {$it_text="Iteration $it_cur";} &printrestrinfo($it_cur,$it_text); &printrestrinfo_spec($it_cur) if ($restrover_spec ne 'off'); } } } sub printpdboverview { if ($pdbover ne 'off') { print "\nWriting PDB overviews...\n"; &print_title("PDB file energies overview",OUTPUT); print OUTPUT ("\n * format: minimum ... <average> ... maximum\n"); for ($it=0;$it<=$#it_dir;$it++) { $it_cur=$it_dir[$it]; next if ($pdbread[$it_cur] == 0 || $itname[$it_cur] =~ /Final/); &printpdbinfo($it_cur); } } } sub printxmgrgraphs { if ($graph eq 'on') { print "\nWriting XMGR graphs...\n"; &makegraph('restraints'); &makegraph('pdb'); } } sub printdetails { print "\nMaking initial matrices...\n"; # Note that the reportdetails subs are necessary for matrix generation! &print_title("Overview of type of restraints .",REPORT); printf REPORT ("\n(ambiguous restraint possibilities divided by total number possibilities for restraint)\n"); foreach $spectrum (@spec_list) { &reportdetails($spectrum,$it_dir[0],'off'); &reportdetails($spectrum,$it_final,$detailprint); } printf REPORT ("\n"); if ($printrestperres eq 'on') { &printrestperrestotal; } } sub calculatematrices { print "\nMaking dependent matrices...\n"; # Make and print 'dependent' contact map. $contact_avg[$it_final]=&makedependentmatrix ($it_final,$contact_all_sum[$it_final],$contact_all_max[$it_final],@contact_all); if ($it_final != $it_dir[0]) { $contact_avg[$it_dir[0]]=&makedependentmatrix ($it_dir[0],$contact_all_sum[$it_dir[0]],$contact_all_max[$it_dir[0]],@contact_all); } } sub reportsuspects { print "\nChecking for suspicious assignments...\n"; # Prints suspected contacts. foreach $spectrum (@spec_list) { &print_title("Spectrum: $spectrum",REPORT); &reportbadass($spectrum,$it_final); } } sub reportunambig{ if ($newasnprint eq 'on') { print "\nChecking for unambiguous assignments...\n"; # Prints unambiguous contacts. foreach $spectrum (@spec_list) { &print_title("Spectrum: $spectrum",NEWUNAMBIG); &reportunambass($spectrum,$it_final); } } } sub reportambig{ if ($newasnprint eq 'on') { print "\nChecking for ambiguous assignments...\n"; # Prints ambiguous contacts. foreach $spectrum (@spec_list) { &print_title("Spectrum: $spectrum",NEWAMBIG); &reportambass($spectrum,$it_final); } } } sub printps { print "\nWriting postscript file...\n"; &print_pspage ($outputfile.".ps","Helvetica",$it_dir[0],$itname[$it_dir[0]],$it_final,$itname[$it_final]); } #################################### # Subroutines for printing data... # #################################### sub printrestperrestotal { my (@type_list) = ('long','med','seq','intra'); print RESTPERRES "All spectra combined\n\n "; foreach $type (@type_list) { printf RESTPERRES ("%8s",$type); } for $res (1...$num_residues) { $line = sprintf ("\n%3d",$res); foreach $type (@type_list) { $line.=sprintf(" %5.1f ",$contact_info_all{$type}[$res]); } print RESTPERRES $line; } print RESTPERRES"\n\n"; } sub printpdbinfo { my ($it)=@_; print OUTPUT ("\nIteration $it (on $pdbcount[$it]/$pdball[$it] structures): \n\n"); foreach $name (@PDBinfolist) { if ($name eq '') {print OUTPUT "\n";next;} ($min,$avg,$max)=&minavgmax($pdb{$name}[$it]); if ($avg == 0) { print OUTPUT (" *All values zero for $name.\n"); next; } ($sd)=&sd($avg,$pdb{$name}[$it]); if ($name=~ /^E/) { printf OUTPUT (" *%-15s: %.2f ... <%.2f> ... %.2f [+/- %.2f]\n", $name,$min,$avg,$max,$sd); } elsif ($name=~ /^rms/) { printf OUTPUT (" *%-13s: %.3f ... <%.3f> ... %.3f [+/- %.3f]\n", $name,$min,$avg,$max,$sd); } elsif ($name=~ /^viol/) { printf OUTPUT (" *%-12s: %d ... <%.1f> ... %d\n", $name,$min,$avg,$max); } } print OUTPUT ("\n"); } sub printrestrinfo { my ($it,$ittext)=@_; my ($ratio,$avg_amb)=0; print OUTPUT ("\n$ittext: \n\n"); if ($restr[$it] == 0) {print OUTPUT " No general info available.\n";return;} if ($unambig_restr[$it] != 0) {$ratio=$ambig_restr[$it]/$unambig_restr[$it];} else {$ratio=0;} if ($ambig_restr[$it] != 0) {$avg_amb=$ambig_restr_poss[$it]/$ambig_restr[$it];} else {$avg_amb=0;} if ($restr_tot[0] != 0) {$del=$restr_tot[0]-$ambig_restr[$it]-$unambig_restr[$it];} else {$del='n/a'}; print OUTPUT (" Unambig ambig n/u ratio ambig_poss avg_amb\n\n"); printf OUTPUT (" %5d %5d %4s) %3.2f %5d %3.1f (%-25s\n\n", $unambig_restr[$it],$ambig_restr[$it],"(".$del, $ratio,$ambig_restr_poss[$it],$avg_amb,"Final)"); printf OUTPUT (" (unamb intra %.1f, seq %.1f, med %.1f, long %.1f)\n",$restrcnt{"intraunamb"}[$it], $restrcnt{"sequnamb"}[$it],$restrcnt{"medunamb"}[$it],$restrcnt{"longunamb"}[$it]); printf OUTPUT (" (ambig intra %.1f, seq %.1f, med %.1f, long %.1f)\n",$restrcnt{"intraamb"}[$it], $restrcnt{"seqamb"}[$it],$restrcnt{"medamb"}[$it],$restrcnt{"longamb"}[$it]); printf OUTPUT (" (all intra %.1f, seq %.1f, med %.1f, long %.1f)\n",$restrcnt{"intra"}[$it], $restrcnt{"seq"}[$it],$restrcnt{"med"}[$it],$restrcnt{"long"}[$it]); if ($graph eq 'on') { $unambig_rest_list.=$it." $unambig_restr[$it]\n"; $ambig_rest_list.=$it." $ambig_restr[$it]\n"; $all_rest_list.=$it." ".($ambig_restr[$it]+$unambig_restr[$it])."\n"; $avg_amb_list.=$it." $avg_amb\n"; $ratio_list.=$it." $ratio\n"; } } sub printrestrinfo_spec { my ($it)=@_; my ($ratio,$avg_amb)=0; print OUTPUT ("\n"); foreach $spectrum (@spec_list) { if ($restr{$spectrum}[$it] == 0) {print OUTPUT (" No info available ($spectrum)\n");next;} if ($unambig_restr{$spectrum}[$it] != 0) {$ratio=$ambig_restr{$spectrum}[$it]/$unambig_restr{$spectrum}[$it];} else {$ratio=0;} if ($ambig_restr{$spectrum}[$it] != 0) {$avg_amb=$ambig_restr_poss{$spectrum}[$it]/$ambig_restr{$spectrum}[$it];} else {$avg_amb=0;} if ($restr_tot{$spectrum}[0] != 0) {$del=$restr_tot{$spectrum}[0]-$ambig_restr{$spectrum}[$it]-$unambig_restr{$spectrum}[$it];} else {$del='n/a'}; printf OUTPUT (" %5d %5d %4s) %3.2f %5d %3.1f (%-25s\n", $unambig_restr{$spectrum}[$it], $ambig_restr{$spectrum}[$it],"(".$del, $ratio,$ambig_restr_poss{$spectrum}[$it], $avg_amb,$spectrum.")"); if ($graph eq 'on') { $unambig_rest_list{$spectrum}.=$it." $unambig_restr{$spectrum}[$it]\n"; $ambig_rest_list{$spectrum}.=$it." $ambig_restr{$spectrum}[$it]\n"; $all_rest_list{$spectrum}.=$it." ".($ambig_restr{$spectrum}[$it]+$unambig_restr{$spectrum}[$it])."\n"; $avg_amb_list{$spectrum}.=$it." $avg_amb\n"; $ratio_list{$spectrum}.=$it." $ratio\n"; } } print OUTPUT ("\n"); } sub printcalibinfo { my ($spectrum)=@_; &print_title("Calibration analysis for spectrum $spectrum",OUTPUT); for ($it=0;$it<=$#it_dir;$it++) { $it_cur=$it_dir[$it]; if ($itname[$it_cur] =~ /Final/) { $it_text="Final data"; } else {$it_text="Iteration $it";} printf OUTPUT ("\n$it_text: %3e %3e.", $calibinfo{$spectrum}[$it_cur][0],$calibinfo{$spectrum}[$it_cur][1]); } print OUTPUT "\n"; } sub printviolationinfo { my ($spectrum)=@_; my ($it,$i,$it_start) =0; my ($violtot)=0; my (@violated); # Prints the violations from the calib file. # Note that the actual number of structures used for violation analysis # is defined by 'assignstruct_?_pro=...'! In run.cns. &print_title("Violation analysis for spectrum $spectrum",OUTPUT); if ($violprintall eq 'off') {$it_start=$#it_dir;} else {$it_start = 0;} for ($i=$it_start;$i<=$#it_dir;$i++) { undef @viol; $it=$it_dir[$i]; if ($itname[$it] =~ /Final/) { $ittext="Final Data"; } else { $ittext="Iteration $it"; } $strongviol_num=0; if (defined $strongviol{$spectrum}[$it][$strongviol_num]) { $violtot=$#{$strongviol{$spectrum}[$it]}+1; print OUTPUT "\n*$ittext: Removed $violtot violated restraints for $spectrum.\n"; print OUTPUT "\nList of consistent violations not used in calculation....\n"; while (defined $strongviol{$spectrum}[$it][$strongviol_num]) { # Uncomment to print short list of violated peaks to screen. #print $strongviol{$spectrum}[$it][$strongviol_num].","; printf OUTPUT ("\n *Peak %-5d [%.2f-%.2f] +++ %.2f (in %d/%d)\n", $strongviol{$spectrum}[$it][$strongviol_num], $strongviol_lol{$spectrum}[$it][$strongviol_num], $strongviol_upl{$spectrum}[$it][$strongviol_num], $strongviol_rms{$spectrum}[$it][$strongviol_num], $strongviol_thres{$spectrum}[$it][$strongviol_num], $strongviol_structurestotal[$it]); $peak_num=$strongviol{$spectrum}[$it][$strongviol_num]; push (@{$viol[0]},$peak_num); push (@{$viol[1]},$origres_ppm1{$spectrum}[$peak_num]); push (@{$viol[2]},$origres_ppm2{$spectrum}[$peak_num]); foreach $el ('i','j') { @{$atname{$el}}=split(',',$strongviol{$el}{$spectrum}[$it][$strongviol_num]); if ($violcompress ne 'off') { @{$atname{$el}}=&mergeatomnames(@{$atname{$el}}); } } for ($line=0;($line<=$#{$atname{'i'}} || $line<=$#{$atname{'j'}});$line++) { @atcomp_i=split(' ',$atname{'i'}[$line]); @atcomp_j=split(' ',$atname{'j'}[$line]); printf OUTPUT (" %3s %-3s %-5s %3s %-3s %-5s\n", $atcomp_i[1],$atcomp_i[0],$atcomp_i[2], $atcomp_j[1],$atcomp_j[0],$atcomp_j[2]); } $strongviol_num++; } # Uncomment to print short list of violated peaks to screen. #print "\n"; # Routine to print chemical shift values with peak nums (sorted & short format) for ($x=0;$x<=$#{$viol[0]};$x++) { for ($y=($x+1);$y<=$#{$viol[0]};$y++) { if ($viol[1][$x] < $viol[1][$y]) { for $k (0...2) { ($viol[$k][$x],$viol[$k][$y])=&swap($viol[$k][$x],$viol[$k][$y]) } } } } &print_title ("Short violation list for spectrum $spectrum",REPORT); for ($x=0;$x<=$#{$viol[0]};$x++) { printf REPORT (" *Peak %4d ppm1 %7.2f ppm2 %7.2f\n",$viol[0][$x],$viol[1][$x],$viol[2][$x]); } print REPORT "\n"; } elsif ($calib{$spectrum}[$it]==0) { print OUTPUT "\n$ittext: No data available.\n"; } else {print OUTPUT "\n$ittext: No consistent strong violations.\n";} } } sub reportunambass { my ($spectrum,$it) = @_; my ($rsdnum1,$rsdnum2); my ($unamblist,$ambiglist,$longlist_ok,$longlist_check); # Loop over all peaks (restraints) for $restr (1...$#{$restr{$spectrum}[$it]}) { # only look at unambiguous restraints. next if ($restr{$spectrum}[$it][$restr] != 1); if ($#{$restr_rsdnum{$spectrum}[$it][$restr][0]} = 1) { $rsdnum1=$restr_rsdnum{$spectrum}[$it][$restr][0][0]; $rsdnum2=$restr_rsdnum{$spectrum}[$it][$restr][1][0]; $peakid=$restr_topeak{$spectrum}[$it][$restr]; $restr_orig=$peak_torestr{$spectrum}[$it[0]][$peakid]; ($res1,$res2)=sort {$a <=> $b} ($rsdnum1,$rsdnum2); $ppm1=$origres_ppm1{$spectrum}[$peakid]; $ppm2=$origres_ppm2{$spectrum}[$peakid]; $ppm3=$origres_ppm3{$spectrum}[$peakid]; # Report only if it was originally not assigned if ($restr{$spectrum}[$it[0]][$restr_orig] != 1) { if ($ppm3 != 0) { $unamblist.=sprintf "*Peak $peakid $spectrum (ppm1 $ppm1, ppm2 $ppm2, ppm3 $ppm3):";} else { $unamblist.=sprintf "*Peak $peakid $spectrum (ppm1 $ppm1, ppm2 $ppm2):";}; $unamblist.=sprintf ("%3d %-5s %3s %-5s\n", $restr_rsdnum{$spectrum}[$it][$restr][0][0], $restr_atname{$spectrum}[$it][$restr][0][0], $restr_rsdnum{$spectrum}[$it][$restr][1][0], $restr_atname{$spectrum}[$it][$restr][1][0]); } } } &print_title("New unambiguous assignments",NEWUNAMBIG); print NEWUNAMBIG $unamblist; } sub reportambass { my ($spectrum,$it) = @_; my ($rsdnum1,$rsdnum2); my ($unamblist,$amblist,$longlist_ok,$longlist_check); # Loop over all peaks (restraints) for $restr (1...$#{$restr{$spectrum}[$it]}) { # only look at ambiguous restraints. if ($#{$restr_rsdnum{$spectrum}[$it][$restr][0]} > 1) { $peakid=$restr_topeak{$spectrum}[$it][$restr]; $restr_orig=$peak_torestr{$spectrum}[$it[0]][$peakid]; ($res1,$res2)=sort {$a <=> $b} ($rsdnum1,$rsdnum2); $ppm1=$origres_ppm1{$spectrum}[$peakid]; $ppm2=$origres_ppm2{$spectrum}[$peakid]; $ppm3=$origres_ppm3{$spectrum}[$peakid]; # Report only if it was originally not assigned if ($restr{$spectrum}[$it[0]][$restr_orig] != 1) { $stop=0; for $k (0...$#{$restr_rsdnum{$spectrum}[$it][$restr][0]}) { if ($ppm3 != 0) { $amblist.=sprintf "*Peak $peakid $spectrum (ppm1 $ppm1, ppm2 $ppm2, ppm3 $ppm3):";} else { $amblist.=sprintf "*Peak $peakid $spectrum (ppm1 $ppm1, ppm2 $ppm2):";}; $amblist.=sprintf ("%3d %-5s %3s %-5s\n", $restr_rsdnum{$spectrum}[$it][$restr][0][$k], $restr_atname{$spectrum}[$it][$restr][0][$k], $restr_rsdnum{$spectrum}[$it][$restr][1][$k], $restr_atname{$spectrum}[$it][$restr][1][$k]); } $amblist.=sprintf ("\n"); } } } &print_title("New ambiguous assignments",NEWAMBIG); print NEWAMBIG $amblist; } sub reportbadass { my ($spectrum,$it) = @_; my ($rsdnum1,$rsdnum2); my ($unamblist,$longlist_ok,$longlist_check); # Loop over all peaks (restraints) # In this loop unambiguous long range contacts that had # an intraresidue or sequential ambiguous possibility are # detected for $restr (1...$#{$restr{$spectrum}[$it]}) { # only look at unambiguous restraints. next if ($restr{$spectrum}[$it][$restr] != 1); $rsdnum1=$restr_rsdnum{$spectrum}[$it][$restr][0][0]; $rsdnum2=$restr_rsdnum{$spectrum}[$it][$restr][1][0]; # Detect an unambiguous medium or long range NOE... if ( (($rsdnum1-$rsdnum2)**2) > 1) { $peakid=$restr_topeak{$spectrum}[$it][$restr]; $restr_orig=$peak_torestr{$spectrum}[$it[0]][$peakid]; ($res1,$res2)=sort {$a <=> $b} ($rsdnum1,$rsdnum2); # If it was originally ambiguously assigned, check if there # was a sequential or intraresidue possibility. if ($restr{$spectrum}[$it[0]][$restr_orig] != 1) { $stop=0; for $k (0...$#{$restr_rsdnum{$spectrum}[$it[0]][$restr_orig][0]}) { $rsdnumt1=$restr_rsdnum{$spectrum}[$it[0]][$restr_orig][0][$k]; $rsdnumt2=$restr_rsdnum{$spectrum}[$it[0]][$restr_orig][1][$k]; if ((($rsdnumt1-$rsdnumt2)**2) <= 1) { $ppm1=$origres_ppm1{$spectrum}[$peakid]; $ppm2=$origres_ppm2{$spectrum}[$peakid]; $unamblist.=sprintf "\n*Peak $peakid in spectrum $spectrum (ppm1 $ppm1, ppm2 $ppm2)\n"; $unamblist.=sprintf (" Final unambiguous assignment: %3d %-5s %3s %-5s\n", $restr_rsdnum{$spectrum}[$it][$restr][0][0], $restr_atname{$spectrum}[$it][$restr][0][0], $restr_rsdnum{$spectrum}[$it][$restr][1][0], $restr_atname{$spectrum}[$it][$restr][1][0]); $unamblist.=sprintf (" Original ambiguous possibility: %3d %-5s %3s %-5s\n", $rsdnumt1, $restr_atname{$spectrum}[$it[0]][$restr_orig][0][$k], $rsdnumt2, $restr_atname{$spectrum}[$it[0]][$restr_orig][1][$k]); $stop=1;last; } last if ($stop==1); } } } } # Here detect restraints that only have assignment possibilities corresponding to # the more unlikely areas of the contact map for $restr (1...$#{$restr{$spectrum}[$it]}) { $contact_maphit=0;$reslist=""; for $k (0...$#{$restr_rsdnum{$spectrum}[$it][$restr][0]}) { $rsdnum1=$restr_rsdnum{$spectrum}[$it][$restr][0][$k]; $rsdnum2=$restr_rsdnum{$spectrum}[$it][$restr][1][$k]; ($res1,$res2) = sort {$a <=> $b} ($rsdnum1,$rsdnum2); # Sometimes additional line for unambiguous ones - quick fix # for that here next if ($res1 eq '' || $res2 eq ''); # Keep track of residue combinations that are one. If so, # skip this restraint if ($contact_dep[$it][$res1][$res2] > $contact_avg[$it]) { $contact_maphit=1; last; } else { if ($reslist !~ /\[$rsdnum1\,$rsdnum2\]\,/) { $reslist.="[$rsdnum1,$rsdnum2],"; } } } if ($contact_maphit!=1) { $peakid=$restr_topeak{$spectrum}[$it][$restr]; $ppm1=$origres_ppm1{$spectrum}[$peakid]; $ppm2=$origres_ppm2{$spectrum}[$peakid]; chop $reslist; $longlist_check .= sprintf (" * Peak %5d %s (ppm1 %.3f, ppm2 %.3f).\n",$peakid,$reslist,$ppm1,$ppm2); } } &print_title("Unambiguous long range with possible intra or sequential assignment",REPORT); print REPORT $unamblist; &print_title("Long range contacts - check",REPORT); print REPORT "\n".$longlist_check; } sub reportdetails { my ($spectrum,$it,$printing) = @_; my ($rsdnum1,$rsdnum2,$double_it); my (@longinfo,%info_count,$count,$res1,$res2)=0; my (%contact,%contact_info,%info); my (@contact_res1,@contact_res2,$detaillist,@contact,$type)=""; my (@type_list)=('long','med','seq','intra'); undef @contact; $distvolmax = $restr_distvolmax{$spectrum}[$it]; #print $distvolmax."\n"; # Loop over all peaks (restraints) for $restr (1...$#{$restr{$spectrum}[$it]}) { undef @longinfo; undef @contact_res1; undef @contact_res2; undef %contact; undef %info_count; $total=0; $info_count{'long'}=0; $addtext=""; if ($it_dir[0] == $it_final) {$double_it=2;} else {$double_it=1;} # Loop over all possibilities for spectrum.tbl list! for $k (0...$#{$restr_rsdnum{$spectrum}[$it][$restr][0]}) { $rsdnum1=$restr_rsdnum{$spectrum}[$it][$restr][0][$k]; $rsdnum2=$restr_rsdnum{$spectrum}[$it][$restr][1][$k]; $atname2=$restr_atname{$spectrum}[$it][$restr][1][$k]; # Make a multiplication factor for this particular contact (see below where $contact[$it][][] is assigned) $atname2=~ s/\%//; chop $atname2; ($res1,$res2) = sort {$a <=> $b} ($rsdnum1,$rsdnum2); push (@contact_res1,$res1); push (@contact_res2,$res2); if ($res2 > $num_residues) {$num_residues=$res2;} # Long range combinations (more than 3 residues apart) if ((($rsdnum1-$rsdnum2)**2) > 9) { $type='long'; $longinfo[$info_count{$type}]=$k; } # Sequential elsif ((($rsdnum1-$rsdnum2)**2) == 1) { $type='seq'; } # Intraresidue elsif (($rsdnum1-$rsdnum2) == 0) { $type='intra'; } # Medium range else { $type='med'; } $info_count{$type}++; $contact{$type}[$rsdnum1]++; $contact{$type}[$rsdnum2]++; } # Calculate results for the matrix. Weighs the contacts on the number of total # assignment possibilities for a restraint! $total=$#{$restr_rsdnum{$spectrum}[$it][$restr][0]}+1; foreach $type (@type_list) { $info{$type}{$spectrum}[$it]+=$info_count{$type}/$total; } # Only run if -sp flag is set to on - default is off. This is an addition # on 13/06/2001 - now uses unambig/ambig.tbl as default to calculate # contact matrix if ($specmatrix eq 'on') { $int_f=($restr_distvol{$spectrum}[$it][$restr])*5/$distvolmax; for $l (0...$#contact_res1) { $res1=$contact_res1[$l]; $res2=$contact_res2[$l]; $contact[$it][$res1][$res2]+=($int_f/$total); } } # Count types of restraint per residue if necessary (only for last iteration read) if ($printrestperres eq 'on' && $it==$it_final) { foreach $type (@type_list) { # Is OK - only have to go to maximal res num found for $res (1...$num_residues) { $contact_info{$type}[$res]+=$contact{$type}[$res]/$total; } } } # Don't look at unambiguous restraints and non-long range in detailed overview. next if ($restr{$spectrum}[$it][$restr] == 1); next if ($info_count{'long'} == 0); # Print out info if long range detected. $peakid=$restr_topeak{$spectrum}[$it][$restr]; if ($info_count{'long'} != $total) { $addtext=sprintf("(other possibilities: %d intrares, %d sequential, %d medium range)", $info_count{'intra'},$info_count{'seq'},$info_count{'med'}); } $detaillist.=sprintf ("*Peak %d $addtext\n",$peakid); for $count (0...($info_count{'long'}-1)) { $rsdnum1=$restr_rsdnum{$spectrum}[$it][$restr][0][$longinfo[$count]]; $atname1=$restr_atname{$spectrum}[$it][$restr][0][$longinfo[$count]]; $rsdnum2=$restr_rsdnum{$spectrum}[$it][$restr][1][$longinfo[$count]]; $atname2=$restr_atname{$spectrum}[$it][$restr][1][$longinfo[$count]]; ($res1,$res2) = sort {$a <=> $b} ($rsdnum1,$rsdnum2); $detaillist.= sprintf (" * Atom %3d %-4s to %3d %-4s (%3d\-%-3d)\n",$rsdnum1,$atname1,$rsdnum2,$atname2,$res1,$res2); } $detaillist .= sprintf "\n"; } printf REPORT "\nDistribution contacts iteration $it, spectrum $spectrum:\n"; printf REPORT "\n Intraresidue: %6.1f",$info{'intra'}{$spectrum}[$it]; printf REPORT "\n Sequential: %6.1f",$info{'seq'}{$spectrum}[$it]; printf REPORT "\n Medium range: %6.1f",$info{'med'}{$spectrum}[$it]; printf REPORT "\n Long range: %6.1f\n",$info{'long'}{$spectrum}[$it]; # Put matrix data in global matrix variables (if -sp flag set to on) if ($specmatrix eq 'on') { $k_max = $#{$contact[$it]}; for $k (1...$k_max) { for $l (1...$#{$contact[$it][$k_max]}) { next if ($k>$l); $contact_all[$it][$k][$l]+=$contact[$it][$k][$l]; # for the full overview (all spectra) $contact_all_sum[$it]+=$contact[$it][$k][$l]; # to get an 'average' $contact_all_max[$it]=$contact_all[$it][$k][$l] if ($contact_all_max[$it]<$contact_all[$it][$k][$l]); } } } if ($printrestperres eq 'on' && $it==$it_final) { print RESTPERRES "Spectrum $spectrum\n\n "; foreach $type (@type_list) { printf RESTPERRES ("%8s",$type); } for $res (1...$num_residues) { $line = sprintf ("\n%3d",$res); foreach $type (@type_list) { $line.=sprintf(" %5.1f ",$contact_info{$type}[$res]); #Keep track of total over all spectra $contact_info_all{$type}[$res]+=$contact_info{$type}[$res]; } print RESTPERRES $line; } print RESTPERRES"\n\n"; } if ($printing eq 'on') { # Print the list of peaks &print_title("Spectrum $spectrum: Long range possibilities in ambiguous restraints.",DETAIL); print DETAIL "\n"; print DETAIL $detaillist; } } #################################### # Subroutines for reading files... # #################################### sub getpdbinfo { my ($it) = @_; my (@pdblist,$gone,@eflds) = ""; my ($pdbnum,$pdbcount,$pdbfile,$max_value) = 0; my (%E); opendir (DIR, "structures/it$it"); @pdblist = grep /\d+.pdb$/, readdir DIR ; closedir DIR; $pdball[$it]=$#pdblist+1; if ($#pdblist < 0) { &send_warning("No .pdb files found for iteration $it.\n",STDOUT); return 0; } $status=open(PDBFILELIST,"structures/it$it/file.nam"); if ($status==0 && $#pdblist >= $pdbbest) { # Get overall energy info from pdb files if no file.nam print "\nGetting info from PDB files!\n"; foreach $pdbfile (@pdblist) { $status=open (PDB,"structures/it$it/$pdbfile"); if ($status == 0) {next;} while (defined ($_ = <PDB>)) { if (/^REMARK\s+overall/) { $_=<PDB>; s/^REMARK energies://; s/\s//g; @eflds=split(','); $E{$pdbfile}=$eflds[0]; $Ecount++; if ($Ecount>$pdbbest) { $max_value=$E{$pdbfile}; $gone=$pdbfile; foreach $pdbkey (keys %E) { if ($E{$pdbkey} > $max_value) { $gone=$pdbkey; $max_value=$E{$pdbkey}; } } delete $E{$gone}; } last; } } } undef @pdblist; @pdblist = keys %E; } elsif ($status==1 && $#pdblist >= $pdbbest) { undef @pdblist; while (defined ($_=<PDBFILELIST>)) { next if (!/.pdb/); chomp; $pdbcount++; push (@pdblist,$_) if ($pdbcount <= $pdbbest); } } close PDBFILELIST; $root_name = $pdblist[0]; $root_name =~ s/_\d+\.pdb$//g; print ("\nReading .pdb files...\n"); foreach $pdbfile (@pdblist) { $status=open (PDB,"structures/it$it/$pdbfile"); if ($status == 0) {&sendwarning("Could not open pdb file: $!.\n",STDOUT);next;} $pdbcount[$it]++; undef @pdb_infolist; $pdbfile =~ s/^$root_name\_//; $pdbfile =~ s/.pdb$//; while (defined ($_ = <PDB>)) { if (/^REMARK\s+[overall|bonds|noe]/) { $_=<PDB>; s/^REMARK\s+.+://; s/\s//g; push(@pdb_infolist,split(',')); #pop(@pdb_infolist); next; } elsif (/^ATOM/) {last;} } for $i (0...$#PDBinfolist) { $pdb{$PDBinfolist[$i]}[$it][$pdbnum]=$pdb_infolist[$i]; } #uncomment next for list of energies for pdb files #printf ("$pdbnum %8.2f %8.2f\n",$pdb_infolist[0],$pdb_infolist[5]); $pdbnum++; close PDB; } return 1; } sub getcalibinfo { # Info here: the list of 'problem' restraints is the same for # each calibration - can basically skip this once you've done # it for an iteration... if not just display the info once! # So reads in data/.../....tbl EACH time anew!!! my ($it) = @_; my (@caliblist,@caliblist_sort) = ""; my ($restr_na,$error_before,$noe_end) = 0; if ($itname[$it] =~ /Final/) { $final=($it-1).'/finaldata'; $it_text="it".($it-1)."/finaldata/"; $itr=$it-1; } else {$final=$it;$it_text="it$it/";$itr=$it;} print ("\n"); foreach $calib (@spec_list) { $status=open (CALIB,"structures/it$final/calib_$calib\_it$itr.out"); if ($status == 0) {&send_warning("No $calib\_it$itr.out file for iteration $it_text: $!.\n",STDOUT);return;} $strongviol_num=0; $restr_na=0; print ("Reading it$final/calib_$calib\_it$itr.out...\n"); if ($restr_prob{$calib} != 1) { &print_title("Spectrum $calib: errors while reading original restraint list",OUTPUT); } while (defined ($_ = <CALIB>)) { if (/reading restraint/ && $restr_prob{$calib} != 1) { s/\s*reading restraint\s+//g; chomp; $restraint=$_; $_=<CALIB>; if (!/restraint successfully/) { $restr_na{$calib}[$restr_na]=$restraint; $restr_na++; $peakid=$origlistid{$calib}[$restraint]; $ppm_list=sprintf ("ppm1 %6.2f ppm2 %6.2f",$origres_ppm1{$calib}[$peakid],$origres_ppm2{$calib}[$peakid]); if ($origres_ppm3{$calib}[$peakid] != 0) {$ppm_list.=sprintf (" ppm3 %6.2f",$origres_ppm3{$calib}[$peakid]);} if (/NOESET\-ERR/) { printf OUTPUT (" *Restraint %-5d (No assignment! Peak %-5d $ppm_list)\n",$restraint,$peakid); } elsif (/SELRPN\-ERR/) { printf OUTPUT (" *Restraint %-5d (Selection error! Peak %-5d $ppm_list)\n",$restraint,$peakid); } else { printf OUTPUT (" *Restraint %-5d (Unknown error! Peak %-5d $ppm_list)\n",$restraint,$peakid); } } next; } elsif (/NOE\-ERR\: unrecognized command\:/ && $restraint != $error_before && $restr_prob{$calib} != 1) { $error_before=$restraint; $peakid=$origlistid{$calib}[$restraint]; printf OUTPUT (" *Error reading after restraint %-5d (Peak %-5d)\n",$restraint,$peakid); next; } elsif (/========== restraint/) { /\d+/; $strongviol{$calib}[$it][$strongviol_num]=$&; $continue=1;$type=''; while ($continue==1) { $_=<CALIB>; s/^\s+//;s/\s+$//;s/\s+/ /g; # strip and unify spaces if (/set-i-atoms/) { $type='i'; } elsif (/set-j-atoms/) { $type='j'; } elsif (/Violations Above/) { s/Violations?|Above|Threshold//g; s/Rms|\,|\=|Target|distance|\(|\)|\-|\+|\///g; chop $strongviol{'i'}{$calib}[$it][$strongviol_num]; chop $strongviol{'j'}{$calib}[$it][$strongviol_num]; @el=split(' '); $strongviol_thres{$calib}[$it][$strongviol_num]=$el[0]; $strongviol_rms{$calib}[$it][$strongviol_num]=$el[1]; $strongviol_lol{$calib}[$it][$strongviol_num]=$el[2]-$el[3]; $strongviol_upl{$calib}[$it][$strongviol_num]=$el[2]+$el[4]; $strongviol_num++; $continue=0; } elsif (/\s*\d+\s+/) { $strongviol{$type}{$calib}[$it][$strongviol_num].=$_."," if ($type =~ /i|j/); } } } elsif (/CHPR: d\^-6, int/) { s/CHPR: d\^-6, int\s+//; @fld=split(' '); $calibinfo{$calib}[$it][0]=$fld[0]; $calibinfo{$calib}[$it][1]=$fld[2]; next; } } $restr_prob{$calib} = 1; close CALIB; } } sub getspectblinfo { my ($it) = @_; my ($restr,$restr_poss,$line)=0; if ($itname[$it] =~ /Final/) { $final=($it-1).'/finaldata'; $it_text="it".($it-1)."/finaldata/"; } else {$final=$it;$it_text="it$it/";} print ("\n"); foreach $spectrum (@spec_list) { $status=open (SPECTBL,"structures/it$final/$spectrum.tbl"); if ($status == 0) {&send_warning("No $spectrum.tbl file for iteration $it: $!.\n",STDOUT);return;} print ("Reading $it_text$spectrum.tbl...\n"); $unambiguous=0;$line=0;$restr=0;$restr_poss=0; while (defined ($_ = <SPECTBL>)) { chomp; next if (/^\s*\!/); next if (/^\s*$/); if (/^ ASSI/) { $restr++; $ambig_restr_poss{$spectrum}[$it]+=$restr_poss; $restr_poss=0; $line=0; $unambiguous=1; $unambig_restr{$spectrum}[$it]++; } elsif (/^ OR/) { $line=0; $restr_poss++; if ($restr_poss == 1) { $ambig_restr_poss{$spectrum}[$it]++; $unambig_restr{$spectrum}[$it]--; $ambig_restr{$spectrum}[$it]++; } } elsif (/^\s*\(/) { # Modification by Kevin Gardner 26/10/2001. ($rsdnum,$atname) = /resid\s+(\d+).+name\s+([\w\%\#]+)[\s\)]/; $restr_rsdnum{$spectrum}[$it][$restr][$line][$restr_poss]=$rsdnum; $restr_atname{$spectrum}[$it][$restr][$line][$restr_poss]=$atname; $atname[$line]=$atname; $atname[$line]=~ s/\%//; chop $atname[$line]; $rsdnum[$line]=$rsdnum; # Each OR 'combination' is stored, not a compressed form where # information is lost about the possible combinations. # In this loop the unambiguous or ambiguous status of the # restraint is checked IF $line==1 and for BOTH lines! if ($line == 1 && $unambiguous==1) { for $j (0...($restr_poss-1)) { if ($restr_rsdnum{$spectrum}[$it][$restr][0][$j] != $rsdnum[0] || $restr_rsdnum{$spectrum}[$it][$restr][1][$j] != $rsdnum[1]) { $unambiguous=0;last; } elsif ($restr_atname{$spectrum}[$it][$restr][0][$j] !~ /$atname[0]/ || $restr_atname{$spectrum}[$it][$restr][1][$j] !~ /$atname[1]/) { $unambiguous=0;last; } } $restr{$spectrum}[$it][$restr]=$unambiguous; } $line++; } elsif (/ peak /) { # This gets volume of peak for calculation contact matrix. Other option is # to use upper distance limit. ($v0,$v1,$v2,$v3,$peaknum,$v5,$v6,$v7,$vol,$v9,$ppm1,$v11,$ppm2)=split(' '); $peak_torestr{$spectrum}[$it][$peaknum]=$restr; $restr_topeak{$spectrum}[$it][$restr]=$peaknum; if ($vol >=0) {$restr_distvol{$spectrum}[$it][$restr]=($vol**0.16667);} else {$restr_distvol{$spectrum}[$it][$restr]=0}; #$restr_intvol{$spectrum}[$it][$restr]=($int**0.16667); $restr_distvolmax{$spectrum}[$it]=$restr_distvol{$spectrum}[$it][$restr] if ($restr_distvol{$spectrum}[$it][$restr] > $restr_distvolmax{$spectrum}[$it]); } } $restr_tot{$spectrum}[$it]=$unambig_restr{$spectrum}[$it]+$ambig_restr{$spectrum}[$it]; close SPECTBL; } } sub getrestrinfo { # BADLY WRITTEN ROUTINE! Lotsa duplicate code - should be smoother. my ($it) = @_; my ($res)=""; my ($dis)=0; my (@type_list) = ('long','med','seq','intra'); my ($tempcnt)=-1; my (@contact_res1,@contact_res2)=0; my ($contact_weight)=0; if ($itname[$it] =~ /Final/) { $final=($it-1).'/finaldata'; $it_text="it".($it-1)."/finaldata/"; } else {$final=$it;$it_text="it$it/";} $restr[$it]=open (AMBIG,"structures/it$final/ambig.tbl"); if ($restr[$it] == 0) {&send_warning("No ambig.tbl file for iteration $it_text: $!.\n",STDOUT);return;} print ("\nReading ".$it_text."ambig.tbl...\n"); while (defined ($_ = <AMBIG>)) { if (/^ ASSI/) { # This misses out on last restraint - added quick patch bottom foreach $type (@type_list) { $restrcnt{$type}[$it]+=$restrcnttemp{$type}/$tempcnt; $restrcnt{$type."amb"}[$it]+=$restrcnttemp{$type}/$tempcnt; $restrcnttemp{$type}=0; } # Add to restraint contact matrix, reset, if -sp to off if ($specmatrix ne 'on') { $contact_weight/=$tempcnt; for $l (0...$#contact_res1) { $res1=$contact_res1[$l]; $res2=$contact_res2[$l]; $contact_all[$it][$res1][$res2]+=$contact_weight; $contact_all_sum[$it]+=$contact_weight; # to get an 'average' $contact_all_max[$it]=$contact_all[$it][$res1][$res2] if ($contact_all_max[$it]<$contact_all[$it][$res1][$res2]); } } undef @contact_res1; undef @contact_res2; # Keep track of numbers... $ambig_restr_poss[$it]+=$tempcnt if ($tempcnt != -1); $ambig_restr[$it]++; $tempcnt=1; $res=""; } elsif (/^ OR/) { $tempcnt++; $res=""; } elsif (/peak/) { ($med,$min,$plus,$v3,$v4,$v5,$wght,$v7,$vol,@void)=split(' '); # By default, use ($med+($plus-$min)/2)^-6 value times $wght*1000 for 'weight' restraint. # Can alternatively use upper limit, volume, ... $contact_weight=(($med+($plus-$min)/2)**-6)*$wght*1000; } elsif (/segid/){ s/.+and resid\s+|\s+and name.+\n//g; # Second line if ($res ne "") { # Count long, intra, ... $restrcnttemp{&checkcontact($res,$_)}++; # Keeps track of residues for contact map ($res1,$res2) = sort {$a <=> $b} ($res,$_); push (@contact_res1,$res1); push (@contact_res2,$res2); } # First line else { $res=$_; } } } # Quick patch for counting last restraint. foreach $type (@type_list) { $restrcnt{$type}[$it]+=$restrcnttemp{$type}/$tempcnt; $restrcnt{$type."amb"}[$it]+=$restrcnttemp{$type}/$tempcnt; $restrcnttemp{$type}=0; } # Add to restraint contact matrix, also quick patch. if ($specmatrix ne 'on') { $contact_weight/=$tempcnt; for $l (0...$#contact_res1) { $res1=$contact_res1[$l]; $res2=$contact_res2[$l]; $contact_all[$it][$res1][$res2]+=$contact_weight; $contact_all_sum[$it]+=$contact_weight; # to get an 'average' $contact_all_max[$it]=$contact_all[$it][$res1][$res2] if ($contact_all_max[$it]<$contact_all[$it][$res1][$res2]); } } close AMBIG; $restr[$it]=open (UNAMBIG,"structures/it$final/unambig.tbl"); if ($restr[$it] == 0) {&send_warning("No unambig.tbl file for iteration $it_text: $!.\n",STDOUT);return;} print ("Reading ".$it_text."unambig.tbl...\n"); while (defined ($_ = <UNAMBIG>)) { if (/^ ASSI/) { $unambig_restr[$it]++; $res=""; } elsif (/^ OR/) { $res="NO"; } elsif (/segid/ && $res ne "NO"){ s/.+and resid\s+|\s+and name.+\n//g; if ($res ne "") { $restrcnttemp{&checkcontact($res,$_)}++; ($res1,$res2) = sort {$a <=> $b} ($res,$_); } $res=$_; } elsif (/peak/) { foreach $type (@type_list) { $restrcnt{$type}[$it]+=$restrcnttemp{$type}; $restrcnt{$type."unamb"}[$it]+=$restrcnttemp{$type}; $restrcnttemp{$type}=0; } if ($specmatrix ne 'on') { ($med,$min,$plus,$v3,$v4,$v5,$wght,$v7,$vol,@void)=split(' '); # By default, use ($med+($plus-$min)/2)^-6 value times $wght*1000 for 'weight' restraint. # Can alternatively use upper limit, volume, ... $contact_weight=(($med+($plus-$min)/2)**-6)*$wght*1000; $contact_all[$it][$res1][$res2]+=$contact_weight; $contact_all_sum[$it]+=$contact_weight; # to get an 'average' $contact_all_max[$it]=$contact_all[$it][$res1][$res2] if ($contact_all_max[$it]<$contact_all[$it][$res1][$res2]); } } } close UNAMBIG; $restr_tot[$it]=$unambig_restr[$it]+$ambig_restr[$it]; } sub checkcontact { my ($res1,$res2)=@_; $dis=sqrt (($res1-$res2)**2); if ($dis==0) { return 'intra'; } elsif ($dis==1) { return 'seq'; } elsif ($dis<=9) { return 'med'; } else { return 'long'; } } #################################### # Other scriptspecific subroutines # #################################### sub readitfiles { my ($it)=@_; &print_title ("Reading \l$itname[$it]"); &getcalibinfo($it) if ($calibprint ne 'off' || $violprintall ne 'off' || $violprintlast ne 'off'); &getspectblinfo($it) if ($restrover_spec ne 'off'); &getrestrinfo($it) if ($restrover ne 'off'); } sub makedependentmatrix { my ($it,$all,$all_max,@pointer) = @_; my ($k_max)=$#{$pointer[$it]}; my ($l_max)=$#{$pointer[$it]->[$k_max]}; my ($dep_max,$avg,$factor)=0; $avg=$all*2/($num_residues**2); $factor=($all_max/$avg)*10/$num_residues; $limit_min=$avg/$factor; for $k (1...$num_residues) { for $l ($k...$num_residues) { $poss=9; $poss-=3 if ($k==$l); $poss-=3 if ($k==1); $poss-=3 if ($l==$num_residues); $poss-=1 if ($l-$k==1); $poss+=1 if ($l==$num_residues && $k==1); for $r (($k-1)...($k+1)) { for $s (($l-1)...($l+1)) { $contact_dep[$it][$k][$l]+=$pointer[$it]->[$r]->[$s]; } } $contact_dep[$it][$k][$l]=($contact_dep[$it][$k][$l])/$poss; $dep_max = $contact_dep[$it][$k][$l] if ($dep_max < $contact_dep[$it][$k][$l]); # Determine 'filtered' map on lower half display. # # 0 if dep. value smaller than $limit_min (depends on max & min value & num residues). # 0 if dep. value smaller than $limit_min*2 AND actual value smaller than $limit_min*2 # 1 in all other cases. if ($contact_dep[$it][$k][$l] >= $limit_min*2) { $contact_map[$it][$k][$l]=1; next; } else { if ($contact_dep[$it][$k][$l] < $limit_min) { $contact_map[$it][$k][$l]=0; } elsif ($contact_dep[$it][$k][$l] > $limit_min && $pointer[$it]->[$r]->[$s] < $limit_min*2) { $contact_map[$it][$k][$l]=0; } else {$contact_map[$it][$k][$l]=1;} } } } $contact_dep_max[$it] = $dep_max; return $avg; } sub print_help { print "\n ".$scriptname_short." recognizes the following options:\n"; print $helptext; } sub mergeatomnames { my (@atomlist) = @_; my ($metlist,$metsep,@atfound,$atom,$atfoundlist,$atfoundsep)=""; my ($atcount)=0; my $sep=""; while (defined ($atomlist[$atcount])) { $atom = $atomlist[$atcount]; if ($atom =~ /$methyls/) { # For methyls if ($atom !~ /$metlist/ || $metlist eq "") { if ($metlist eq "") {$sep="";} else {$sep="|";} $metlist.=$sep.$atom; chop $metlist; $atomlist[$atcount] =~ s/\d$/\#/; } else { splice(@atomlist,$atcount,1); $atcount--; } } else { # For recurring same name @atfound=grep (/$atom/,@atomlist); if ($#atfound > 0) { splice(@atomlist,$atcount,1); $atcount--; } # For methylenes elsif ($atom =~ /\d$/) { chop $atom; @atfound=grep (/$atom/,@atomlist); if ($atom =~ /TRP/ && $atom !~ /HB/) {;} elsif ($#atfound == 1) { splice(@atomlist,$atcount,1); $atcount--; if ($atfoundlist eq "") {$sep="";} else {$sep="|";} $atfoundlist.=$sep.$atom; } elsif ($atom =~ /$atfoundlist/ && $atfoundlist ne "") { $atomlist[$atcount] =~ s/\d$/\#/; } } } $atcount++; if ($atcount < 0) { print "Bad ".$peak_num." $el $atcount\n"; } } return @atomlist; } ####################### # General use subs... # ####################### sub swap { my ($swap1,$swap2)=@_; return $swap2,$swap1; } sub print_title { my ($title,$output) = @_; if (!defined $output) {$output=STDOUT}; $title_cap=$title="# ".$title." #"; $title_cap =~ s/./#/g; print $output "\n".$title_cap."\n".$title."\n".$title_cap."\n"; } sub send_warning { my ($warning,$output) = @_; print $output $warning; return; } sub send_fatal { my ($warning) = @_; print "Program exiting!\n"; die $warning; return; } sub minavgmax { my ($arraypointer) = @_; my ($min,$avg,$max,$num,$tot); return if ($#{$arraypointer} < 0); $min = $arraypointer->[0]; for $num ( 0 .. $#{$arraypointer}) { if ($max < $arraypointer->[$num]) { $max = $arraypointer->[$num]; } if ($min > $arraypointer->[$num]) { $min = $arraypointer->[$num]; } $tot += $arraypointer->[$num]; } $avg = $tot / ($#{$arraypointer} + 1); return ($min,$avg,$max); } sub sd { my ($avg,$arraypointer) = @_; my ($sd,$num,$sse); for $num ( 0 .. $#{$arraypointer}) { $sse += (($arraypointer->[$num])-$avg)**2; } if ($#{$arraypointer} !=0) { $sd = sqrt ($sse/$#{$arraypointer}); } else { $sd = sqrt ($sse/($#{$arraypointer}+1)); } return ($sd); } sub gettime { my ($d,$dt); # Get time and date $dt = `date '+%d/%m/%Y at %H:%M:%S'`; $d = `date '+%d.%m.%Y'`; chomp $d;chomp $dt; return $dt,$d; } sub iduser { # ID person calling script... $_ = `who am i`; $_ = substr($_,0,index($_," ")); s/^.+\!//; return $_; } sub getprojectinfo { my ($cur_dir,@dir_comp); # Get project name info... $cur_dir=`pwd`; chomp $cur_dir; @dir_comp=split(/\//,$cur_dir); return $dir_comp[($#dir_comp-1)],$dir_comp[$#dir_comp]; } sub openoutputfile { my ($file,$handle,$text)=@_; # Open output file with header line open ($handle,">$file") || die "Unable to open file $file: $!"; print $handle $text; } ####################################### # Subroutines for making xmgr file... # ####################################### sub makegraph { my ($type) =@_; my (@coord,@graph_text,$ydiv,$yadd,$spectrum,$tot_spec,$min,$avg,$max) = 0; my ($name) =""; open(GRAPHFILE,">$outputfile.$type.xmgr")|| die "Unable to make xmgr file: $!"; &print_graph_initialize(GRAPHFILE); if ($type eq 'restraints') { # Print overall stats. Is fixed. # Info in lists comes from print_restr subroutines. chomp $unambig_rest_list; chomp $ambig_rest_list; chomp $all_rest_list; @coord=(0.11,0.44,0.65,0.88); @graph_text=("Restraint progression",'Iterations (overall)','Number of restraints'); &print_graph('xydydy',GRAPHFILE,0,@coord,@graph_text,'All',$all_rest_list,4,'Unambiguous',$unambig_rest_list,3,'Ambiguous',$ambig_rest_list,2); @coord=(0.59,0.92,0.65,0.88); @graph_text=("Relative progression",'Iterations (overall)','Ratio'); &print_graph('xydydy',GRAPHFILE,1,@coord,@graph_text,'Poss amb',$avg_amb_list,2,'Amb/unamb',$ratio_list,100); # Print graph stats $tot_spec = $#spec_list+1; if ($tot_spec > 1 ) {$ydiv=(0.6/$tot_spec);$yadd=0.05;$ysub=0.07;} else {$ydiv=0.28;$yadd=0.3;$ysub=0.04;} for $i (0...$#spec_list) { $spectrum=$spec_list[$i]; @coord=(0.11,0.44,($i*$ydiv)+$yadd,($i*$ydiv)+$ydiv+$yadd-$ysub); @graph_text=("","Iterations (detail $spectrum)",'Number of restraints'); &print_graph('xydydy',GRAPHFILE,2+($i*2),@coord,@graph_text,'All',$all_rest_list{$spectrum},4,'Unambiguous',$unambig_rest_list{$spectrum},3,'Ambiguous',$ambig_rest_list{$spectrum},2); @coord=(0.59,0.92,($i*$ydiv)+$yadd,($i*$ydiv)+$ydiv+$yadd-$ysub); @graph_text=("","Iterations (detail $spectrum)",'Ratio'); &print_graph('xydydy',GRAPHFILE,3+($i*2),@coord,@graph_text,'Poss amb',$avg_amb_list{$spectrum},2,'Amb/unamb',$ratio_list{$spectrum},100); } } if ($type eq 'pdb') { # No new pdb files for final data - don't print that. if ($it_dir[$#it_dir]+1==$it_final) { $it_last=$it_final-1; } else { $it_last=$it_final; } if ((($it_last)-3)<0) {$detail_start=0;} else {$detail_start=$it_last-3;} @rest=&makegraphpdbdata('xydydy',$it_dir[0],$it_last,@Egraphlist); @coord=(0.11,0.44,0.65,0.88); @graph_text=('Energy progression','','Energy'); &print_graph('xydydy',GRAPHFILE,0,@coord,@graph_text,@rest); @rest=&makegraphpdbdata('xydydy',$detail_start,$it_last,@Egraphlist); @coord=(0.59,0.92,0.65,0.88); @graph_text=('Energy progression detail','','Energy'); &print_graph('xydydy',GRAPHFILE,1,@coord,@graph_text,@rest); @rest=&makegraphpdbdata('xydydy',$it_dir[0],$it_last,@rmsgraphlist); @coord=(0.11,0.44,0.37,0.55); @graph_text=('Rms progression','','rms'); &print_graph('xydydy',GRAPHFILE,2,@coord,@graph_text,@rest); @rest=&makegraphpdbdata('xydydy',$detail_start,$it_last,@rmsgraphlist); @coord=(0.59,0.92,0.37,0.55); @graph_text=('Rms progression detail','','rms'); &print_graph('xydydy',GRAPHFILE,3,@coord,@graph_text,@rest); @rest=&makegraphpdbdata('xydydy',$it_dir[0],$it_last,@violgraphlist); @coord=(0.11,0.44,0.07,0.25); @graph_text=('Violations progression','Iterations','Number'); &print_graph('xydydy',GRAPHFILE,4,@coord,@graph_text,@rest); @rest=&makegraphpdbdata('xydydy',$detail_start,$it_last,@violgraphlist); @coord=(0.59,0.92,0.07,0.25); @graph_text=('Violations progression detail','Iterations','Number'); &print_graph('xydydy',GRAPHFILE,5,@coord,@graph_text,@rest); } } sub makegraphpdbdata { my ($gtype,$start,$stop,@list) = @_; my (%list); my (@datalist); for ($it=$start;$it<=$stop;$it++) { foreach $name (@list) { ($min,$avg,$max)=&minavgmax($pdb{$name}[$it]); next if (!defined $avg); #($sd)=&sd($avg,$pdb{$name}[$it]); if ($gtype eq 'xydydy') { $list{$name}.="$it $avg ".($avg-$min)." ".($max-$avg)."\n"; } elsif ($gtype eq 'xy') { $list{$name}.="$it $avg \n"; } } } $color=2; foreach $name (@list) { chomp $list{$name}; next if $list{$name} eq ""; $short=$name; $short=~ s/\w+\(//; $short=~ s/\)$//; push (@datalist,$short,$list{$name},$color); $color++; } return @datalist; } sub print_graph { my ($graphtype,$graphfile,$graphnum,$xmin,$xmax,$ymin,$ymax,$title,$xlabel,$ylabel,@set) = @_; my ($xaxislabel,$yaxislabel,$graphname)=""; my ($i)=0; $graphname="g$graphnum"; $xaxislabel="\@ xaxis label \"$xlabel\"\n" if ($xlabel ne ""); $yaxislabel="\@ yaxis label \"$ylabel\"" if ($ylabel ne ""); $title="\@ title \"$title\"\n" if ($title ne ""); print $graphfile <<GRAPH_TEXT; @ $graphname on @ $graphname label off @ $graphname hidden false @ $graphname type xy @ $graphname autoscale type AUTO @ $graphname fixedpoint off @ $graphname fixedpoint type 0 @ $graphname fixedpoint xy 0.000000, 0.000000 @ $graphname fixedpoint format general general @ $graphname fixedpoint prec 2, 2 @ with $graphname @ world xmin 0 @ world xmax 200 @ world ymin 0 @ world ymax 10000000 @ stack world 0, 0, 0, 0 tick 0, 0, 0, 0 @ view xmin $xmin @ view xmax $xmax @ view ymin $ymin @ view ymax $ymax @ title font 4 @ title size 0.9 @ title color 1 @ title linewidth 1 @ xaxis tick on @ xaxis tick major 1 @ xaxis tick minor 0.5 @ xaxis label layout para @ xaxis label place spec @ xaxis label place 0.00, 0.03 @ xaxis label char size 0.850000 @ xaxis label font 4 @ xaxis label color 1 @ xaxis label linewidth 1 @ xaxis ticklabel on @ xaxis ticklabel char size 0.790000 @ xaxis ticklabel font 4 @ xaxis tick major on @ xaxis tick minor off @ xaxis tick size 0.500000 @ xaxis tick minor size 0.250000 @ xaxis tick op bottom @ yaxis tick on @ yaxis tick major 100 @ yaxis tick minor 50 @ yaxis label layout para @ yaxis label place auto @ yaxis label char size 0.800000 @ yaxis label font 4 @ yaxis label color 1 @ yaxis label linewidth 1 @ yaxis ticklabel on @ yaxis ticklabel char size 0.790000 @ yaxis ticklabel font 4 @ yaxis tick major on @ yaxis tick minor on @ yaxis tick size 0.500000 @ yaxis tick minor size 0.250000 @ yaxis tick op left @ legend char size 0.600000 @ legend font 4 @ legend vgap 2 @ legend hgap 1 @ legend length 2 $title$subtitle$xaxislabel$yaxislabel @ legend on GRAPH_TEXT for ($i=0;$i<=$#set;$i+=3) { if ($graphtype =~ /xydydy/) { print $graphfile "\@ s".($i/3)." errorbar length 0.30000\n"; } print $graphfile "\@ legend x1 ".($xmax-0.15)."\n"; print $graphfile "\@ legend y1 ".($ymax-0.02)."\n"; print $graphfile "\@ legend string ".($i/3)." \"$set[$i]\"\n"; print $graphfile "\@ s".($i/3)." color $set[$i+2]\n"; print $graphfile "\@ s".($i/3)." symbol 2\n"; print $graphfile "\@ s".($i/3)." symbol color $set[$i+2]\n"; print $graphfile "\@ s".($i/3)." symbol fill 1\n"; print $graphfile "\@ s".($i/3)." symbol linewidth 1\n"; print $graphfile "\@ s".($i/3)." symbol size 0.480000\n"; print $graphfile "\@ s".($i/3)." symbol center false\n"; print $graphfile "\@ s".($i/3)." symbol char 0\n"; print $graphfile "\@ s".($i/3)." on\n"; print $graphfile "\@ s".($i/3)." type $graphtype\n"; if ($set[$i+1] ne '') {print $graphfile "$set[$i+1]\n\&\n\@ autoscale\n";} } } sub print_graph_initialize { my ($graphfile) = @_; print $graphfile <<GRAPH_TEXT; @ version 40100 @ page layout free @ ps linewidth begin 1 @ ps linewidth increment 2 @ hardcopy device 2 @ page 5 @ page inout 5 @ page layout portrait @ link page off @ auto redraw on @ default linestyle 1 @ default linewidth 1 @ default color 1 @ default char size 0.6000000 @ default font 3 @ default font source 0 @ default symbol size 0.0500000 @ with string @ string on @ string loctype view @ string 0.5, 0.96 @ string font 4 @ string just 2 @ string char size 1.1 @ string def "Project: $projectname, $run ($pdbbest lowest E structs)" GRAPH_TEXT } ############################## # Subs for printing PS graph # ############################## sub print_pspage { my ($psoutput,$font,$it_firstgraph,$firsttext,$it_secondgraph,$secondtext)=@_; open (PSFILE,">$psoutput")|| die "Could not open $psoutput: $!"; &print_psheader($font); &print_psgraphtext("Weighed ambiguous contact map for project: $projectname ($run)",18,60,400); &print_psgraphtext("(matrices based on $matrixtype)",10,570,100); &print_psgraph($firsttext,14,150,500,40,390,10,$it_firstgraph); &print_psgraph($secondtext,14,150,500,420,770,10,$it_secondgraph); &print_psgraphtext("By user $user on $date.",10,570,630); &print_pstail(); close PSFILE; } sub print_psgraph { my ($text,$text_font_point,$xmin,$xmax,$ymin,$ymax,$font_point,$it) = @_; my ($max,$maxl,$minl) =0; &print_psgraphtext($text,$text_font_point,$xmin-$text_font_point*2,$ymin+($ymax-$ymin)/2); &print_psgraphbox($xmin,$xmax,$ymin,$ymax); &print_psyaxis($xmin,$xmax,$ymin,$ymax,$font_point,@{$contact_dep[$it]}); &print_psxaxis($xmin,$xmax,$ymin,$ymax,$font_point,@{$contact_dep[$it]}); ($minl,$maxl,$max)=&print_psboxes_upper($xmin,$xmax,$ymin,$ymax,$contact_dep_max[$it],$contact_avg[$it],@{$contact_dep[$it]}); &print_psboxes_lower($xmin,$xmax,$ymin,$ymax,@{$contact_map[$it]}); &print_psgraphtext("(max value $max, lower limit $minl, upper shading limit $maxl)",$font_point,$xmax+$font_point*1.5,$ymin+($ymax-$ymin)/2); } sub print_psyaxis { # Attention! Landscape printing so y axis is x coordinates! my ($xmin,$xmax,$ymin,$ymax,$font_point,@pointer) = @_; my ($lineprint,$textprint) =""; print PSFILE ("%% This part does the residues and x axis ticks\n\n"); $xbu = ($xmax-$xmin)/$num_residues; $tick_length = $font_point/5; for ($i=0;$i<$num_residues;$i++) { $xtext=$xmin+$xbu*($i+0.5); $xtick=$xmin+$xbu*($i+1); if ((($i+1)/5) == int (($i+1)/5)) { $lineprint.=sprintf("%3.2f %3.2f %3.2f %3.2f L\n",$xtick,($ymin-$tick_length*1.5),$xtick,$ymax); } else {$lineprint.=sprintf ("%3.2f %3.2f %3.2f %3.2f L\n",$xtick,($ymin-$tick_length),$xtick,$ymin);} if (($i+1) ==1 || $num_residues <= 40 || ($num_residues > 40 && (($i+1)/5) == int (($i+1)/5))) { $textprint.=sprintf ("\n%3.2f %3.2f moveto\n",$xtext,$ymin-$font_point*1.5); $textprint.=sprintf ("(%d) %.1f CenterRot90 Rot90\n(%d) %.1f Print\ngrestore\n",($i+1),$font_point,($i+1),$font_point); } } print PSFILE ("\n0.10 setlinewidth\nCol03\n$lineprint\n"); print PSFILE ("\nCol01\n$textprint\n"); } sub print_psxaxis { # Attention! Landscape printing so x axis is y coordinates! my ($xmin,$xmax,$ymin,$ymax,$font_point,@pointer) = @_; my ($lineprint,$textprint) =""; print PSFILE ("%% This part does the residues and y axis ticks\n\n"); $ybu = ($ymax-$ymin)/$num_residues; $tick_length = $font_point/5; for ($i=0;$i<$num_residues;$i++) { $ytext=$ymin+$ybu*($i+0.5); $ytick=$ymin+$ybu*($i+1); $text=($i+1)-(int(($i+1)/10)*10); if ((($i+1)/5) == int (($i+1)/5)) { $lineprint.=sprintf ("%3.2f %3.2f %3.2f %3.2f L\n",($xmin-$tick_length*1.5),$ytick,$xmax,$ytick); } else {$lineprint.=sprintf ("%3.2f %3.2f %3.2f %3.2f L\n",($xmin-$tick_length),$ytick,$xmin,$ytick);} if ($text == 0) {$text=$i+1;} if ($i == 0 || $num_residues <= 40 || ($num_residues > 40 && ($text/5) == int ($text/5))) { $textprint.=sprintf ("\n%3.2f %3.2f moveto\n",$xmin-$font_point,$ytext); $textprint.=sprintf ("(%d) %.1f CenterRot90 Rot90\n(%d) %.1f Print\ngrestore\n",$text,$font_point,$text,$font_point); } } print PSFILE ("\n0.10 setlinewidth\nCol03\n$lineprint\n"); print PSFILE ("\nCol01\n$textprint\n"); } sub print_psboxes_upper { my ($xmin,$xmax,$ymin,$ymax,$max,$avg,@pointer) = @_; my ($xbu,$ybu,$limit_max,$limit_min,$factor,$scale,$log_min)=0; $xbu = ($xmax-$xmin)/($num_residues); $ybu = ($ymax-$ymin)/($num_residues); if ($max==0) {$max=0.5;} $factor=($max/$avg)*10/$num_residues; $limit_max=$max/$factor; $limit_min=$avg/$factor; $log_min=log($limit_min+1); $scale=(log($limit_max+1)-$log_min)/9; print PSFILE ("\n0.10 setlinewidth\n"); for $k (1...$num_residues) { for $l ($k...$num_residues) { $xlow=$xmin + ($xbu*($k-1)); $xhigh=$xmin + ($xbu*$k); $ylow=$ymin + ($ybu*($l-1)); $yhigh=$ymin + ($ybu*$l); if ($pointer[$k]->[$l] > $limit_max) {$color=9;} elsif ($pointer[$k]->[$l] < $limit_min) {next;} else {$color = sprintf ("%.0f",((log($pointer[$k]->[$l]+1)-$log_min)/$scale));} print PSFILE ("Col1$color\n"); printf PSFILE ("%5.2f %5.2f %5.2f %5.2f ",$xlow,$ylow,$xlow,$yhigh); printf PSFILE ("%5.2f %5.2f %5.2f %5.2f Pline4\n",$xhigh,$yhigh,$xhigh,$ylow); } } return ((sprintf "%.2f",$limit_min),(sprintf "%.2f",$limit_max),(sprintf "%.1f",$max)); } sub print_psboxes_lower { my ($xmin,$xmax,$ymin,$ymax,@pointer) = @_; $xbu = ($xmax-$xmin)/$num_residues; $ybu = ($ymax-$ymin)/$num_residues; print PSFILE ("\n0.10 setlinewidth\n"); for $k (1...$num_residues) { for $l (1...($k-1)) { $xlow=$xmin + ($xbu*($k-1)); $xhigh=$xmin + ($xbu*$k); $ylow=$ymin + ($ybu*($l-1)); $yhigh=$ymin + ($ybu*($l)); if ($pointer[$l]->[$k] < 1) {next;} else {$color = 5;} print PSFILE ("Col0$color\n"); printf PSFILE ("%5.2f %5.2f %5.2f %5.2f ",$xlow,$ylow,$xlow,$yhigh); printf PSFILE ("%5.2f %5.2f %5.2f %5.2f Pline4\n",$xhigh,$yhigh,$xhigh,$ylow); } } print PSFILE ("\nCol01\n"); } # Prints the ps header. sub print_psheader { my ($font)=@_; print PSFILE <<END_PS; %!PS-Adobe-3.0 %%Creator: $user %%DocumentNeededResources: font $font Symbol %%BoundingBox: (atend) %%Pages: 1 %%Title: $scriptname_short graph %%EndComments %%BeginProlog /L { moveto lineto stroke } bind def /Col { setrgbcolor } bind def /Col01 {gsave 0.0000 0.0000 0.0000 setrgbcolor } def /Col02 {gsave 1.0000 1.0000 1.0000 setrgbcolor } def /Col03 {gsave 0.2000 0.2000 0.2000 setrgbcolor } def /Col05 {gsave 0.5000 0.5000 0.5000 setrgbcolor } def /Col09 {gsave 0.9000 0.9000 0.9000 setrgbcolor } def /Col10 {gsave 1.0000 1.0000 1.0000 setrgbcolor } def /Col11 {gsave 1.0000 0.9000 0.9000 setrgbcolor } def /Col12 {gsave 1.0000 0.8000 0.8000 setrgbcolor } def /Col13 {gsave 1.0000 0.7000 0.7000 setrgbcolor } def /Col14 {gsave 1.0000 0.6000 0.6000 setrgbcolor } def /Col15 {gsave 1.0000 0.5000 0.5000 setrgbcolor } def /Col16 {gsave 1.0000 0.4000 0.4000 setrgbcolor } def /Col17 {gsave 1.0000 0.3000 0.3000 setrgbcolor } def /Col18 {gsave 1.0000 0.2000 0.2000 setrgbcolor } def /Col19 {gsave 1.0000 0.1000 0.1000 setrgbcolor } def /Poly4 { moveto lineto lineto lineto fill grestore } bind def /Pline4 { 8 copy Poly4 moveto lineto lineto lineto closepath stroke } bind def /Print { /$font findfont exch scalefont setfont show } bind def /Gprint { /Symbol findfont exch scalefont setfont show } bind def /Center { dup /$font findfont exch scalefont setfont exch stringwidth pop -2 div exch -3 div rmoveto } bind def /CenterRot90 { dup /$font findfont exch scalefont setfont exch stringwidth pop -2 div exch 3 div exch rmoveto } bind def /UncenterRot90 { dup /$font findfont exch scalefont setfont exch stringwidth } bind def /Rot90 { gsave currentpoint translate 90 rotate } bind def %%EndProlog %%BeginSetup 1 setlinecap 1 setlinejoin 1 setlinewidth 0 setgray [ ] 0 setdash newpath %%EndSetup %%Page: 1 1 1 1 scale 0 0 translate /psoverview save def closepath gsave 1.0000 setgray fill grestore stroke gsave 0.10 setlinewidth END_PS } sub print_psgraphbox { my ($xmin,$xmax,$ymin,$ymax)=@_; print PSFILE <<END_PS; %% This part draws the graph box 0.80 setlinewidth Col09 $xmin $ymin $xmin $ymax $xmax $ymax $xmax $ymin Pline4 END_PS } sub print_psgraphtext { my ($Graph_title,$font_point,$x,$y)=@_; print PSFILE <<END_PS; %% This part does the text $x $y moveto ($Graph_title) $font_point CenterRot90 Rot90 ($Graph_title) $font_point Print grestore END_PS } sub print_pstail { print PSFILE <<END_PS; psoverview restore showpage %%Trailer %%BoundingBox: 0 0 612 792 %%EOF END_PS }