tools/gc_dep_plot/old/gc_vs_depth.pl

   1 #!/usr/bin/perl -w
   2 #use strict;
   3 use Time::HiRes qw ( gettimeofday tv_interval );
   4 use Getopt::Long;
   5 my ($depth_file,$fa_file,$out_file,$win_size,$help);
   6 GetOptions(
   7         "depth:s"=>\$depth_file,
   8         "fa:s"=>\$fa_file,
   9         "output:s"=>\$out_file,
  10         "win:i"=>\$win_size,
  11         "help"=>\$help
  12 );
  13 if(!$depth_file || !$fa_file || !$out_file || !$win_size || $help)
  14 {
  15         print "\t -depth: the depth file\n";
  16         print "\t -fa: the .fa file\n";
  17         print "\t -output: the output file\n";
  18         print "\t -win: the window size\n";
  19         print "\t -help: if you need the help";
  20         die;
  21 }
  22 #set time
  23 my $str_time=[gettimeofday];
  24 #getDepth
  25 print "Begin to count average of depth\n";
  26 my %depth=();
  27 my $cnt=0;
  28 my $sum=0;
  29 my $chr;
  30 open(IN,"<$depth_file")or die "can not open $depth_file :$!\n";
  31 while(<IN>)
  32 {
  33         chomp;
  34         if($_!~/^>/)
  35         {
  36                 while($_=~/(\d+)/g)
  37                 {
  38                         my $num=$1;
  39                         if($num==65535){$num=0;}
  40                         $sum+=$num;
  41                         $cnt++;
  42                         if($cnt==$win_size)
  43                         {
  44                                 my $aver=$sum/$cnt;
  45                                 push @{$chr},$aver;
  46                                 $sum=0;
  47                                 $cnt=0;
  48                         }
  49                 }
  50         }
  51         else
  52         {
  53                 $chr=$_;
  54                 $depth{$chr}=$chr;
  55                 print "$chr\n";
  56         }
  57 }
  58 close IN;
  59 close OUT;
  60 print "Finish the depth counting\n";
  61
  62 #countGC
  63 $chr=();
  64 $cnt=0;
  65 $sum=0;
  66 print "Begin to count the gc\n";
  67 open(IN,"<$fa_file")or die "can not open $fa_file:$!";
  68 while(<IN>)
  69 {
  70         chomp;
  71         if($_!~/^>/)
  72         {
  73                 my $num=s/[gc]/Y/ig;
  74                 $sum+=$num;
  75                 $cnt+=length($_);
  76                 if($cnt==$win_size)
  77                 {
  78                         my $aver=int($sum/$cnt*100);
  79                         push @{$chr},$aver;
  80                         $sum=0;
  81                         $cnt=0;
  82                 }
  83                 elsif($cnt>$win_size)
  84                 {
  85                         my $dis=$cnt-$win_size;
  86                         my @chr=split(//,$_);
  87                         my $len=@chr;
  88                         my $count=0;
  89                         for(my $i=$len-$dis;$i<$len;$i++)
  90                         {
  91                                 if($chr[$i]=~/[gc]/i)
  92                                 {
  93                                         $count++;
  94                                 }
  95                         }
  96                         my $aver=int(($sum-$count)/$win_size*100);
  97                         push @{$chr},$aver;
  98                         $sum=$count;
  99                         $cnt=$dis;
 100                 }
 101                 else {next;}
 102         }
 103         else
 104         {
 105                 $chr="gc$_";
 106                 print "$_\n";
 107         }
 108 }
 109 close IN;
 110
 111 foreach(keys %depth)
 112 {
 113         my $chrname=$_;
 114         foreach(@{"gc$chrname"})
 115         {
 116                 my $dep=shift @{$chrname};
 117                 push @{$_},$dep;
 118         }
 119 }
 120
 121 my $str=0;
 122 my $end=80;
 123 my $na=0;
 124 my $line=$end-$str;
 125 my $max_depth=0;
 126 my $cnt_depth=0;
 127 open(OUT,">$out_file")or die "$!\n";
 128 for(my $i=$str;$i<=$end;$i++)
 129 {
 130         print OUT "$i\t";
 131 }
 132 print OUT "\n";
 133 while($na<=$line)
 134 {
 135         $na=0;
 136         for(my $j=$str;$j<=$end;$j++)
 137         {
 138                 my $dep=shift @{$j};
 139                 if(!defined $dep){$dep="NA";$na++;}
 140                 else{$max_depth+=$dep;$cnt_depth++;}
 141                 print OUT "$dep\t";
 142         }
 143         print OUT "\n";
 144 }
 145 close OUT;
 146 #R plot
 147 $max_depth/=$cnt_depth/2;
 148 $max_depth=(int($max_depth/50)+1)*50;
 149 my $area=0;
 150 foreach(1..$max_depth/50)
 151 {
 152         my $ev_depth=$_*50;
 153         $area.=",$ev_depth";
 154 }
 155 open(OUT,">$out_file.R")or die;
 156 print OUT 'a=read.delim("',"$out_file",'");',"\n";
 157 print OUT 'pdf("',"$out_file.pdf",'",width=8,height=6,paper="special");',"\n";
 158 print OUT "boxplot(a,ylim=c(0,$max_depth),xlim=c(0,90),axes=FALSE,outline=FALSE,xlab='GC%',ylab='Depth');\n";
 159 print OUT "axis(1,at=c(0,10,20,30,40,50,60,70,80),labels=c(0,10,20,30,40,50,60,70,80),las=0);\n";
 160 print OUT "axis(2,at=c($area),labels=c($area),las=2);\nq();\nn\n";
 161 close OUT;
 162 `Rscript $out_file.R`;
 163 my $end_time=[gettimeofday];
 164 my $time=tv_interval($str_time,$end_time);