perl/etc/justonce/tmp/cmp2bam.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4 #use IO::Unread;
   5
   6 die "Usage: $0 <bigger_bam> <sorted_bam> <out>\n" if @ARGV < 3;
   7 my ($bam1,$bam2,$out)=@ARGV;
   8
   9 open IN1,"-|","samtools view $bam1" or die "Error opening $bam1: $!\n";
  10 open IN2,"-|","samtools view $bam2" or die "Error opening $bam2: $!\n";
  11
  12 my %Dat;
  13 while (<IN1>) {
  14         chomp;
  15         my ($id,$flag,$chr,$pos) = split /\t/;
  16         #$Dat{$id} = [$chr,$pos];
  17         $Dat{$id} = 1;
  18 }
  19 close IN1;
  20
  21 my ($lastCP,@CurrArray)=('');
  22 open OUT,'>',$out or die "Error opening $out: $!\n";
  23 while (<IN2>) {
  24         #chomp;
  25         my ($id,$flag,$chr,$pos) = split /\t/;
  26         my $t = "$chr,$pos";
  27         if ($lastCP eq $t) {
  28                 push @CurrArray,[$id,$_];
  29         } else {
  30                 $lastCP = $t;
  31                 my $flag=0;
  32                 for my $t (@CurrArray) {
  33                         unless (exists $Dat{$t->[0]}) {
  34                                 $flag=1;
  35                                 $t->[1] = '*'.$t->[1];
  36                                 #print $t->[0],"\n";
  37                         }
  38                 }
  39                 if ($flag==1) {
  40                         for my $t (@CurrArray) {
  41                                 print OUT $t->[1];
  42                         }
  43                         print OUT "\n";
  44                 }
  45                 @CurrArray=([$id,$_]);
  46         }
  47 }
  48 close IN2;
  49 close OUT;
  50
  51 __END__
  52 ./cmp.pl Tiger_aln_bam/pti096_clean_aln_pe.bam Tiger_aln_rmdup/pti096_clean_aln_pe_rmdup.bam pti096.cmp.sam &
  53
  54 # http://www.nature.com/nature/journal/v399/n6737/full/399682a0.html
  55 # Cultures in chimpanzees