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

rnaSpliceCounts.pl - Summarize counts of read that overlap known splice junctions.

=head1 SYNOPSIS

 rnaSpliceCounts.pl --chrom=chr1.fa <chr1.fa_splice_sorted.txt

Options:

=over 4

=item *

--chrom=s

Chromosome name

=item *

--totalReads=i

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

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

=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

Roman Petrovski

=cut

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

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

use lib '@CASAVA_FULL_PERL_LIBDIR@';
use Casava::PostAlignment::Sequencing::GenomicIOLib qw(parseSpliceName);

my $help = 0;
my $man = 0;
my $chrom;
my $totalReads;

GetOptions(
    "chrom|c=s"      => \$chrom,
    'totalReads=i'   => \$totalReads,
    "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: --totalReads not specified.\n") unless defined $totalReads;
pod2usage("$0: --totalReads is not valid.\n") unless ($totalReads>=0);

my $rpkmFactor = ($totalReads==0 ? 0. : (10**9/$totalReads));

my @junctions;

while(<>) {
    chomp;
    my @s = split("\t");
    my $info = { 'coverage' => $s[1] };
    parseSpliceName($s[0],%$info);
    push @junctions, $info;
}

for my $jref (sort {$a->{start} <=> $b->{start}} @junctions) {
    my $cov = $jref->{coverage};
    my $acov = ( $jref->{size} > 0 ? $cov/$jref->{size} : 0 );
    printf "%s\t%d\t%d\t%s\t%.5f\t%d\n",
      $chrom,$jref->{start},$jref->{end},$jref->{name},($acov*$rpkmFactor),$cov;
}

1;
__END__
