#!/usr/bin/env perl

=head1 NAME

alicat.pl

=head1 SYNOPSIS

alicat.pl [-a] [-c chr] [-g] [-l loOffset] [-h hiOffset] [-i pats] [-o file] [-p [-q qTh]] [-s] [-t len] [-v] [-x pats] [-z] file1 file2 ...

=head1 DESCRIPTION

Alignment file equivalent to "cat -n". The input files are as produced
by the SGRP pipeline, e.g. layout.gz, sorted.gz, imputed.gz, under the
"match" part of the hierarchy. 

If -l and/or -h are specified, only positions within that region are
output.

If -a is given, offsets refer to positions in the alignment,
i.e. including reference-sequence gaps; otherwise (the default) gaps
do not add to the offset and we show such lines with offsets NNN.001,
NNN.002 etc.

If -z is given, start from offset 0, else offset 1.

If -o is given, output to that file; default is standard output. If file name
ends in .gz, gzip to it.

If -p is given, only polymorphic data lines (those with more than one
different nucleotide) are printed.

If -q is given (as well as -p), only count values with at least that quality
(as a number) when deciding polymorphism.

If -s is given, separate nucleotides and qualities in output (see below).

If -g is given, leave numeric gap characters (created during SGRP
imputation process) as they are; default is to translate them to '-'.

If -i str1,str2,... is given, include only strains matching one of those names
    (Perl-style regular expressions allowed) except any matching -x list.

If -v or -t is given, print strain names vertically at the start.

If -t is given, turn the whole alignment to be horizontal, using line
length given by -t.  For this to work, you will need to have the
script "aliturn.pl" in the same directory as this one.

If -x str1,str2,... is given, exclude strains matching any of those names.

Each input line starting with ">" is printed as is, except that for
the first line, if -v is given, we print the strain names on it vertically.
Otherwise, the value of "-c" (if any) is put in after the
">". Thus if we say "alicat.pl -c foo bar", and the first ">" line is
">xyz", the output will be ">foo xyz".

Other input lines are taken to be alternating nucleotide and quality
characters. If "-l" and/or "-h" are given, only lines in that range of
offsets (counting from 1, or from 0 if "-z" is given) are printed, and
if "-p" is given, only polymorphic lines are printed.  For such
inputs, the columns in the output are as follows:

1-7:   the offset in the chromosome
8-11:  if the first character in the input line is a gap symbol
       ("-" or a digit), cols 8-11 are of the form ".NNN", and
        cols 1-7 are the same as in the previous line; otherwise they
        are blank. Thus MMM.NNN means "the NNN'th inserted position
        after offset MMM".
13-15: Consensus nucleotide value. If consensus is unanimously X,
       field is " X "; if not unanimous, i.e. if the line is polymorphic, 
       it's "[X]". A lower-case consensus means there is a gap in some 
       positions and nucleotides in others.
16:    always blank
17:    always "|"

For the remaining columns: N is the number of strains including
the reference, i.e. each non-">" input line has 2*N characters:

If "-s" ("separate") switch is given, then remaining columns are:

18 to 17+N: the nucleotide values from the input (with digits translated
            to "-" if "-t" was given).

18+N to 20+N: always "| |"
21+N to 20+2*N: the quality characters from the input.
21+2*N:       always "|"

If "-s" is not given, then:

18 to 17+2*N: an exact copy of the input, i.e. with nucleotides and qualities
              still interleaved (but digits are translated to "-" if "-t" was 
              given).
18+2*N:       always "|"
       
=head1 AUTHOR

David Carter, C<dmc@sanger.ac.uk>

=head1 TERMS AND CONDITIONS

Copyright (c) 2007 Wellcome Trust Sanger Institute (WTSI)

This program is free software; you can redistribute it and/or modify
it under the same terms as Perl itself. See the Artistic License file
in the main Perl distribution for specific terms and conditions of
use.  In addition, the following disclaimers apply:

WTSI makes no representations whatsoever as to the SOFTWARE contained
herein.  It is experimental in nature and is provided WITHOUT WARRANTY
OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE OR ANY OTHER
WARRANTY, EXPRESS OR IMPLIED.  CSHL MAKES NO REPRESENTATION OR
WARRANTY THAT THE USE OF THIS SOFTWARE WILL NOT INFRINGE ANY PATENT OR
OTHER PROPRIETARY RIGHT.

By downloading this SOFTWARE, your Institution hereby indemnifies WTSI
against any loss, claim, damage or liability, of whatsoever kind or
nature, which may arise from your Institution's respective use,
handling or storage of the SOFTWARE.

If publications result from research using this SOFTWARE, we ask that
WTSI be acknowledged and/or credit be given to WTSI scientists, as
scientifically appropriate.

=head1 SEE ALSO

       

=cut

use strict;

use File::Basename;

use FileHandle;
use Getopt::Std;
use strict;

my %OPT;

getopts('ac:gi:l:h:o:pq:st:vx:z',\%OPT);

my $PERL_DIR = dirname($0);

# Bit positions are:
#  A 1
#  C 2
#  G 4
#  T 8
#  - 16

my @DNA_SYM;
my %DNA_BITS; 

setDNABitsAndSymbols();
my $gh = openToWrite($OPT{o},$OPT{t});
foreach my $f (@ARGV) {
    alicat($f,$gh,$OPT{c},$OPT{l},$OPT{h},$OPT{p},$OPT{q},$OPT{s},!$OPT{g},$OPT{i},$OPT{x},$OPT{v} || $OPT{t},$OPT{t});
}

sub setDNABitsAndSymbols {
    %DNA_BITS = qw(A  1 C  2 G  4 T  8
		   M  3 R  5 W  9 
		   S  6 Y 10 K 12
		   B 14 D 13 H 11 V  7
		   N 15 - 16 = 16
		   a 17 c 18 g 20 t 24
		   m 19 r 21 w 25 
		   s 22 y 26 k 28
		   b 30 d 29 h 27 v 23
		   n 31);
    @DNA_SYM = ();
    while (my ($k,$v) = each %DNA_BITS) {
	# Make sure it's been "used as a number" so 
	# bitwise "|" works numerically!
	$DNA_BITS{$k} += 0;
	$DNA_SYM[$v] = $k;
    }
}	      

sub alicat {
    my ($f,$gh,$chr,$lo,$hi,$polyOnly,$qTh,$sep,$trans,$incPats,$excPats,$vertHdr,$turn) = @_;
    # Offset:
    my $n = $OPT{z} ? -1 : 0;
    # Delta-offset (increment for each nucleotide, unless -a)
    my $d=0;
    my $fh = openToRead($f);
    chomp (my $hdr = <$fh>);
    my @strains = ('REF',split(" ",substr($hdr,1)));
    my @accept = applyIncAndExc($incPats,$excPats,\@strains);
    if ($vertHdr) {
	printHeader(\@strains,\@accept,$sep,$gh);
    } else {
	print $gh ">";
	if ($chr) {
	    print $gh "$chr ";
	}
	my $sep = '';
	for (my $k=0; $k<=$#strains; $k++) {
	    if ($accept[$k]) {
		print $gh "$sep$strains[$k]";
		$sep = ' ';
	    }
	}
	print $gh "\n";
    }
    my $incrementOffsetAtGaps = $OPT{a};
    my $expectedLineLength = scalar(@strains)*2;
    while (<$fh>) {
	chomp;
	my $dl = $expectedLineLength-length($_);
	$_ .= ' ' x $dl if ($dl > 0);
	if (/^>/) {
	    if ($chr) {
		# Smuggle in dir name
		printf($gh ">%s %s\n",$chr,substr($_,1));
		$chr = undef; # stop it happening again
	    } elsif (/^>.*\S/) {
		# Print if we have something other than a space after the ">"
		print $gh "$_\n";
	    }
	} else {
	    $_ = pruneDataLine(\@accept,$_);
	    if (defined $lo && $n+1 < $lo) {
		$n++ unless (/^[-_0-9]/);
		next;
	    }
	    last if (defined $hi && $n+1 > $hi);
	    my $cons = consensusValue($_,$qTh);
	    next if ($cons eq '-'); # can happen if we have exclusions
	    my $kMax = length($_)-2;
	    if ($trans) {
		for (my $k=0; $k<=$kMax; $k+= 2) {
		    substr($_,$k,1) =~ s/[0-9]/-/;
		}
	    }
	    if ($sep) {
		my ($nucs,$quals);
		for (my $k=0; $k<=$kMax; $k+= 2) {
		    $nucs .= substr($_,$k,1);
		    $quals .= substr($_,$k+1,1);
		}
		$_ = "$nucs| |$quals";
	    }		
	    if (/^[-=_0-9]/) {
		unless ($cons eq ' - ') { # _all_ gaps!
		    if ($incrementOffsetAtGaps) {
			$n++;
		    } else {
			$d++;
		    }
		    my $dd = sprintf("%3d",$d);
		    $dd =~ tr/ /0/;
		    printf($gh "%7d.%s %s |%s|\n",$n,$dd,$cons,$_);
		}
	    } else {
		$d=0;
		$n++;
		unless ($polyOnly && substr($cons,0,1) eq ' ') {
		    printf($gh "%7d     %s |%s|\n",$n,$cons,$_);
		}
	    }
	}
    }
}

sub applyIncAndExc {
    my ($incPats,$excPats,$pStrains) = @_;
    my @accept;
    my @incPats = split(',',$incPats);
    my @excPats = split(',',$excPats);
    foreach my $str (@$pStrains) {
	my $got = 0;
	if (@incPats && $str ne 'REF') {
	    foreach my $pat (@incPats) {
		if ($str eq $pat || $str =~ /^$pat$/) {
		    $got = 1;
		    last;
		}
	    }
	} else {
	    $got = 1;
	}
	if (@excPats && $str ne 'REF') {
	    foreach my $pat (@excPats) {
		if ($str eq $pat || $str =~ /^$pat$/) {
		    $got = 0;
		    last;
		}
	    }
	}
	push @accept,$got;
    }
    return @accept;
}

sub pruneDataLine {
    my ($pAccept,$line) = @_;
    my $ans;
    foreach my $i (0 .. $#$pAccept) {
	if ($pAccept->[$i]) {
	    $ans .= substr($line,$i*2,2);
	}
    }
    return $ans;
}

sub printHeader {
    my ($pStrains,$pAccept,$sep,$gh) = @_;
    my $k=0;
    my $line;
    my $going = 1;
    my $cons = 'Cons';
    while ($going) {
	$going = 0;
	$line = sprintf("%13s%1s%3s","",$k >= 4 ? ' ' : substr($cons,$k,1),"");
	$line .= printStrainRow($k,$pStrains,$pAccept,$sep);
	if ($sep) {
	    $line .= '   ';
	    $line .= printStrainRow($k,$pStrains,$pAccept,$sep);
	}
	$line =~ s/  +$//;
	print $gh "$line\n";
	if ($line) {
	    $going = 1;
	}
	$k++;
    }
}

sub printStrainRow {
    my ($k,$pStrains,$pAccept,$sep) = @_;
    my $line;
    my @acc = @$pAccept;
    foreach my $str (@$pStrains) {
	next unless (shift @acc);
	if (length($str) > $k) {
	    $line .= substr($str,$k,1);
	} else {
	    $line .= ' ';
	}
	$line .= ' ' unless ($sep);
    }
    return $line;
}

sub consensusValue {
    my ($line,$qTh) = @_;
    my $kMax = length($line)-2;
    my $cons = substr($line,0,1);
    $cons =~ tr/_a-z/-A-Z/;
    $cons =~ s/[0-9]/-/g;		    
    my $cons1 = $cons;
    my $qc;
    $qc = chr($qTh+33) if (defined $qTh);
    for (my $k=2; $k<=$kMax; $k+=2) {
	if (defined $qc) {
	    my $qk = substr($line,$k+1,1);
	    next if ($qk ne ' ' && $qk lt $qc);
	}
	$cons = consensusValue1($cons,substr($line,$k,1));
	#printf("'%s' '%s'\n",substr($line,$k,1),$cons);
    }
    return $cons1 eq $cons ? " $cons " : "[$cons]";
}

sub consensusValue1 {
    my ($v1,$v2) = @_;
    if ($v2 eq ' ') {
	return $v1;
    } else {
	$v2 =~ tr/_a-z/-A-Z/;
	$v2 =~ s/[0-9]/-/g;	   
	return $DNA_SYM[$DNA_BITS{$v1} | $DNA_BITS{$v2}];
    }
}

sub openToRead {
    my ($f) = @_;
    if (!$f) {
	return *STDIN;
    } elsif ($f =~ /\.gz$/) {
	return new FileHandle("gunzip -c $f |");
    } else {
	return new FileHandle($f,'r');
    }
}
    
sub openToWrite {
    my ($g,$turn) = @_;
    my $turner = "$PERL_DIR/aliturn.pl $turn";
    if (!$g) {
	if ($turn) {
	    return new FileHandle("| $turner");
	} else {
	    return *STDOUT;
	}
    } elsif ($g =~ /\.gz$/) {
	return new FileHandle(($turn ? "| $turner " : 0) . "| gzip -c > $g");
    } elsif ($turn) {
	return new FileHandle("| $turner > $g");
    } else {
	return new FileHandle($g,'w');
    }
}
    
