YARM samples:

Statistical analysis: per residue

Statistically compare the NOE volumes of a simulated data set from a Model and an experimentally measured NOE volume data set per residue. This script returns four statistical value, RMS, R-factor, Q-factor and Q^(1/6)-factor, see the Stats page for a description of the statistical functions.

#!/usr/local/bin/perl 
# yarm Yet Another Relaxation Matrix program

# Read in the library file:
require "/usr/local/yarm/yarm_lib.pl";

#############################################
# Define variables
#############################################
$pdb_file = "dick_xtal.pdb";
$sfreq = 600;
$vol0 = 100;
$tmix = 0.2;

$tl = 2.5;
$ts = 6;
$tc = 3;

$f95_vol_file = "d12_30s.vols";
$f95_peak_file = "d12_30s.peaks";

# Set to 0 for no debugging messages
# Set to 1 for a few debugging messages
# set to 2 for TONS of debugging messages (lots!)
$debug = 0; 

#############################################
# Call YARM subroutines
#############################################
print "$yarm_version\n";

# Get non-exchangeable nucleic acid protons from a PDB file
%xyz = &Pdb_Read_All( $pdb_file );
%xyz = &Get_Atom_Type( \%xyz, \%nonX_NA );
%xyz = &Pseudo_Methyl( \%xyz );

# Measure distances of all atom pairs between 
# 0 to 10 angstroms
%rij  = &Rij_Hash( \%xyz, 0, 10 );

# Read in experimental volumes
%vol_exp = &F95_Read_Merge( $f95_vol_file, $f95_peak_file );
%vol_exp = &Make_Symm_Molecule( A, B, \%vol_exp );

# ANISOTROPIC ROTATION
# Calculate the principal axis of rotation vector
print "Simulating NOE volumes using anisotropic-rigid...\n";
( $Ax, $Ay, $Az ) = Principal_Axis( \%xyz );
printf ("Principal axis vector components Ax=%4.2f Ay=%4.2f Az=%4.2f\n", $Ax, $Ay, $Az);
%vol_sim = &Sim_Vol( $sfreq, $tmix, $vol0, \%xyz, \%rij, $tl, $ts, $Ax, $Ay, $Az );

# ISOTROPIC ROTATION
# print "Simulating NOE volumes using isotropic-rigid...\n";
# %vol_sim = &Sim_Vol( $sfreq, $tmix, $vol0, \%xyz, \%rij, $tc );

# Normalize the experimental volumes by comparing h5-h6 volumes
%vol_exp = &Norm_Hash( \%vol_exp, \%vol_sim);

print "Results are save in temp.report in a more xmgr-able form\n\n";
open ( REPORT, "> temp.report" );

# This is assuming 2 segids named A and B for a DNA with 12 basepairs
print "Segid Num Noe  RMS     R      Q    Q^(1/6)\n";
print "----- --- --- ------ ------ ------ ------\n";
foreach $segid ( A, B ) {
	foreach $res_num ( 1 .. 12 ) {
		# This one returns only intranucleotide NOE info
		# %residue = &Subset( "^$segid $res_num ", "^$segid $res_num ", \%vol_exp );

		# This one returns intra and inter-nucleotide NOE info
		# The "period" for the second match means "anything"
		%residue = &Subset( "^$segid $res_num ", ".", \%vol_exp );

		($rms, $r, $q, $q6) = &Stats( \%vol_exp, \%vol_sim, \%residue );

		@protons = keys %residue;
		
		printf (" %4s %3d %3d %5.4f %5.4f %5.4f %5.4f\n",
			$segid, $res_num, $#protons+1, $rms, $r, $q, $q6);
		
		print REPORT "$res_num $rms $r $q $q6\n";
	}
}
exit;

close REPORT;