#!/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 $last_aln = "";
my $last_filter_flag = 0;

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

    my $score;
    for ($j=11; $j<=$#lls; $j++) { if ($lls[$j] =~ /AS:i:(.+)/) { $score = $1; last; } }
    if (defined($score)) {
      $this_flag = 1 if ($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 );

    if    ($id    ne $last_id) {
      if ($last_aln) {
        print $last_aln unless ($last_filter_flag);
      }
      $last_aln = "";
      $last_filter_flag = 0;
    }
    $last_aln .= "$ll\n";
    $last_id = $id;
    $last_filter_flag = $this_flag if ($this_flag);

  } #### alignment section
}

    #### last reads
    if    ($id    ne $last_id) {
      if ($last_aln) {
        print $last_aln unless ($last_filter_flag);
      }
      $last_aln = "";
      $last_filter_flag = 0;
    }
close(OUT) if ($output);

sub usage {
<<EOD
Given a sam file generated by 
  bwa mem host(human) R1.fa_or_fq R2.fa_or_fq 
This script takes the sam from STDIN, removes alignment where either R1 or R2 aligned with
reference, print other mappings in to sam

usage:
  $script_name  -T alignment_cutoff_score

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

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

