released/pIRS.old/stat_illumina_reads/plotdat.pl

   1 #!/bin/env perl
   2 =pod
   3 Author: Hu Xuesong @ BGI <huxuesong@genomics.org.cn>
   4 Version: 0.1.1 @ 20111111
   5 =cut
   6 use strict;
   7 use warnings;
   8
   9 die "Usage: $0 <datfile> [max-depth-to-plot]\n" unless @ARGV;
  10 my $name=$ARGV[0];
  11 my $maxdep=$ARGV[1];
  12
  13 my $WinSize=0;
  14 my (%GCdat,@GCvalues);
  15 open IN,'<',$name or die "Error openimg $name: $!\n";
  16 while(<IN>) {
  17     if (/^#WinSize=(\d+)/) {
  18         $WinSize=$1;
  19         next;
  20     }
  21     next if /^(#|$)/;
  22     chomp;
  23     my ($gc,undef,$depcnt,$mean)=split /\s+/;
  24     $gc -= 0.5;
  25     $GCdat{$gc}=[$depcnt,$mean];
  26     #print "{$gc}=[$depcnt,$mean]\n";
  27     push @GCvalues,$gc;
  28 }
  29 close IN;
  30
  31 open OUT,'>',$name.'.gc' or die "Error openimg ${name}.gc: $!\n";
  32 print OUT "#WinSize: $WinSize\n",join("\t",'#GC%(left border)','MeanDepth'),"\n";
  33 my $step=5; # 5%
  34 my $lastGC=0;
  35 my ($summer,$counter)=(0,0);
  36 for my $gc (0..100) {
  37     if ($gc-$lastGC>=$step) {
  38         print OUT join("\t",$lastGC,$summer/$counter),"\n" if $counter;
  39         $lastGC=$gc;
  40         ($summer,$counter)=(0,0);
  41     }
  42     next unless exists $GCdat{$gc};
  43     $summer += $GCdat{$gc}->[1] * $GCdat{$gc}->[0];
  44     $counter += $GCdat{$gc}->[0];
  45 }
  46 #print OUT join("\t",$GCvalues[-1],$summer/$counter),"\n" if $counter;
  47 close OUT;
  48
  49 chomp(my $user=`id -nru`);
  50 my ($gnuplot,$font);
  51 if ($user eq 'galaxy') {
  52     $gnuplot='gnuplot';
  53     $font='/home/galaxy/fonts/arial.ttf';
  54 } else {
  55     $gnuplot='/opt/blc/genome/biosoft/gnuplot-4.4.0/bin/gnuplot';
  56     $font='/ifs1/ST_ASMB/USER/huxuesong/public/fonts/arial.ttf';
  57 }
  58 open P,'>',$name.'.dem' or die "Error openimg $name.dem: $!\n";
  59 my $yrange='#set yrange [0:10]';
  60 if ($maxdep) {
  61         $yrange="set yrange [0:$maxdep]";
  62 }
  63 print P <<"CODE";
  64 reset
  65 set xrange [-1:101]
  66 $yrange
  67 #set y2range [0:100]
  68 set border 11 lw 2
  69 set xtics 0,5 nomirror
  70 set ytics nomirror format "% 4.0f"
  71 set y2tics nomirror
  72
  73 #set term png font "$font" 64 size 4204,3154 truecolor linewidth 3
  74 set term png font "$font" 16 size 1600,1200
  75 set output "$name.png"
  76 set multiplot layout 2,1
  77 set title "GC-dep plot of Window_Size=$WinSize"
  78 set xlabel "GC %"
  79 set ylabel "Depth"
  80 set y2label "Max Depth"
  81 set boxwidth 0.75
  82 set style fill empty
  83 set origin 0,0.21
  84 set size 1,0.79
  85 set logscale y2
  86 plot '$name' \\
  87      using 1:11 with lines axis x1y2 lt rgb "#006400" lw 2 title 'max depth', \\
  88      ''      using 1:6:5:9:8 with candlesticks lt rgb "navy" lw 2 title 'depth box' whiskerbars 0.5, \\
  89      ''      using 1:7:7:7:7 with candlesticks lt -1 lw 2 notitle, \\
  90      ''      using 1:4 with points lt rgb "red" lw 3 title 'mean depth'
  91
  92 set y2label " "
  93 set y2tics textcolor rgb "#FFFFFF"
  94 set title ""
  95 set xtics out nooffset format ''
  96 set xlabel ""
  97 set ylabel ""
  98 set yrange [*:*]
  99 set logscale y
 100 set ytics nomirror format "% 4.0f" 1,10,1e5
 101 set size 1,0.21
 102 plot '$name' using 1:3 with boxes fs solid 0.62 lt rgb "navy" lw 1.5 title 'window count'
 103 unset multiplot
 104 #pause -1 "Hit return to continue"
 105 CODE
 106 close P;
 107
 108 system($gnuplot,$name.'.dem');
 109 __END__
 110