#!/usr/bin/env perl
# $Id: prepare_data,v 1.7 2007/03/26 17:38:06 c4chris Exp $

################################################################################
#
# prepare_data
# ------------
#
# Claudio Lottaz, SIB-ISREC, Claudio.Lottaz@isb-sib.ch
# Christian Iseli, LICR ITO, Christian.Iseli@licr.org
#
# Copyright (c) 2006 Swiss Institute of Bioinformatics. All rights reserved.
#
################################################################################

use strict;
use FASTAFile;
use Symbol;

# global variables
my $verbose = 1;
my $dosplit = 0;
my $forcedminmask = undef;

my $datadir = '.';
my $filestem = '';

require "build_model_utils.pl";

################################################################################
#
#   Check command-line for switches
#

my $usage = "Usage: prepare_data [options] <conf-files>\n" .
    "       where options are:\n" .
    "        -q       don't log on terminal\n" .
    "        -e       split extracted data into training and test sets\n" .
    "        -m <int> force minimal mask, overwrites entry in config-files\n" .
    "More information is obtained using 'perldoc prepare_data'\n";

while ($ARGV[0] =~ m/^-/) {
    if ($ARGV[0] eq '-q') { shift; $verbose = 0; next; }
    if ($ARGV[0] eq '-e') { shift; $dosplit = 1; next; }
    if ($ARGV[0] eq '-m') { shift; $forcedminmask      = shift; next; }
    die "Unrecognized switch: $ARGV[0]\n$usage";
}
if ($#ARGV < 0) { die "No configuration file specified\n$usage"; }

################################################################################
#
#   Main-loop through all specified config-files
#

my $parFile;
while($parFile = shift) {
    my($organism, $hightaxo, $dbfiles, $ugdata, $estdata, $datadir2, $filestem2,
       $rnafile, $estfile, $estcdsfile, $estutrfile, $trainingfile, $testfile,
       $utrfile, $cdsfile, $tuplesize, $minmask, $pseudocounts, $minscore,
       $startlength, $startpreroll, $stoplength, $stoppreroll, $smatfile,
       $estscanparams, $nb_isochores, $isochore_borders) =
	   readConfig($parFile, undef, $forcedminmask, undef,
		      undef, undef, undef,
		      undef, undef, $verbose);
    log_open("readconfig.log");
    $datadir = $datadir2; $filestem = $filestem2;
    showConfig($parFile, $organism, $hightaxo, $dbfiles, $ugdata, $estdata, $datadir,
	       $rnafile, $estfile, $estcdsfile, $estutrfile, $trainingfile,
	       $testfile, $utrfile, $cdsfile, $tuplesize, $minmask,
	       $pseudocounts, $minscore, $startlength, $startpreroll,
	       $stoplength, $stoppreroll, $smatfile, $estscanparams,
	       $nb_isochores, $isochore_borders);
    log_close();

    # prepare data
    log_open("prepare_data.log");
    my $gc_histogram = analyzeGC($rnafile);
    my $isochores = computeIsochores($nb_isochores, $isochore_borders, $gc_histogram);

    if ($dosplit) {
      splitTraining($rnafile, $trainingfile, $testfile, $isochores);
    } else {
      symlink $rnafile, $trainingfile;
      symlink $rnafile, $testfile;
      splitIsochores($rnafile, $isochores);
    }
    maskRedundancy($minmask, $isochores);
    log_print("\nGenerating evaluation mRNA data....");
    split_mRNAs($testfile, $utrfile, $cdsfile);
    log_close();

    print "$parFile done.\n";
}

exit(0);


################################################################################
#
#   Analyze GC-content distribution
#

sub analyzeGC {
    # crawl through the RNAs collected mRNA and establish a GC-contents histogram

    my($rnafile) = @_;
    my $datfh = gensym;

    log_print("\nAnalyzing GC contents...");

    my $nb_seqs = 0;
    my @gc_histogram;
    if (-s "$datadir/Report/gc.dat") {
	log_print(" - loading existing GC-content histogram...");
	open($datfh, "$datadir/Report/gc.dat");
	while(<$datfh>) {
	    my($index, $count) = split;
	    $gc_histogram[$index] = $count;
	    $nb_seqs += $count;
	}
	close($datfh);
	log_print("   read $nb_seqs sequences");
    } else {
	log_print(" - generating GC-content histogram...");
	my $e;
	my $src = FASTAFile->new("$rnafile"); $src->openStream;
	while(defined ($e = $src->getNext)) {
	    my $gc = int(gc_content($e->{_seq}));
	    $gc_histogram[$gc] += 1;
	    $nb_seqs += 1;
	}
	close($src->{_BTFfile});
	log_print("   read $nb_seqs sequences");

	# write datafile
	open($datfh, ">$datadir/Report/gc.dat");
	for (my $i = 0; $i < 100; $i++) {
	    if (!defined $gc_histogram[$i]) { print $datfh $i, " 0\n"; }
	    else { print $datfh $i, " ", $gc_histogram[$i], "\n"; }
	}
	close($datfh);
    }

    # write the gnu-script for the gc-histogram
    my $scriptfh = gensym;
    open($scriptfh, ">$datadir/Report/gc.gplot");
    print $scriptfh <<EndOfBuffer;
set title "GC-content Histogram"
set xlabel "percentage"
set ylabel "frequency"
set style data boxes
unset key
plot 'gc.dat'
EndOfBuffer
    close($scriptfh);

    log_print(" - written datafile and gnuplot script");
    return \@gc_histogram;
}

sub gc_content {
    # expectes a nucleotide sequence as entry and computes its
    # gc-content in percentage
    my($seq) = @_;
    my $a = $seq =~ s/A/A/g;
    my $c = $seq =~ s/C/C/g;
    my $g = $seq =~ s/G/G/g;
    my $t = $seq =~ s/T/T/g;
    return 100*($c + $g) / ($a + $c + $g + $t);
}

################################################################################
#
#   Split mRNA data into isochores, test and training sequences
#

sub computeIsochores {
    # according to the GC-contents histogram compute the isochore borders
    # such that close to equal numbers of sequences are contained in the
    # $nb_isochores isochores. If $nb_isochores is zero use the isochore
    # borders defined explicitely.

    my(@isochores, $i);
    my ($nb_isochores, $isochore_borders_ref, $gc_histogram_ref) = @_;
    my @isochore_borders = @{$isochore_borders_ref};
    my @gc_histogram = @{$gc_histogram_ref};
    my $nb_seqs = 0;
    for ($i = 0; $i<= $#gc_histogram; $i++) { $nb_seqs += $gc_histogram[$i]; }

    # compute isochore borders if to be done automatically
    log_print(" - computing isochore borders...");
    if ($nb_isochores == 0) { log_print("   skipped"); }
    else {
	undef @isochore_borders;
	$isochore_borders[0] = 0;
	my $cumulator = 0;
	my $isoIncrement = $nb_seqs/$nb_isochores;
	my $currentTop = $isoIncrement;
	for (my $i = 0; $i < 100; $i++) {
	    $cumulator += $gc_histogram[$i];
	    if ($cumulator > $currentTop) {
		push(@isochore_borders, $i);
		$currentTop += $isoIncrement;
	    }
	}
	push(@isochore_borders, 100);
    }
    my($buffer) = "   isochores used: ";
    for ($i = 0; $i < $#isochore_borders; $i++) {
	$isochores[$i] = $isochore_borders[$i] . "-" . $isochore_borders[$i+1];
	$buffer .= $isochores[$i] . ", ";
    }
    chop($buffer);chop($buffer);
    log_print("$buffer");
    return \@isochores;
}

sub splitTraining {
    # Splits the mRNA data into training set and test set. Moreover,
    # generates the isochore partitionning for the training set.

    my($rnafile, $trainingfile, $testfile, $isochores_ref) = @_;
    my @isochores = @{$isochores_ref};
    my $isodir = "$datadir/Isochores";

    log_print("\nSplit mRNA data into isochores, test and training data...");

    log_print(" - writing isochores...");
    my($isochore, @isochorefhs, $isochorefh);

    # check whether isochores and testfile have already been computed
    my($isochoresReady) = 1;
    foreach $isochore (@isochores) {
	if (!(-s "$isodir/mrna$isochore.seq")) { $isochoresReady=0; next; }
	log_print("   $isodir/mrna$isochore.seq already exists");
    }
    if (-s "$trainingfile") { log_print("   $trainingfile already exists"); }
    else { $isochoresReady = 0; }
    if (-s "$testfile") { log_print("   $testfile already exists"); }
    else { $isochoresReady = 0; }

    # write trainingfile, testfile and isochores (CDS)
    if ($isochoresReady == 1) { log_print("   all files exist, skipped"); }
    else{
	# open all isochore files for writing
	foreach $isochore (@isochores) {
	    $isochorefh = gensym;
	    open($isochorefh, ">$isodir/mrna$isochore.seq");
	    push(@isochorefhs, $isochorefh);
	}
	my $trainingfh = gensym; open($trainingfh, ">$trainingfile");
	my $testfh = gensym; open($testfh, ">$testfile");

	# read mRNAs and write the isochores
	my $e;
	my $testSeqs = 0;
	my %isoSeqs;
	my $src = FASTAFile->new("$rnafile");
	$src->openStream;
	while(defined ($e = $src->getNext)) {
	    my $gc = gc_content($e->{_seq});
	    for (my $i = 0; $i <= $#isochores; $i++) {
		my($low, $high) = ($isochores[$i] =~ m/^([^\-]+)\-(.*)$/);
		if ($gc < $high) {
		    my($id, $x, $begin, $end, $desc)
		      = split(/ /, $e->{_seqHead}, 5);
		    $isoSeqs{$isochores[$i]}++;
		    if ($isoSeqs{$isochores[$i]}%2) {
		      $e->printFASTA($testfh);
		      $testSeqs++;
		    } else {
			$e->printFASTA($trainingfh);
			my $fh = $isochorefhs[$i];
			print $fh &genRNAEntry($e->ac, $desc, $begin,
					       $end, $e->{_seq});
		    }
		    last;
		}
	    }
	}
	close($src->{_BTFfile});

	# close the isochore files
	foreach $isochorefh (@isochorefhs) { close($isochorefh); }
	close($trainingfh);
	close($testfh);
	foreach (sort(keys(%isoSeqs))) {
	    log_print("   ", $isoSeqs{$_}, " sequences found in isochore $_");
	}
	log_print("   $testSeqs of these written into $testfile");
    }
}

sub splitIsochores {
    # Split the mRNAs into isochores (CDS).

    my($rnafile, $isochores_ref) = @_;
    my @isochores = @{$isochores_ref};
    my($isochore, @isochorefhs, $isochorefh);
    my $isodir = "$datadir/Isochores";

    log_print("\nSplit mRNA data into isochores...");

    # check whether isochores have already been computed
    my($isochoresReady) = 1;
    foreach $isochore (@isochores) {
	if (!(-s "$isodir/mrna$isochore.seq")) { $isochoresReady=0; next; }
	log_print("   isochore $isodir/mrna$isochore.seq already exists");
    }

    # write isochores (CDS)
    if ($isochoresReady == 1) { log_print("   all files exist, skipped"); }
    else{
	# open all isochore files for writing
	foreach $isochore (@isochores) {
	    $isochorefh = gensym;
	    open($isochorefh, ">$isodir/mrna$isochore.seq");
	    push(@isochorefhs, $isochorefh);
	}

	# read mRNA and write the isochores
	my $e;
	my %isoSeqs;
	my $src = FASTAFile->new("$rnafile");
	$src->openStream;
	while(defined ($e = $src->getNext)) {
	    my $gc = gc_content($e->{_seq});
	    for (my $i = 0; $i <= $#isochores; $i++) {
		my($low, $high) = ($isochores[$i] =~ m/^([^\-]+)\-(.*)$/);
		if ($gc < $high) {
		    $isoSeqs{$isochores[$i]}++;
		    my($id, $x, $begin, $end, $desc)
		      = split(/ /, $e->{_seqHead}, 5);
		    my $fh = $isochorefhs[$i];
		    print $fh &genRNAEntry($e->ac, $desc, $begin,
					   $end, $e->{_seq});
		    last;
		}
	    }
	}
	close($src->{_BTFfile});

	# close the isochore files
	foreach $isochorefh (@isochorefhs) { close($isochorefh); }
	foreach (sort(keys(%isoSeqs))) {
	    log_print("   ", $isoSeqs{$_}, " sequences found in isochore $_");
	}
    }
}

sub genRNAEntry {
    my($id, $desc, $cdsBegin, $cdsEnd, $seq) = @_;
    $seq =~ s/(.{80})/$1\n/g;
    return ">tem|$id CDS: $cdsBegin $cdsEnd $desc\n$seq\n";
}

sub maskRedundancy {
    # replace redundant regions by 'N'

    my($minmask, $isochores_ref) = @_;
    my @isochores = @{$isochores_ref};
    log_print("\nMasking redundancy from isochores...");

    my $isochore;
    my $isodir = "$datadir/Isochores";
    foreach $isochore (@isochores) {
	my $maskedFile = "$isodir/mrna$isochore\_mr$minmask.seq";
	if (-s $maskedFile) {
	  log_print(" - $maskedFile already exists, skipped");
	} else {
	  my $infile = "$isodir/mrna$isochore.seq";
	  log_print(" - masking redundancy in isochore $isochore");
	  system("maskred -m $minmask < $infile > $maskedFile");
	}
	my $masked = `tail -1 $maskedFile`;
	$masked =~ m/^>masked nucleotides: (\d+)/;
	my $m = $1;
	my $nts = `grep -v '^>' $maskedFile | wc -cl`;
	$nts =~ s/^\s*(\d+)\s+(\d+)\s*$/$2-$1/e;
	my $pct = int(10000*$m/$nts + 0.5) / 100;
	log_print("   masked $m of $nts nucleotides ($pct%)");
    }
}

sub split_mRNAs {
    # reads testfile and splits the entries into untranslated and
    # coding sequences according to the annotation expected in the
    # FASTA headers.

    my($testfile, $utrfile, $cdsfile) = @_;

    log_print(" - splitting mRNAs into UTRs and CDSs...");
    if ((-s $utrfile) && (-s $cdsfile)) { 
	log_print("   UTR and CDS files already exist, skipped");
	return;
    }

    my $e;
    my $utrfh = gensym; open($utrfh, ">$utrfile");
    my $cdsfh = gensym; open($cdsfh, ">$cdsfile");
    my $src = FASTAFile->new($testfile); $src->openStream;
    while(defined($e = $src->getNext)) {
	if ($e->{_seqHead} =~ m/CDS: (\S+) (\S+)/) {
	    my $cdsStart = $1; my $cdsEnd = $2;
	    my $cdsLen = $cdsEnd - $cdsStart + 1;
	    my $head = $e->{_seqHead}; my $seq = $e->{_seq};
	    $e->{_seqHead} =~ s/CDS: (\S+) (\S+)/CDS: 1 $cdsLen/;
	    $e->{_seq} = substr($e->seq, $cdsStart - 1, $cdsLen);
	    $e->printFASTA($cdsfh);
	    $e->{_seqHead} = "$head, 5'UTR";
	    $e->{_seq} = substr($seq, 0, $cdsStart - 1);
	    $e->printFASTA($utrfh);
	    $e->{_seqHead} = "$head, 3'UTR";
	    $e->{_seq} = substr($seq, $cdsEnd);
	    $e->printFASTA($utrfh);
	}
    }
    close($src->{_BTFfile});
    close($cdsfh);
    close($utrfh);

    my $nbUTRs = `grep -c '^>' $utrfile`; chop($nbUTRs);
    my $nbCDSs = `grep -c '^>' $cdsfile`; chop($nbCDSs);
    log_print("   found $nbUTRs untranslated regions and $nbCDSs coding sequences");
}

################################################################################
#
#   Documentation
#

=head1 NAME

prepare_data - prepare training and test data for ESTScan

=head1 SYNOPSIS

prepare_data [options] <config-files...>

=head1 DESCRIPTION

prepare_data prepares training and test data to generate codon usage
tables for ESTScan.
The script reads
configuration files on the comannd line and performs the following
steps for each of them:

 - Split the data into isochores and testsets
 - Mask redundant pieces of sequences
 - Extract untranslated regions from test mRNA

Files which already exist are reused. If an existing file is to be
recomputed, it must be deleted before the script is run again. If a
particular collection of mRNA is to be used instead of data extracted
by extract_mRNA, providing these in FASTA format under the name of the
mRNA file (where extract_mRNA would store the extracted data) is
enough.

If the '-e' switch is provided, prepare_data splits the data extracted
from the given databases or provided by the user into training and
test data.

=head1 DIRECTORY STRUCTURE

extract_data uses the following directory structure, the root of which
is given in the configuration file. From this root it contains the
following subdirectories:

 - Isochores: data split in isochores, some with redundancy masked
 - Report: contains all log files

mRNA data, test and training data files, is deposited in the data-root
directory if not otherwise specified in the configuration file.

=head1 OPTIONS

 -q quiet
    Do not log on terminal.

 -e split data into test and training
    If this switch is provided extracted mRNA data is split into
    training and test set. As sequences are extracted they are
    alternately deposited in the training and the test set..

 -m <minmask>
    Minimum consecutive runs of nucleotides masked from the original
    sequences in order to limit data redundancy. Pieces of sequence
    are only masked, if all of their nucleotides are part of
    reoccuring 12-tuples which overlap by at least 4 nucleotides. This
    switch overwrites the variable $minmask from paramter files.

=head1 CONFIGURATION FILE

The parameters defined in the configuration file have the
following meaning:

 * $organism (mandatory)
    The desired organism as it is given in EMBL "OS" or RefSeq
   "ORGANISM" lines.

 * $dbfiles
   Files from where full-length mRNA sequences are to be extracted,
   tries to guess whether the files come from EMBL or RefSeq. If
   this is not specified, expects a collection of mRNA in $rnafile.

 * $datadir (mandatory)
   Base directory where all of the above files are located
   and the temporary result files are stored

 * $rnafile (default is "$datadir/mrna.seq")
   Name of the file with the extracted mRNA entries.

 * $smatfile (default is "$datadir/Matrices/$filestem.smat")
   Name of the file where the HMM-model is to be written

 * $nb_isochores (default is 0)
   Number of isochores, when isochores are to be determined
   automatically from the GC-content distributaion as equal-sized
   groups. 0 means,no automatic detection.

 * @isochore_borders (default is (0, 43, 47, 51, 100)
   Array of GC percentages where isochores are split, first entry is
   usually 0 and last 100. This is overwritten when $nb_isochores is
   not zero.

 * $tuplesize (default is 6)
   Size of tuples counted for codon statistics. This is overwritten by
   the -t switch.

 * $minmask (default is 30)
   Minimal run of consecutive nucleotides masked as redundant. This is
   overwritten by the -m switch.

 * $pseudocounts (default is 1)
   pseudocount to be added when generating the codon usage tables,
   overwritten by -m

 * $minscore (default is -100)
   minimum score attributed to log-odds and log-probabilities,
   overwritten by -s.

 * $startlength, $startpreroll (default 2+ceil(tuplesize/3) and 2)
   number of nucleotide triplets contained in the start profile and
   how many of these are contained in the 5' untranslated region,
   overwritten by -l (length) and -r (preroll)

 * $stoplength, $stoppreroll (default 2+ceil(tuplesize/3) and 2)
   number of nucleotide triplets contained in the stop profile and
   how many of these are contained in the coding sequence,
   overwritten by -L (length) and -R (preroll)

 * $estscanparams (default is "-m -50 -d -50 -i -50 -N 0")
   parameters passed to ESTScan during evaluation

$filestem is used to generate many filenames. It is generated
automatically according to the tuplesize, the minmask and the
pseudocounts applied to generate them.

=head1 AUTHOR

Claudio Lottaz, SIB-ISREC, Claudio.Lottaz@isb-sib.ch

=cut

#
#   End of file
#
################################################################################