#!/usr/bin/perl
## =========================== NGS tools ==========================================
## NGS tools for metagenomic sequence analysis
## May also be used for other type NGS data analysis
##
##                                      Weizhong Li, UCSD
##                                      liwz@sdsc.edu
## http://weizhongli-lab.org/
## ================================================================================

use Getopt::Std;
getopts("T:o:",\%opts);
die usage() unless ($opts{T});

my $map_score_cutoff = $opts{T};
my $output           = $opts{o};

my ($i, $j, $k, $ll, $cmd);

my $fh = "STDOUT";
if ($output) {
  open(OUT, ">$output") || die "Can not write to $OUT";
  $fh = "OUT";
}

my $last_id = "";
my $R1_pri_aln = "";
my $R2_pri_aln = "";

while($ll=<>){
  if ($ll =~ /^\@/) { #### headers
    print $fh $ll;
    next;
  }
  else { #### alignment
    chop($ll);
    my @lls = split(/\t/,$ll);
    my $id  = $lls[0];
    my $rid = $lls[2];    if ($rid eq "*") {  next; }

    my $score;
    for ($j=11; $j<=$#lls; $j++) { if ($lls[$j] =~ /AS:i:(.+)/) { $score = $1; last; } }
    next unless (defined($score));
    next unless ($score >= $map_score_cutoff);

    my $flag_info = <<EOD;
        Chr     Flag    Description
        p       0x0001  the read is paired in sequencing
        P       0x0002  the read is mapped in a proper pair
        u       0x0004  the query sequence itself is unmapped
        U       0x0008  the mate is unmapped
        r       0x0010  strand of the query (1 for reverse)
        R       0x0020  strand of the mate
        1       0x0040  the read is the first read in a pair
        2       0x0080  the read is the second read in a pair
        s       0x0100  the alignment is not primary
        f       0x0200  QC failure
        d       0x0400  optical or PCR duplicate
EOD

    my $FLAG = $lls[1];
    next unless ($FLAG & 0x0001 );
    next unless ($FLAG & 0x0002 );
    next if     ($FLAG & 0x0004 );
    next if     ($FLAG & 0x0008 );

    if    ($id    ne $last_id) {
      my $print_flag = 1;
      if ($R1_pri_aln and $R2_pri_aln) {
        my @t_R1 = split(/\t/, $R1_pri_aln);
        my @t_R2 = split(/\t/, $R2_pri_aln);
        $print_flag = 0 unless ($t_R1[2] eq $t_R2[2]);
        $print_flag = 0 unless ($t_R1[6] eq "=");
        $print_flag = 0 unless ($t_R2[6] eq "=");
      }
      else {
        $print_flag = 0;
      }
      if ($print_flag) {
        print $fh "$R1_pri_aln\n$R2_pri_aln\n";
      }
      $R1_pri_aln = "";
      $R2_pri_aln = "";
    }

    if    ( $FLAG & 0x0040 ) { $R1_pri_aln = $ll unless ($R1_pri_aln); }
    elsif ( $FLAG & 0x0080 ) { $R2_pri_aln = $ll unless ($R2_pri_aln); }
    $last_id = $id;

  } #### alignment section
}

   #### last reads
    if    ($id    ne $last_id) {
      my $print_flag = 1;
      if ($R1_pri_aln and $R2_pri_aln) {
        my @t_R1 = split(/\t/, $R1_pri_aln);
        my @t_R2 = split(/\t/, $R2_pri_aln);
        $print_flag = 0 unless ($t_R1[2] eq $t_R2[2]);
        $print_flag = 0 unless ($t_R1[6] eq "=");
        $print_flag = 0 unless ($t_R2[6] eq "=");
      }
      else {
        $print_flag = 0;
      }
      if ($print_flag) {
        print $fh "$R1_pri_aln\n$R2_pri_aln\n";
      }
      $R1_pri_aln = "";
      $R2_pri_aln = "";
    }

close(OUT) if ($output);

sub usage {
<<EOD
Given a sam file generated by 
  bwa mem reference R1.fa_or_fq R2.fa_or_fq 
This script takes the sam from STDIN and prints out paired mapped reads pass alignment score cutoff
usage:
  $script_name  -T alignment_cutoff_score

  options
    -o output file, default STDOUT
    -T alignment cutoff score, defined in sam as AS:i:score

EOD
}
########## END usage

