projects/recombineMap/sperm_by_fw/snp_calling/matrix_pirs2soapsnp.pl

   1 #!/usr/bin/perl
   2
   3 =head1 Name
   4
   5 matrix_pirs2soapsnp.pl   --  convert pirs matrix to soapsnp format
   6
   7 =head1 Description
   8
   9
  10
  11 =head1 Version
  12
  13   Author: Fan Wei, fanw@genomics.org.cn
  14   Version: 1.0,  Date: 2006-12-6
  15   Note:
  16
  17 =head1 Usage
  18
  19   matrix_pirs2soapsnp.pl <*pirs.matrix>
  20   --help      output help information to screen
  21
  22 =head1 Exmple
  23
  24
  25
  26 =cut
  27
  28 use strict;
  29 use Getopt::Long;
  30 use FindBin qw($Bin $Script);
  31 use File::Basename qw(basename dirname);
  32 use Data::Dumper;
  33 use File::Path;  ## function " mkpath" and "rmtree" deal with directory
  34
  35 ##get options from command line into variables and set default values
  36 my ($Verbose,$Help);
  37 GetOptions(
  38         "verbose"=>\$Verbose,
  39         "help"=>\$Help
  40 );
  41 die `pod2text $0` if (@ARGV == 0 || $Help);
  42
  43 my $Ref_num = 0;
  44 my $Qual_num = 0;
  45 my $Cyc_num = 0;
  46 my $Base_num = 0;
  47 my $head;
  48 my %data;
  49
  50 my %BaseNum = ("A",0,"C",1,"G",2,"T",3);
  51
  52 while (<>) {
  53         if (/^\#/) {
  54                 $head .= $_;
  55                 $Ref_num = $1 if (/Ref_base_number (\d+)/);
  56                 $Qual_num = $1 if (/Quality_number (\d+)/);
  57                 $Cyc_num = $1 if (/Cycle_number (\d+)/);
  58                 $Base_num = $1 if (/Seq_base_num (\d+)/);
  59         }
  60
  61         if (/^\w/) {
  62                 my @t = split /\s+/;
  63                 my $ref = shift @t;
  64                 my $cycle = shift @t;
  65                 pop @t;
  66                 $ref = $BaseNum{$ref};
  67                 $cycle -= 1;
  68
  69                 if (@t % $Qual_num != 0) {
  70                         print "error #####################\n";
  71                 }
  72                 for (my $i=0; $i<@t; $i++) {
  73                         my $qual = $i%$Qual_num;
  74                         my $base = ($i - $i%$Qual_num)/$Qual_num;
  75                         $data{$qual}{$cycle}{$ref}{$base} = $t[$i];
  76                 }
  77
  78         }
  79         if(/^</){
  80                 last;
  81         }
  82 }
  83
  84 print "#qual\tcycle\tA->A\tA->C\tA->G\tA->T\tC->A\tC->C\tC->G\tC->T\tG->A\tG->C\tG->G\tG->T\tT->A\tT->C\tT->G\tT->T\n";
  85
  86 #print Dumper \%data;
  87
  88 foreach my $qual (sort {$a<=>$b} keys %data) {
  89         my $qual_p = $data{$qual};
  90         foreach my $cyc (sort {$a<=>$b} keys %$qual_p) {
  91                 print "$qual\t$cyc";
  92                 my $cyc_p = $qual_p->{$cyc};
  93                 foreach my $ref (sort {$a<=>$b} keys %$cyc_p) {
  94                         my $ref_p = $cyc_p->{$ref};
  95                         my $total = 0;
  96                         foreach my $base (sort {$a<=>$b} keys %$ref_p) {
  97                                 my $freq = $ref_p->{$base};
  98                                 $freq = 1 if($freq == 0);
  99                                 $total += $freq;
 100                         }
 101                         foreach my $base (sort {$a<=>$b} keys %$ref_p) {
 102                                 my $freq = $ref_p->{$base};
 103                                 $freq = 1 if($freq == 0);
 104                                 printf("\t%e", $freq/$total);
 105
 106                         }
 107                 }
 108                 print "\n";
 109         }
 110 }
 111