perl/etc/itilestator.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <read1.fq.gz> <out>\n" if @ARGV != 2;
   6 my ($fq1,$out)=@ARGV;
   7
   8 sub openfile($) {
   9     my ($filename)=@_;
  10     my $infile;
  11     if ($filename=~/.bz2$/) {
  12             open( $infile,"-|","bzip2 -dc $filename") or die "Error opening $filename: $!\n";
  13     } elsif ($filename=~/.gz$/) {
  14         open( $infile,"-|","gzip -dc $filename") or die "Error opening $filename: $!\n";
  15     } else {open( $infile,"<",$filename) or die "Error opening $filename: $!\n";}
  16     return $infile;
  17 }
  18 sub readfq($) {
  19         my $fh=$_[0];
  20         defined(my $a=<$fh>) or return [];
  21         chomp($a);
  22         chomp(my $b=<$fh>) or return [];
  23         chomp(my $c=<$fh>) or return [];
  24         chomp(my $d=<$fh>) or return [];
  25         return [$a,$b,$c,$d];
  26 }
  27 sub getQ($) {
  28     my @Qstr=split //,$_[0];
  29     my @Qvalue=();
  30     push @Qvalue,ord($_)-64 for @Qstr;
  31     return \@Qvalue;
  32 }
  33 sub statQ($) {
  34     my $Qvalues=$_[0];
  35     my $Qlen=scalar @$Qvalues;
  36     die "[x]Read Length must >= 50." if $Qlen<50;
  37     my ($Q1,$Q5,$Q10,$Q15,$Qs15,$Q30,$Q50,$Qall);
  38     $Q1=$$Qvalues[0];
  39     $Q5+=shift @$Qvalues for (1..5);
  40     #$Q10=$Q15=$Qs15=$Q30=$Q50=$Qall=$Q5;
  41     $Q10=$Q5;
  42     $Q10+=shift @$Qvalues for (1..5);
  43     $Q15=$Q10;
  44     $Q15+=shift @$Qvalues for (1..5);
  45     $Qs15+=shift @$Qvalues for (1..15);
  46     $Q30=$Q15+$Qs15;
  47     $Q50=$Q30;
  48     $Q50+=shift @$Qvalues for (1..20);
  49     $Qall=$Q50;
  50     $Qall+=shift @$Qvalues while (@$Qvalues);
  51     return [$Q1,$Q5/5,$Q10/10,$Q15/15,$Qs15/15,$Q30/30,$Q50/50,$Qall/$Qlen];
  52 }
  53 my @statQlabel=qw/Q1 Q5 Q10 Q15 Q16t30 Q30 Q50 Qall/;
  54
  55 my $FQa=openfile($fq1);
  56 #my $FQb=openfile($fq2);
  57
  58 my (%ReadsinTile,%Xrange,%Yrange,%statQ);
  59
  60 my $fqitem=&readfq($FQa);
  61 while (@$fqitem > 0) {
  62     my ($id,$Q)=@$fqitem[0,3];
  63     my ($Tile,$X,$Y)=(split /[:#]/,$id)[2,3,4];
  64     my $Qvalues=&getQ($Q);
  65     my $Qstat=&statQ($Qvalues);
  66     #my ($Surface,$Swath,$subTile)=(split //,$Tile)[0,1,3];
  67     #print "$Surface,$Swath,$subTile\t$Tile,$X,$Y\n";
  68     ++$ReadsinTile{$Tile};
  69     $statQ{$Tile}=[] unless exists $statQ{$Tile};
  70     my $i=0;
  71     for my $Qv (@$Qstat) {
  72         ${$statQ{$Tile}}[$i] += $Qv;
  73         ++$i;
  74     }
  75     ++$Xrange{$X}; ++$Yrange{$Y};
  76     $fqitem=&readfq($FQa);
  77 }
  78 close $FQa;
  79
  80 open OUT,'>',$out or die "Error opening $out: $!\n";
  81 for my $tile (sort keys %ReadsinTile) {
  82     print OUT "Tile\t$tile\t$ReadsinTile{$tile}\n";
  83 }
  84 print OUT "#Qstat\tTile\t",join("\t",@statQlabel),"\n";
  85 for my $tile (sort keys %statQ) {
  86     $_ /=$ReadsinTile{$tile} for @{$statQ{$tile}};
  87     print OUT "Qstat\t$tile\t",join("\t",@{$statQ{$tile}}),"\n";
  88 }
  89 print OUT "X\t$_\t$Xrange{$_}\n" for sort {$a<=>$b} keys %Xrange;
  90 print OUT "Y\t$_\t$Yrange{$_}\n" for sort {$a<=>$b} keys %Yrange;
  91 close OUT;
  92
  93
  94