projects/recombineMap/sperm_by_fw/snp_calling/make_SE_pirs_matrix.pl

   1 #!/usr/bin/perl
   2
   3 =head1 Name
   4
   5 snp_dist_filter.pl   --  filter snp by neigboring distance
   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   snp_dist_filter.pl <*.cns.snp>
  19   --distance <int>  set the minimum distance between two neighboring SNPs, default=5
  20   --verbose   output running progress information to screen
  21   --help      output help information to screen
  22
  23 =head1 Exmple
  24
  25
  26
  27 =cut
  28
  29 use strict;
  30 use Getopt::Long;
  31 use FindBin qw($Bin $Script);
  32 use File::Basename qw(basename dirname);
  33 use Data::Dumper;
  34 use File::Path;  ## function " mkpath" and "rmtree" deal with directory
  35
  36 ##get options from command line into variables and set default values
  37 my ($Dist_cutoff, $Verbose,$Help);
  38 GetOptions(
  39         "distance:i"=>\$Dist_cutoff,
  40         "verbose"=>\$Verbose,
  41         "help"=>\$Help
  42 );
  43 $Dist_cutoff ||= 5;
  44 die `pod2text $0` if (@ARGV == 0 || $Help);
  45
  46 my $head;
  47 my %data;
  48
  49 while (<>) {
  50         if (/^\#/) {
  51                 $head .= $_;
  52         }
  53
  54         if (/^\w/) {
  55                 my @t = split /\s+/;
  56                 my $ref = shift @t;
  57                 my $cycle = shift @t;
  58                 $data{$ref}{$cycle} = \@t;
  59         }
  60         if(/^</){
  61                 last;
  62         }
  63 }
  64
  65 $head =~ s/Cycle_number 200/Cycle_number 100/;
  66
  67 print $head;
  68
  69 #print Dumper \%data;
  70
  71 my $outstr;
  72 foreach my $ref (sort keys %data) {
  73         my $ref_p = $data{$ref};
  74         for (my $i=1; $i<=100; $i++) {
  75                 $outstr .= "$ref\t$i";
  76                 my $read1_p = $ref_p->{$i};
  77                 my $read2_p = $ref_p->{$i+100};
  78                 for (my $j=0; $j<@$read1_p; $j++) {
  79                         my $sum = $read1_p->[$j] + $read2_p->[$j];
  80                         $outstr .= "\t$sum";
  81                 }
  82                 $outstr .= "\n";
  83         }
  84 }
  85
  86 print $outstr;
  87