projects/radseq/novogene/splitbam.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <in bam> <out path>\n" if @ARGV<2;
   6
   7 my $in=shift;
   8 my $outp=shift;
   9 my $infilename = (split /\//,$in)[-1];
  10 $infilename =~ s/\.bam$//i;
  11
  12 my $SAMTOOLS='samtools';
  13 system($SAMTOOLS) == -1 and $SAMTOOLS='/usr/local/bin/samtools';
  14
  15 my %Chr2ID;
  16 open H, '-|', "$SAMTOOLS view -H $in" or die "Error opening $in: $!\n";
  17 my @HEADS = <H>;
  18 close H;
  19 for (@HEADS) {
  20         my ($type,$chr,$len)=split /\t/;
  21         next unless $type eq '@SQ';
  22         $chr=(split(':',$chr))[1];
  23         my $str = $chr;
  24         if ($str =~ /gb\|([\w.]+)\|/) {
  25                 $str=$1;
  26                 $str=~s/\.\d+$//;
  27         } else {
  28                 $str =~ s/\W/-/g;
  29                 $str=~s/-+$//;
  30         }
  31         print "$str, $chr, $len";
  32         die if exists $Chr2ID{$chr};
  33         $Chr2ID{$chr} = $str;
  34 }
  35
  36 my %Chr2FH;
  37 for my $k (keys %Chr2ID) {
  38         my $out = "$outp/$infilename.$Chr2ID{$k}.bam";
  39         open $Chr2FH{$k},'|-', "$SAMTOOLS view -bS - >$out" or die "Error opening $out: $!\n";
  40         print { $Chr2FH{$k} } @HEADS;
  41 }
  42 open I, '-|', "$SAMTOOLS view $in" or die "Error opening $in: $!\n";
  43 while (<I>) {
  44         my @items = split /\t/;
  45         next if $items[2] eq '*';
  46         my $fh = $Chr2FH{$items[2]} or die;
  47         print $fh $_;
  48 }
  49 close I;
  50 close $Chr2FH{$_} for keys %Chr2FH;
  51 __END__
  52