CD task10 protocol

From Bioinformatikpedia
Revision as of 16:33, 14 August 2012 by Vorbergs (talk | contribs) (Significant Differences of RMSF)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)



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>