#!/usr/bin/env perl

=head1 LICENSE

Copyright (c) 2007-2009 Illumina, Inc.

This software is covered by the "Illumina Genome Analyzer Software
License Agreement" and the "Illumina Source Code License Agreement",
and certain third party copyright/licenses, and any user of this
source file is bound by the terms therein (see accompanying files
Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
Illumina_Source_Code_License_Agreement.pdf and third party
copyright/license notices).

This file is part of the Consensus Assessment of Sequence And VAriation
(CASAVA) software package.

=head1 NAME

nonoverlappingExonCoords.pl - produce non-redundant exon set.

=head1 SYNOPSIS

 nonoverlappingExonCoords.pl [options] < /path/to/refFlat.txt >exon_coords.txt 

Options:

=over 4

=item *

--annotationChromPrefix=s

Prefix to be replaced by --referenceChromPrefix when constructing reference
chromosome names from annotation file chromosome names

=item *

--referenceChromPrefix=s

Prefix to be replace the --annotationChromPrefix when constructing reference
chromosome names from annotation file chromosome names

=item *

--chromosomeSuffix=s

If specified, the value is appended to annotation chromosome names when constructing
reference chromosome names. Default is '.fa' if --chromNameSource=fileName.

=item *

--chromNameSource=contigName|fileName

Affects only --chromosomeSuffix behavior.

=back

=head1 DESCRIPTION

TheNon-redundant exon set is used by CASAVA to compute the gene expression 
counts. For every gene, overlapping exons from the same gene are combined to
form one larger region used for counting. Regions where there are overlapping 
exons from other genes are excluded from counting and are not present in the 
non-redundant exon set file.

=head1 DIAGNOSTICS

=head2 Exit status

0: successful completion
1: abnormal completion
2: fatal error

=head2 Errors

All error messages are prefixed with "ERROR: ".

=head2 Warnings

All warning messages generated by CASAVA are prefixed with "WARNING: ".

=head1 CONFIGURATION AND ENVIRONMENT

=head1 DEPENDENCIES

=over 4

=item Standard perl modules

strict, warnings, Getopt::Long, Pod::Usage

=item External perl modules

=item CASAVA perl modules

Casava::Common::Log

=back

=head1 BUGS AND LIMITATIONS

There are no known bugs in this module.

All documented features are fully implemented.

Please report problems to Illumina Technical Support (support@illumina.com)

Patches are welcome.

=head1 AUTHOR

Jerry Kakol

=cut
use strict;
use warnings;
use Getopt::Long qw(:config gnu_compat);
use Pod::Usage;

use lib '/home/psgendb/local/pkg/CASAVA_v1.8.2-build/lib/CASAVA-1.8.2/perl';

use Casava::Common::Log qw(errorExit logInfo initLog);

# Input as per refFlat table is defined as follows:
#
# CREATE TABLE `refFlat` (
#   `geneName` varchar(255) NOT NULL default '',
#   `name` varchar(255) NOT NULL default '',
#   `chrom` varchar(255) NOT NULL default '',
#   `strand` char(1) NOT NULL default '',
#   `txStart` int(10) unsigned NOT NULL default '0',
#   `txEnd` int(10) unsigned NOT NULL default '0',
#   `cdsStart` int(10) unsigned NOT NULL default '0',
#   `cdsEnd` int(10) unsigned NOT NULL default '0',
#   `exonCount` int(10) unsigned NOT NULL default '0',
#   `exonStarts` longblob NOT NULL,
#   `exonEnds` longblob NOT NULL);
#
# Warning:
#    `exonStarts` and `exonEnds` fileds end with spurious commas.

my $annotationChromPrefix = '';
my $referenceChromPrefix = '';
my $chromosomeSuffix = '';
my $chromNameSource = '';
my $help = 0;
my $man = 0;
my $verbose = 4;

GetOptions('annotationChromPrefix=s'    => \$annotationChromPrefix,
           'referenceChromPrefix=s'     => \$referenceChromPrefix,
           'chromosomeSuffix=s'  => \$chromosomeSuffix,
           'chromNameSource=s'      => \$chromNameSource,
           'verbose|v=i'           => \$verbose,
           'help|?'           => \$help,
           'man'              => \$man) or pod2usage(2);

pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2, -input => $1) if ($man and $0 =~ /(.*)/);
pod2usage("$0: --chromNameSource must be contigName or fileName.\n") unless grep /^$chromNameSource$/, ('contigName', 'fileName');

$chromosomeSuffix = '.fa' if (!$chromosomeSuffix and 'fileName' eq $chromNameSource);

initLog( undef, $verbose, 0 ); #log everything to stderr until we get projectDir and verbose level

my %exons;

while (<>)
{
  chomp;
  
  my ($symbol, $accession, $chromosome, $strand, $tx_start, $tx_end, $cds_start, $cds_end, $exon_count, $exon_starts, $exon_ends) = split (/\t/);

  $exon_starts =~ s/,$//;
  $exon_ends   =~ s/,$//;

  my @exon_starts = split (/,/, $exon_starts);
  my @exon_ends   = split (/,/, $exon_ends);
  for (my $i = 0; $i < @exon_starts; $i++) { ++$exon_starts[$i]; }
  
  $symbol = $accession if (!$symbol);
  if ($symbol !~ /^[[:alnum:]]+$/)
  {
    $symbol =~ s/[^[:alnum:]]+/-/g;
  }
    
  $chromosome =~ s/^$annotationChromPrefix/$referenceChromPrefix/;
  my $l = @exon_starts;
  
  for (my $i = 0; $i < $l; $i++)
  {
    if (exists $exons{$chromosome}{$exon_starts[$i]})
    {
      my $prev_end = $exons{$chromosome}{$exon_starts[$i]}{"end"};

      if ($exons{$chromosome}{$exon_starts[$i]}{"symbol"} ne $symbol)
      {
        $exons{$chromosome}{$exon_starts[$i]}{"symbol"} = "MULTIPLE_GENES";
        if ($exon_ends[$i] > $prev_end)
        {
          %{$exons{$chromosome}{$prev_end + 1}} = ("symbol" =>$symbol, "end" => $exon_ends[$i]);
        }
      }
    }
    else
    {
      %{$exons{$chromosome}{$exon_starts[$i]}} = ("symbol" =>$symbol, "end" => $exon_ends[$i]);
    }
  }
}

logInfo "Finished reading and consolidating exon coordinates\nConsolidating partially overlapping exons...", 1;

foreach my $chr (keys %exons)
{
  my $k = 0;
  my $l = scalar (keys (%{$exons{$chr}})) - 1;

  while ($k < $l)
  {
    my @exon_starts = sort { $a <=> $b } keys (%{$exons{$chr}});
    $l = @exon_starts - 1;
 
    for (my $i = $k; $i < $l; $i++, $k++)
    {
      my $prev_start = $exon_starts[$i];
      my $prev_end = $exons{$chr}{$prev_start}{"end"};
      my $next_start = $exon_starts[$i + 1];
      my $next_end = $exons{$chr}{$next_start}{"end"};

      my $prev_symbol = $exons{$chr}{$prev_start}{"symbol"};
      my $next_symbol = $exons{$chr}{$next_start}{"symbol"};

      if ($prev_end > $next_start)
      {
        # This check can see only single overlap, that is prev_exon overlaping next_exon, 
        # while in fact it may be possible that the two are overlapped also by next_next_exon
        # and so on. That's why we fix only boundaries of the first two and restart the
        # scan with the next iteration of the 'while ($k < $l)' loop

        if ($next_end > $prev_end)
        {
          if ($prev_symbol ne $next_symbol)
          {
            # Next exon partially overlap previous (current $i) exon
            # Adjust prev_end and next_start accordingly in %final_exons
            if ($prev_start < $next_start - 1) { $exons{$chr}{$prev_start}{"end"} = $next_start - 1; }
            if (($prev_start + 1) < $next_end) 
            { 
              %{$exons{$chr}{$prev_end + 1}} = %{$exons{$chr}{$next_start}};
              $l++;
            }
          }
          else
          {
            $exons{$chr}{$prev_start}{"end"} = $next_end;
          }
          delete $exons{$chr}{$next_start};
          $l--;
        }
        else
        {
          # Next exon is fully contained in the current exon and needs to be skipped
          # while prev exon has to be divided in two non-overlapping parts
          if ($prev_symbol ne $next_symbol)
          {
            # Symbols are different, therefore non-overlapped part should be retained
            if ($prev_start < ($next_start - 1)) { $exons{$chr}{$prev_start}{"end"} = $next_start - 1; }
            if (($next_end + 1) < $prev_end) 
            { 
              %{$exons{$chr}{$next_end + 1}} = ("symbol" => $prev_symbol, "end" => $prev_end);
              $l++;
            }
            logInfo "Deleting exon $chr:$next_start-$next_end of gene $next_symbol as fully contained in exon $chr:$prev_start-$prev_end of gene $prev_symbol", 3;
          }
          else
          {
            logInfo "Deleting exon $chr:$next_start-$exons{$chr}{$next_start}{end} from gene $next_symbol as a smaller exon contained in a bigger one", 3;
          }
          delete $exons{$chr}{$next_start};
          $l--;
        }
        last;
      }
    }
  }

  logInfo "Finished chromosome $chr", 2;
}

logInfo "Finished gathering and consolidating splice sites coordinates\nExtracting sequences...", 1;

foreach my $chr (keys %exons)
{
  next if $chr =~ /_/;

  foreach my $exon_start (keys %{$exons{$chr}})
  {
    my $symbol = $exons{$chr}{$exon_start}{"symbol"};
    my $exon_end = $exons{$chr}{$exon_start}{"end"};

    next if $symbol eq "MULTIPLE_GENES";

    # A sanity check:
    errorExit "Exon start $exon_start is greater then exon end $exon_end for gene $symbol!" if $exon_start > $exon_end;

    print "${chr}${chromosomeSuffix}\t$exon_start\t$exon_end\t$symbol\n";
  }
}
