#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use Pod::Usage;


# Options definitions (name,format,default,required)
my @Options = ( ['help',undef,undef,0],
		['man',undef,undef,0],
		['MAPQ','i',undef,0],
	      );
# defaults values
my %Values=();
my %Required=();
map {$Values{$_->[0]}=$_->[2];$Required{$_->[0]}=$_->[3];} @Options;
# build options list
my %Options=();
map {$Options{defined($_->[1])?$_->[0].'='.$_->[1]:$_->[0]}=\$Values{$_->[0]};} @Options;
# retrieve options
GetOptions(%Options) || pod2usage(2);
pod2usage(1) if ($Values{'help'});
pod2usage(-exitstatus => 0, -verbose => 2) if ($Values{'man'});

# check required
map {pod2usage("$_ is required") if ($Required{$_}&&!defined($Values{$_}));} keys %Values;

# check input file
$#ARGV==0 || die "no input file given as parameter";

my $FH= $ARGV[0] eq '-'?'STDIN':'FIN';
if ($FH eq 'FIN') {
  open(FIN,$ARGV[0])||die "Can't open file $ARGV[0]";
}

while (my $line=<$FH>) {
  if ($line=~/^@/ || $line!~/XA:Z:/) { # header line or alignment line without XA tag
    print $line;
    next;
  }
  # build TAG index
  chomp $line;
  my @F=split/\t/,$line;
  my %TAG_index=();
  for (my $i=11;$i<=$#F;$i++) {
    $F[$i]=~/^(..)/;
    $TAG_index{$1}=$i;
  }

  # print main alignment
  $F[1]= $F[1] & 2041; # remove read mapped in proper pair and read unmapped flags
  printf "%s\n",join("\t",@F[0..$TAG_index{XA}-1],@F[$TAG_index{XA}+1..$#F]);
  $F[$TAG_index{XA}]=~s/^XA:Z://;

  my $MAPQ = $Values{MAPQ}?$Values{MAPQ}:$F[4];

  # print secondary alignments
  while ($F[$TAG_index{XA}]=~/([^,;]+),([-+]\d+),([^,]+),(\d+);/g) {
    my ($seq,$qual)=@F[9..10];
    if (($F[1]&0x10 && $2>0)||(!($F[1]&0x10) && $2<0)) { # secondary and main alignments have oposite strand
      $seq =~ tr/atgcrymkswATGCRYMKSW/tacgyrkmswTACGYRKMSW/;
      $seq = join('',reverse(split //,$seq));
      $qual = join('',reverse(split //,$qual));
    }
    printf "%s\n",join("\t",$F[0], 
		       ($F[1]& 1769)|($2<0?0x10:0), # remove previous tags plus read reversed and secondary alignment, don't add secondary alignment (bowtie don't do that)
		       $1, abs($2), $MAPQ, $3, @F[6..7], 0, $seq, $qual, "NM:i:$4");
  }
}

close (FIN) if ($FH eq 'FIN');


#============================================================

=head1 NAME

samNoXA.pl - duplicates alignment lines with XA tags

=head1 SYNOPSIS

samNoXA.pl [options] <SAM file>

=head1 OPTIONS

=over 8

=item B<-help>

Print a brief help message and exits.

=item B<-man>

Prints the manual page and exits.

=back

=item B<-MAPQ val>

Set the MAPing Quality for secondary alignments (default same as main alignment MAPQ)

=head1 INPUTS

=over 8

=item B<SAM file>

SAM file to be computed, or - for STDIN

=back

=head1 DESCRIPTION

For each alignment line containing an XA tag, several equivalent lines are generated without XA tag.

=head1 AUTHORS

Patrice DEHAIS

=head1 CONTACT

Questions can be posted to the sigenae mailing list:
sigenaesupport@jouy.inra.fr

=cut
