bamReads2

#!/usr/bin/perl
use Getopt::long;
use strict;

my ($bam, $fq, $NMleft, $NMright, $multi);
GetOptions("b=s"=>\$bam, "f=s"=>\$fq, "l=s"=>\$NMleft, "r=s"=>\$NMright, "m=s"=>\$multi);

chomp(my $NML = `gunzip -c $NMleft | wc -l`);
chomp(my $NMR = `gunzip -c $NMright | wc -l`);
chomp(my $MM = `gunzip -c $multi | wc -l`);


system "samtools view $file | cut -f1,2,3,10 > bamreads_temp_$$";
my (%unique, %aligns);
open IN, "bamreads_temp_$$";
my ($idx, $flags, $chr, $seq);
while (<IN>) {
    ($idx, $flags, $chr, $seq) = split /\t/, $_;
    $unique{TI}{$idx} = 1;
    $unique{TS}{$seq} = 1;
    next if $chr eq '*';
    $unique{AI}{$idx} = 1;
    $unique{AS}{$seq} = 1;
    if ($flags & 128) {   # second read in pair
	$aligns{R}++;
    } else {
	$aligns{L}++;
    }
    $aligns{P}++ unless $flags & 8;
}
my $readlen = length($seq)-1;  # -1 due to newline
close IN;
system "rm -f bamreads2_temp_$$*";
my $TIunique = scalar keys %{ $unique{TI} };
my $TSunique = scalar keys %{ $unique{TS} };
my $AIunique = scalar keys %{ $unique{AI} };
my $ASunique = scalar keys %{ $unique{AS} };
$aligns{P} /= 2;   # these get counted twice
print "ReadLen\tE1.Alns\tE1.NM\tE1.Multi\tE2.Alns\tE2.NM\tE2.Multi\tPE.Alns\tTot.Unq.Rd\tAln.Unq.Rd\tTot.Unq.Seq\tAln.Unq.Seq\n";
print "$file\t$readlen\t$aligns{L}\t$aligns{R}\t$aligns{P}\t$TIunique\t$AIunique\t$TSunique\t$ASunique\n";
exit;