#!/usr/bin/perl 

#===============================================================================
#
#         FILE:  get_input_for_LOX.pl
#
#        USAGE:  ./get_input_for_LOX.pl 
#
#  DESCRIPTION:  This is a utility to create the input file for LOX software.
#                To create the input file, the script uses the output files from 
#                the ELAND software. 
#                ELAND is a short-read aligner part of the Illumina/Solexa suite.
#                This script only considers the unique exact matches from Eland 
#                output.
#                After the name of the Perl script, the user needs to indicate 
#                a list of input files.
#                The first file should be the CDS of the genome (fasta format).
#                This is the file used by the Perl script to census 
#                the short reads.
#                The rest of files should be the ELAND output files.
#                Each ELAND output file should correspond to a different 
#                treatment.
#                Although the application only considers the unique exact 
#                matches from ELAND, it does not discard the possibility of 
#                multiple matching against CDSs.
#                Multiple matching can occur if ELAND matched the short reads 
#                against a genome sequence and CDS for that genome has 
#                overlapping CDSs (like yeast).
#                Reads with multiple matches will be excluded to produce the 
#                LOX input.
#
#      OPTIONS:  ---
# REQUIREMENTS:  For portability, the code is designed to expressly avoid 
#                the use of packages.
# 
#        NOTES:  Logging is done to STDOUT 
#       AUTHOR:   (Francesc Lopez-Giraldez), <francesc.lopez@yale.edu>
#      VERSION:  0.1
#      CREATED:  12/23/2009 07:33:56 PM EST
#===============================================================================

use strict;
use warnings;

my $USAGE = <<"END_USAGE";


Usage: After the name of the Perl script, the user needs to indicate 
       a list of input files.
       The first file should be the CDS of the genome (fasta format).
       This is the file used by the Perl script to census 
       the short reads.
       The rest of files should be the ELAND output files.
       Each ELAND output file should correspond to a different 
       treatment.

       EX:
       $0 some_yeast_CDS.txt s1_ELAND.txt s2_ELAND.txt s3_ELAND.txt s4_ELAND.txt

       or:
       perl $0 some_yeast_CDS.txt s1_ELAND.txt s2_ELAND.txt s3_ELAND.txt s4_ELAND.txt


END_USAGE

my $initial_time = time;
my $CDS_filename    = shift @ARGV or die $USAGE;
my @ELAND_filenames = @ARGV or die $USAGE;
my %census_results_for;

# Logging
print "\nInput CDS file    ::: $CDS_filename\n";
print "Input ELAND files ::: @ELAND_filenames\n";
#--

my %CDS_sequences_for = %{ get_CDS_sequences($CDS_filename) };

# Logging
my $total_CDS = scalar keys %CDS_sequences_for;
print "Total CDS ::: $total_CDS\n";
#--

foreach my $ELAND_file (@ELAND_filenames) {
    print "\tReading ELAND file [$ELAND_file]\n";
    censing_reads_from($ELAND_file);
}


# Printing LOX input file #####
print "\nPrinting LOX input file.....\n";
my $LOX_infile = 'LOXinputfile.txt';

open my $LOX, '>', $LOX_infile
    or die "$0: failed to open output file '$LOX_infile': $!\n";

printing_header_rows_for($LOX);

#Printing all CDSs including the ones with 0 hits in all treatments.
foreach my $CDS (keys %CDS_sequences_for) {
    print {$LOX} "\n$CDS\t#N/A";
    foreach my $treatment (@ELAND_filenames) {
        if (exists $census_results_for{$CDS}{$treatment}) {
            print {$LOX} "\t$census_results_for{$CDS}{$treatment}";
        }
        else {
            print {$LOX} "\t0";
        }
    }
}

close $LOX
    or warn "$0: failed to close output file '$LOX_infile': $!\n";



print "\nDONE!\n";
my $elapsed_time = time - $initial_time;
printf "Total time ::: %d days, %d hours, %d minutes and %d seconds.\n\n",
                (gmtime $elapsed_time)[7,2,1,0];

exit;



############################################################################
#---------------------------------------------------------------------------
#  SUBROUTINES
#---------------------------------------------------------------------------

# Read file and get sequences in a hash.
# Key=Sequence ID. The first word of the fasta file (usually the gene id). 
# Value=Sequence.
sub get_CDS_sequences {
    my ($CDS_filename) = @_;
    my $file_content;
    my %CDS_sequences_for;
    my $fasta_separator = q{>};

    open my $filehandle, '<', $CDS_filename
        or die "$0: failed to open input file '$CDS_filename': $!\n";
    
    LINE:
    while (my $line = <$filehandle>) {
        next LINE if $line !~ /\S/xm; #read next line if empty
        $file_content .= $line;
    }

    close $filehandle
        or warn "$0: failed to close input file '$CDS_filename': $!\n";

    my @fasta_sequences = split $fasta_separator, $file_content;
    #first entry is empty
    shift @fasta_sequences;

    FASTA:
    foreach my $fasta_sequence (@fasta_sequences) {
        my ($sequence_ID) = $fasta_sequence =~ /^\s*(\S+)/xm;
        #Remove fasta title line, keep only sequence
        $fasta_sequence =~ s/^([^\r\n]+)//xm;
        #Remove white spaces and new lines/return
        $fasta_sequence =~ s/\s+//xmg;
        #Check if sequence already exists        
        die "Sequence_ID [$sequence_ID] is duplicated in CDS file.\n"
                    if (exists $CDS_sequences_for{$sequence_ID});
        #Check sequence exists
        if (!$fasta_sequence) {
            warn "Sequence_ID [$sequence_ID] doesn't have a sequence"
                 . " in CDS file.\n"
                 ;
            next FASTA;
        }

       #Store sequence in hash       
       $CDS_sequences_for{$sequence_ID} = uc $fasta_sequence;
    }

    return \%CDS_sequences_for;
}   # ----------  end of subroutine get_CDS_sequences  ----------



# Open and reads content of ELAND output files
# ELAND single exact reads are matched against the CDSs 
sub censing_reads_from {
    my ( $ELAND_filename ) = @_;
    my $total_reads = 0;
    my %match_type_count_for = (
        0 => 0, #no_matches_count
        1 => 0, #single_matches_count
        2 => 0, #multiple_matches_count
    );
    my $total_U0_reads = 0; # unique exact matches in ELAND file
    local $|=1;

    open my $FH, '<', $ELAND_filename
        or die "$0: failed to open input file '$ELAND_filename': $!\n";

    LINE:
    while (my $line = <$FH>) {
        next LINE if $line !~ /^\s*>/xm; #read next line if no read
        ++$total_reads;
        
        my @fields = split /\t/xm, $line;
        # Third field is type of match.
        # U0—Best match found was a unique exact match
        next LINE if $fields[2] ne 'U0';
        ++$total_U0_reads;

        #The second field in ELAND output file is the sequence/read (see below)
        my $query_seq = uc $fields[1];
        my $reverse_query_seq = reverse $query_seq;
        (my $comp_rev_query_seq = $reverse_query_seq) =~ tr/ACTG/TGAC/;
       
        my $match_type = get_matches_for($query_seq, $comp_rev_query_seq, $ELAND_filename);
        ++$match_type_count_for{$match_type};

        if ($total_reads %100 == 0) {
            print "\r\t\t$total_reads total reads.";
        }
    }

    close $FH
        or warn "$0: failed to close input file '$ELAND_filename': $!\n";

    #LOGging
    print "\n\t\tTotal ELAND reads: $total_reads\n";
    print "\t\tTotal U0 reads:    $total_U0_reads\n";
    print "\t\t\tReads with no matches      : $match_type_count_for{0}\n";
    print "\t\t\tReads with single matches  : $match_type_count_for{1}\n";
    print "\t\t\tReads with multiple matches: $match_type_count_for{2}\n";
    my $unmatched_reads = $total_reads - $match_type_count_for{0}
                        - $match_type_count_for{1} - $match_type_count_for{2};
    print "\t\t\tUNMATCHED READS            : $unmatched_reads\n";
    #####

    return;
}   # ----------  end of subroutine censing_reads_from  ----------


# Match the read and its reverse complement against each CDS
sub get_matches_for {
    my ( $query_seq,  $comp_rev_query_seq, $ELAND_filename) = @_;
    my $query_regex = qr/$query_seq/;             #Precompiled to gain speed
    my $cr_query_regex = qr/$comp_rev_query_seq/; #Precompiled to gain speed
    my $number_of_matches = 0;
    my $single_matched_CDS;

    foreach my $CDS (keys %CDS_sequences_for) {
        if ($CDS_sequences_for{$CDS} =~ $query_regex or $CDS_sequences_for{$CDS} =~ $cr_query_regex) {
            ++$number_of_matches;
            $single_matched_CDS = $CDS;
        }
        # If read matches multiple CDSs exit process earlier
        return 2 if $number_of_matches > 1; #Multiple matches
    }

    # if single match
    if ($single_matched_CDS) {
        ++$census_results_for{$single_matched_CDS}{$ELAND_filename};
        return 1;
    }
    
    return 0; #0 matches
}   # ----------  end of subroutine get_matches_for  ----------


# This subroutine prints the first three rows of the LOX file.
# However, second and third row should be modified manually
# at the end of the process.
sub printing_header_rows_for {
    my ($LOX) = @_;
    my $header_row_1;
    my $header_row_2;
    my $header_row_3;

    foreach my $treatment (1..scalar @ELAND_filenames) {
        $header_row_1 .= "\t$ELAND_filenames[$treatment-1]";
        $header_row_2 .= "\tx$treatment";
        $header_row_3 .= "\tx0$treatment";
    }

    print {$LOX} "\t$header_row_1\n";
    print {$LOX} "\t$header_row_2\n";
    print {$LOX} "ID\tCommon Name" . "$header_row_3";
	
    return ;
}   # ----------  end of subroutine printing_header_rows_for  ----------


#Unfiltered ELAND alignment output
#Each line of the output file contains the following fields:
#1. Sequence name (derived from file name and line number if format is not fasta)
#2. Sequence
#3. Type of match codes:
#• NM—No match found
#• QC—No matching done: QC failure (too many Ns)
#• RM—No matching done: repeat masked (may be seen if repeatFile.txt was
#specified)
#• U0—Best match found was a unique exact match
#• U1—Best match found was a unique 1-error match
#• U2—Best match found was a unique 2-error match
#• R0—Multiple exact matches found
#• R1—Multiple 1-error matches found, no exact matches
#• R2—Multiple 2-error matches found, no exact or 1-error matches
#4. Number of exact matches found
#5. Number of 1-error matches found
#6. Number of 2-error matches found
#7. The following fields are only used if a unique best match was found:
#8. Genome file in which match was found
#9. Position of match (bases in file are numbered starting at 1)
#10. Direction of match (F=forward strand, R=reverse)
#11. How N characters in read were interpreted (“.”=not applicable, “D”=Detection,
#“I”=Insertion)
#The following field is only used in the case of a unique inexact match:
#12. Position and type of first substitution error (A numeral refers to a run of matching
#bases, an upper case base or N refers to a base in the reference that differs from the
#read. For example, 11A: after 11 matching bases, base 12 is A in the reference but
#not in the read)

