#!/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

rnaExonCounts.pl - Summarize read bases counts for non-overlapping exons.

=head1 SYNOPSIS

 rnaExonCounts.pl --chrom=<chrom name> --QVCutOff=<number> --featureFile=<path> [--geneCountsFile=<path>] <sort.count files>

Options:

=over 4

=item *

--chrom=s

Chromosome name

=item *

--QVCutOff=i

Reads with single read alignment score less than QVCutOff are skipped

=item *

--featureFile=s

Path to exon_coords.txt file

=item *

--totalBases=i

Total number of mapped bases in the experiment. This value is used for RPKM normalization.

=item *

--geneCountsFile=s

Path to the file in which the gene coverage counts are stored.

=back

=head1 DESCRIPTION

Prints the splice site coverage report in the following format:

chromosomeName\tfeatureStartPos\tfeatureEndPos\tfeatureName\tRPKM\trawValue\n

=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
Casava::PostAlignment::Plugins::GenomicIOLib
Casava::PostAlignment::Sequencing::RnaSeqLib

=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

Chris Saunders

=cut

use warnings FATAL => 'all';
use strict;

use Pod::Usage;
use Getopt::Long;

use lib '@CASAVA_FULL_PERL_LIBDIR@';
use Casava::Common::Log qw(errorExit);
use Casava::PostAlignment::Sequencing::GenomicIOLib qw(printFeatureCounts readFeaturesUCSCB);
use Casava::PostAlignment::Sequencing::RnaSeqLib qw(getGenesFromExons);

my $help = 0;
my $man = 0;
my $chrom;
my $featureFile;
my $geneCountsFile;
my $totalBases = 0;

GetOptions(
    "chrom|c=s"         => \$chrom,
    'featureFile=s'     => \$featureFile,
    'geneCountsFile=s'  => \$geneCountsFile,
    'totalBases=i'      => \$totalBases,
    "help|?"            => \$help,
    'man'               => \$man) or pod2usage(2);

pod2usage(1) if $help;
pod2usage(-exitstatus => 0, -verbose => 2, -input => $1) if ($man and $0 =~ /(.*)/);
pod2usage("$0: --chrom not specified.\n") unless defined $chrom;
pod2usage("$0: --featureFile not specified.\n") unless defined $featureFile;

my $cmdGunzip = "gunzip";

my @sitesFiles = @ARGV;

my %exonStorage = ();

readFeaturesUCSCB( $chrom, $featureFile, %exonStorage );

# TODO -- pick these up from COLUMNS meta-data:
my $posCol = 1;
my $usedCol = 2;
my $filtCol = 3;

my $storageFeature = $exonStorage{feature};
my $pos = 0;
my $readBases = 0;

my $SITESFH = undef;

sub getNextLine() {
    while(1) {
        if(not defined $SITESFH) {
            last if(not scalar(@sitesFiles));
            my $sitesFile = shift @sitesFiles;
            if(-e $sitesFile) {
                open($SITESFH,"$cmdGunzip -c $sitesFile |")
                  or errorExit("ERROR: Can't open file: $sitesFile\n");
            }
        } else {
            my $line=<$SITESFH>;
            return $line if(defined $line);
            close($SITESFH);
            $SITESFH = undef;
        }
    }
    return undef;
}


for my $exon (sort {$storageFeature->{$a}->{end} <=> $storageFeature->{$b}->{end}} keys(%$storageFeature) )
{
    next if ($storageFeature->{$exon}->{end} < $pos);

    if ( $storageFeature->{$exon}->{start} <= $pos
      && $storageFeature->{$exon}->{end} >= $pos )
    {
        $storageFeature->{$exon}->{coverage} += $readBases;
    }

    while (1)
    {
        my $line=getNextLine();
        last if(not defined $line);
        next if ( $line =~ /^#/ );
        chomp($line);
        my @spl = split( '\t', $line );
        my ($new_pos,$used,$filt) = (int($spl[$posCol]),int($spl[$usedCol]),int($spl[$filtCol]));
        next if($new_pos < 1);

        errorExit("ERROR: unexpected line format at: $. in line: $line\n") if($new_pos <= $pos);

        $pos = $new_pos;
        $readBases = ($used+$filt);

        next if($storageFeature->{$exon}->{start} > $pos);
        last if($storageFeature->{$exon}->{end} < $pos);
        $storageFeature->{$exon}->{coverage} += $readBases;
    }
}

printFeatureCounts( *STDOUT, $chrom, %exonStorage, 0, $totalBases);

if ($geneCountsFile)
{
    my $featurefile = IO::File->new( ">$geneCountsFile" )
      || errorExit "ERROR: Couldn't create/open file handle for $geneCountsFile $!";
    my %geneStorage = ();
    getGenesFromExons( %exonStorage, %geneStorage );
    printFeatureCounts( $featurefile, $chrom, %geneStorage, 0, $totalBases);
    close $featurefile;
}

1;
__END__
