Difference between revisions of "CD task10 protocol"

From Bioinformatikpedia
(→‎Significant Differences of RMSF)
 
(13 intermediate revisions by the same user not shown)
Line 1: Line 1:
−
==Script with all analysis steps==
  
−
##!/bin/bash
 
   
−
#wildtype
  
−
path=/opt/SS12-Practical/gromacs/bin/
  
−
EDR="2O4H_chainA_md"
  
−
MUT="wt"
  
   
   
−
#gmxcheck
  
−
$path gmxcheck -c $EDR.tpr
 
−
gmxcheck -f $EDR.xtc
  
−  
−
#pymol visualization
  
−
echo 1 0| trjconv -s $EDR.tpr -f $EDR.xtc -o $MUT_vis.pdb -pbc nojump -dt 10
  
   
  +
==Script with all analysis steps==
   
  +
In this script, all the required commands for the analysis methods are listed.
   
  +
[https://www.dropbox.com/s/jse4hqskne4t9tz/analysis.sh analysis.sh]
−
#energy
  
−
echo 12 0| g_energy -f $EDR.edr -o $MUT.temp.xvg
  
−
echo 13 0 | g_energy -f $EDR.edr -o $MUT.pressure.xvg
  
−
echo 9 0| g_energy -f $EDR.edr -o $MUT.potential.xvg
  
−
echo 11 0 | g_energy -f $EDR.edr -o $MUT.totenergy
  
   
  +
==Pymol visualisation ==
   
  +
simple Pymol visalization script
−
#min dist between periodic boundaries
  
−
echo 1 0| g_mindist -f $EDR.xtc -s $EDR.tpr -od $MUT.minimal-periodic-distance.xvg -pi
  
−
echo 3 0| g_mindist -f $EDR.xtc -s $EDR.tpr -od $MUT.minimal-periodic-distance_calpha.xvg -pi
  
   
  +
[https://www.dropbox.com/s/5kn68awfdfwl5xe/create_pymol_vis.py create_pymol_vis.py]
−
#fluctuations
  
   
−
echo 1 0| g_rmsf -f $EDR.xtc -s $EDR.tpr -o $MUT.rmsf-per-residue.xvg -ox $MUT.average.pdb -oq $MUT.bfactors.pdb -res
 
   
−
#rewrite trajectories
  
−
trjconv -f $EDR.xtc -o $MUT.traj_nojump.xtc -pbc nojump
 
   
  +
mencoder "mf://*.png" -mf w=640:h=480:type=png -ofps 25 -ovc lavc
−
#calculate rmsd
  
  +
-lavcopts vcodec=mpeg4:vpass=1:vbitrate=5000000:mbd=2:keyint=132:v4mv:vqmin=3:lumi_mask=0.07:dark_mask=0.2:scplx_mask=0.1:tcplx_mask=0.1:naq
−
g_rms -f $MUT.traj_nojump.xtc -s $EDR.tpr -o $MUT.rmsd-all-atom-vs-start.xvg
 
  +
-o wt_movie4.avi
   
  +
== Significant Differences of RMSF==
−
#rmsd only calpha
  
−
echo 4 4 0 | g_rms -f $MUT.traj_nojump.xtc -s $EDR.tpr -o $MUT.rmsd-backbone-vs-start.xvg
 
   
  +
Provided by Marc
  +
<source lang="perl">
   
  +
#!/usr/bin/perl -w
−
#rmsds vs average structure
  
  +
$swExit = 0;
−
#####
  
−
g_rms -f $MUT.traj_nojump.xtc -s $EDR.tpr -o $MUT.protein.xtc
 
−
##select1
  
−
##select4
  
   
  +
if(@ARGV < 2){
−
#radius of gyration
  
  +
print STDERR "statTD.pl <RMSF FILE1> <RMSF FILE2> [START] [END]\n";
−
g_gyrate -f $EDR.xtc -s $EDR.tpr -o $MUT.radius-of-gyration.xvg
 
  +
exit;
  +
}
  +
else{
  +
if(! -e "$ARGV[0]"){
  +
print STDERR "File $ARGV[0] does not exist!\n";
  +
$swExit = 1;
  +
}
  +
if(! -e "$ARGV[1]"){
  +
print STDERR "File $ARGV[1] does not exist!\n";
  +
$swExit = 1;
  +
}
  +
}
  +
if(defined $ARGV[2]){
   
  +
if(not defined $ARGV[3]){
  +
print STDERR "END value missing!\n";
  +
$swExit = 1;
  +
}
  +
elsif($ARGV[2] !~ m/[0-9]+/ || $ARGV[3] !~ m/[0-9]+/){
  +
print STDERR "Problem with the ranges please check!\n";
  +
$swExit = 1;
  +
}
  +
}
  +
if($swExit){exit;}
  +
if(defined $ARGV[2] && defined $ARGV[3] && $ARGV[2] > $ARGV[3]){
  +
$tmp = $ARGV[2];
  +
$ARGV[2] = $ARGV[3];
  +
$ARGV[3] = $tmp;
  +
}
   
   
  +
open(FILE,$ARGV[0]);
  +
while(<FILE>){
  +
chomp($_);
  +
if(/^\s*([0-9]+)\s+([0-9]+\.[0-9]+)/){
  +
#print "$1 $2\t";
  +
if(defined $ARGV[2] && defined $ARGV[3] && $1 >= $ARGV[2] && $1 <= $ARGV[3]){
  +
$in1{$1} = $2;
  +
#print "1\n";
  +
}
  +
elsif(not defined $ARGV[2]){
  +
$in1{$1} = $2;
  +
#print "$ARGV[2] $ARGV[3] 3\n";
  +
}
  +
  +
  +
}
  +
}
  +
close FILE;
  +
open(FILE,$ARGV[1]);
  +
while(<FILE>){
  +
chomp($_);
  +
if(/^\s*([0-9]+)\s+([0-9]+\.[0-9]+)/){
  +
#print "$1 $2\n";
  +
if(defined $ARGV[2] && defined $ARGV[3] && $1 >= $ARGV[2] && $1 <= $ARGV[3]) {
  +
$in2{$1} = $2;
  +
}
  +
elsif(not defined $ARGV[2]){
  +
$in2{$1} = $2;
  +
}
  +
}
  +
}
  +
close FILE;
  +
if((keys %in1) != (keys %in2)){
  +
print STDERR "Number of datapoints not identical!\n";
  +
exit;
  +
}
   
  +
$value = &getPV(\%in1,\%in2);
−
#a305e
  
   
  +
print "$value";
−
EDR="A305E_scwrl_mini_md"
  
  +
if($value <= 0.001){
−
MUT="A305E_scwrl"
  
  +
print "**\n";
  +
}
  +
elsif($value <= 0.05){
  +
print "*\n";
  +
}
  +
else{
  +
print "-\n";
  +
}
  +
sub getPV{
  +
my (%data1) = %{$_[0]};
  +
my (%data2) = %{$_[1]};
  +
my $avg = 0;
  +
my $cnt = 0;
  +
my $std = 0;
  +
my $ste = 0;
  +
my $ratio = 0;
  +
  +
foreach $key (sort keys %data1){
  +
#print "$key\n";
  +
$avg += ($data1{$key} - $data2{$key});
  +
$cnt++;
  +
}
  +
$avg /= $cnt;
  +
foreach $key (sort keys %data1){
  +
$std += (($data1{$key} - $data2{$key}) - $avg)**2;
  +
}
  +
$std /= $cnt;
  +
$std = sqrt($std);
  +
$ste = $std/sqrt($cnt);
  +
$ratio = abs($avg/$ste);
  +
open(OUT,">/tmp/rscript.tmp");
  +
print OUT "ratio <- c($ratio);\n";
  +
print OUT "dfree <- c($cnt);\n";
  +
print OUT "dt(ratio,dfree);\n";
  +
close OUT;
  +
$value = `Rscript /tmp/rscript.tmp`;
  +
chomp($value);
  +
`rm -f /tmp/rscript.tmp`;
  +
return (split/\s+/,$value)[1];
  +
}
   
  +
</source>
−
echo 12 0| g_energy -f $EDR.edr -o $MUT.temp.xvg
  
−
echo 13 0 | g_energy -f $EDR.edr -o $MUT.pressure.xvg
  
−
echo 9 0| g_energy -f $EDR.edr -o $MUT.potential.xvg
  
−
echo 11 0 | g_energy -f $EDR.edr -o $MUT.totenergy.xvg
  

Latest revision as of 16:33, 14 August 2012



Script with all analysis steps

In this script, all the required commands for the analysis methods are listed.

analysis.sh

Pymol visualisation

simple Pymol visalization script

create_pymol_vis.py 


mencoder "mf://*.png" -mf w=640:h=480:type=png -ofps 25 -ovc lavc 
-lavcopts vcodec=mpeg4:vpass=1:vbitrate=5000000:mbd=2:keyint=132:v4mv:vqmin=3:lumi_mask=0.07:dark_mask=0.2:scplx_mask=0.1:tcplx_mask=0.1:naq
-o wt_movie4.avi

Significant Differences of RMSF

Provided by Marc <source lang="perl">

  1. !/usr/bin/perl -w

$swExit = 0;

if(@ARGV < 2){ print STDERR "statTD.pl <RMSF FILE1> <RMSF FILE2> [START] [END]\n"; exit; } else{ if(! -e "$ARGV[0]"){ print STDERR "File $ARGV[0] does not exist!\n"; $swExit = 1; } if(! -e "$ARGV[1]"){ print STDERR "File $ARGV[1] does not exist!\n"; $swExit = 1; } } if(defined $ARGV[2]){

if(not defined $ARGV[3]){ print STDERR "END value missing!\n"; $swExit = 1; } elsif($ARGV[2] !~ m/[0-9]+/ || $ARGV[3] !~ m/[0-9]+/){ print STDERR "Problem with the ranges please check!\n"; $swExit = 1; } } if($swExit){exit;} if(defined $ARGV[2] && defined $ARGV[3] && $ARGV[2] > $ARGV[3]){ $tmp = $ARGV[2]; $ARGV[2] = $ARGV[3]; $ARGV[3] = $tmp; }


open(FILE,$ARGV[0]); while(<FILE>){ chomp($_); if(/^\s*([0-9]+)\s+([0-9]+\.[0-9]+)/){ #print "$1 $2\t"; if(defined $ARGV[2] && defined $ARGV[3] && $1 >= $ARGV[2] && $1 <= $ARGV[3]){ $in1{$1} = $2; #print "1\n"; } elsif(not defined $ARGV[2]){ $in1{$1} = $2; #print "$ARGV[2] $ARGV[3] 3\n"; }


} } close FILE; open(FILE,$ARGV[1]); while(<FILE>){ chomp($_); if(/^\s*([0-9]+)\s+([0-9]+\.[0-9]+)/){ #print "$1 $2\n"; if(defined $ARGV[2] && defined $ARGV[3] && $1 >= $ARGV[2] && $1 <= $ARGV[3]) { $in2{$1} = $2; } elsif(not defined $ARGV[2]){ $in2{$1} = $2; } } } close FILE; if((keys %in1) != (keys %in2)){ print STDERR "Number of datapoints not identical!\n"; exit; }

$value = &getPV(\%in1,\%in2);

print "$value"; if($value <= 0.001){ print "**\n"; } elsif($value <= 0.05){ print "*\n"; } else{ print "-\n"; } sub getPV{ my (%data1) = %{$_[0]}; my (%data2) = %{$_[1]}; my $avg = 0; my $cnt = 0; my $std = 0; my $ste = 0; my $ratio = 0;

foreach $key (sort keys %data1){ #print "$key\n"; $avg += ($data1{$key} - $data2{$key}); $cnt++; } $avg /= $cnt; foreach $key (sort keys %data1){ $std += (($data1{$key} - $data2{$key}) - $avg)**2; } $std /= $cnt; $std = sqrt($std); $ste = $std/sqrt($cnt); $ratio = abs($avg/$ste); open(OUT,">/tmp/rscript.tmp"); print OUT "ratio <- c($ratio);\n"; print OUT "dfree <- c($cnt);\n"; print OUT "dt(ratio,dfree);\n"; close OUT; $value = `Rscript /tmp/rscript.tmp`; chomp($value); `rm -f /tmp/rscript.tmp`; return (split/\s+/,$value)[1]; }

</source>

Retrieved from "https://i12r-studfilesrv.informatik.tu-muenchen.de/wiki/index.php?title=CD_task10_protocol&oldid=27453"