perl/etc/bam2chrs.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <sorted_bam> <outprefix>\n" if @ARGV != 2;
   6 my ($in,$out)=@ARGV;
   7
   8 open( HEAD,"-|","samtools view -H $in") or die "Error opening $in: $!\n";
   9 my ($LenT,@Heads,@ChrID,%ChrLen,%ChrFlag,%ChrStat);
  10 while (<HEAD>) {
  11         chomp;
  12         if (/\@SQ\tSN:(\S+)\tLN:(\d+)/) {
  13                 push @ChrID,$1;
  14                 $LenT += $2;
  15                 $ChrLen{$1} = $2;
  16                 $ChrFlag{$1} = 0;
  17                 $ChrStat{$1} = 0;
  18                 print ">$1: $2\n";
  19         } else {
  20                 push @Heads,$_;
  21         }
  22 }
  23 $ChrLen{'*'} = -$LenT;
  24 close HEAD;
  25
  26 open LOG,'>',"$out.log" or die "Error opening $out.log: $!\n";
  27 sub initOUT($) {
  28         my $outChr = $_[0];
  29         open( OUT,"|-","samtools view -bS - > $out.$outChr.bam") or die "Error opening $out.$outChr.bam: $!\n";
  30         print OUT "\@SQ\tSN:$outChr\tLN:$ChrLen{$outChr}\n",join("\n",@Heads),"\n";
  31 }
  32
  33 my $outChr = $ChrID[0];
  34 $ChrFlag{$outChr} = 1;
  35 initOUT($outChr);
  36
  37 open( IN,"-|","samtools view $in") or die "Error opening $in: $!\n";
  38 while (<IN>) {
  39         my $chr=(split /\t/)[2];
  40         ++$ChrStat{"\t"};
  41         if ($chr eq '*') {
  42                 ++$ChrStat{$chr};
  43                 next;   # unmap are skipped outputing
  44         }
  45         if ($chr eq $outChr) {
  46                 print OUT $_;
  47         } else {
  48                 $ChrFlag{$outChr} = 1;
  49                 close OUT;
  50                 die "File not sorted:[$chr] was before [$outChr].\n" if $ChrFlag{$chr};
  51                 initOUT($chr);
  52                 $outChr = $chr;
  53                 print OUT $_;
  54         }
  55         ++$ChrStat{$chr};
  56 }
  57 close IN;
  58 close OUT;
  59 print LOG "ChrID\tLen\tLenRatio\tCount\tRatio\tReletiveRatio\n";
  60 for my $chr (@ChrID,'*') {
  61         print LOG join("\t",$chr,$ChrLen{$chr},$ChrLen{$chr}/$LenT,$ChrStat{$chr},$ChrStat{$chr}/$ChrStat{"\t"},($ChrStat{$chr}/$ChrStat{"\t"})/($ChrLen{$chr}/$LenT)),"\n";
  62 }
  63 close LOG;