use strict;

use warnings;

use FindBin;

use Pod::Usage;

use Getopt::Long qw(:config posix_default no_ignore_case bundling pass_through);

use Data::Dumper;

use List::Util qw (min max);

use File::Basename;

use File::Spec;

use Carp;

use Digest::MD5;

our $ROOT_DIR;

BEGIN {

$ROOT_DIR = File::Spec->rel2abs(File::Spec->catdir($FindBin::RealBin, File::Spec->updir()));

}

use lib ("$ROOT_DIR/PerlLib");

use POSIX qw(ceil);

use Gene_obj;

use Nuc_translator;

use Fasta_reader;

use Longest_orf;

use Pipeliner;

use DelimParser;

use Cwd;

my $RETAIN_LONG_ORFS_MIN_LENGTH = 1000000; # so essentially, off by default

srand(1234);

my $genetic_code = "Universal";

my $genetic_code_options = join("\n", &Nuc_translator::get_genetic_codes());

my $usage = <<__EOUSAGE__;

;

my $UTIL_DIR = "$ROOT_DIR/util";

$ENV{PATH} = "$UTIL_DIR/bin:$ENV{PATH}";

my ($transcripts_file,$train_file);

my $top_ORFs_train = 500;

my $max_prot_length_for_training = 5000;

my $help;

my $verbose;

my $search_pfam = "";

my ($reuse,$pfam_out);

my $RETAIN_LONG_ORFS_MODE = 'dynamic';

my @GC_TO_MIN_LONG_ORF_LENGTH = ( [25, 519],

[30, 552],

[35, 573],

[40, 615],

[45, 663],

[50, 726],

[55, 804],

[60, 927],

[65, 1083],

[70, 1347],

[75, 1746],

[80, 2436] );

my $retain_pfam_hits_file;

my $retain_blastp_hits_file;

my $cpu = 1;

my $MPI_DEBUG = 1;

my $single_best_flag = 0; # deprecated: single best is the default

my $retain_multiple_orfs_flag = 0;

my $start_contrast_alpha = 0.4;

my $start_kozak_weight = 25;

my $short_orf_length = 100;

my $short_orf_min_contrast = 200;

my $NO_REFINE_START_CODONS_FLAG = 0;

my $output_dir = &Pipeliner::ensure_full_path(cwd()); # current working directory by default.

&GetOptions( 't=s' => \$transcripts_file,

'h' => \$help,

'v' => \$verbose,

'T=i' => \$top_ORFs_train,

'search_pfam=s' => \$search_pfam,

'reuse' => \$reuse,

'retain_long_orfs_mode=s' => \$RETAIN_LONG_ORFS_MODE,

'retain_long_orfs_length=i' => \$RETAIN_LONG_ORFS_MIN_LENGTH,

'debug' => \$MPI_DEBUG,

'retain_pfam_hits=s' => \$retain_pfam_hits_file,

'retain_blastp_hits=s' => \$retain_blastp_hits_file,

'cpu=i' => \$cpu,

'single_best_only' => \$single_best_flag,

'retain_multiple_orfs' => \$retain_multiple_orfs_flag,

'start_contrast_alpha=f' => \$start_contrast_alpha,

'start_kozak_weight=f' => \$start_kozak_weight,

'short_orf_length=i' => \$short_orf_length,

'short_orf_min_contrast=f' => \$short_orf_min_contrast,

'genetic_code|G=s' => \$genetic_code,

'no_refine_starts' => \$NO_REFINE_START_CODONS_FLAG,

'output_dir|O=s' => \$output_dir,

);

if ($help) {

die $usage;

}

if (@ARGV) {

die "Error, don't understand options: @ARGV";

}

our $SEE = $verbose;

unless ($transcripts_file && -s $transcripts_file) {

die $usage;

}

if ($genetic_code) {

$genetic_code = " --genetic_code $genetic_code";

}

main: {

if ($transcripts_file =~ /\.gz$/) {

# create an uncompressed version locally and use that instead

my $uncompressed_transcripts_file = basename($transcripts_file);

$uncompressed_transcripts_file =~ s/\.gz$//;

if (! -s $uncompressed_transcripts_file) {

&process_cmd("gunzip -c $transcripts_file > $uncompressed_transcripts_file");

}

$transcripts_file = $uncompressed_transcripts_file;

}

my $workdir = &Pipeliner::ensure_full_path("$output_dir/" . basename($transcripts_file) . ".transdecoder_dir");

unless (-d $workdir) {

die "Error, cannot find directory: $workdir, be sure to first run TransDecoder.LongOrfs before TransDecoder.Predict\n\n";

}

my $checkpoints_dir = "$workdir/__checkpoints_TDpredict";

if (! -d $checkpoints_dir) {

mkdir($checkpoints_dir);

}

my $pipeliner = new Pipeliner('-verbose' => 2, '-checkpoint_dir' => $checkpoints_dir);

my $prefix = "$workdir/longest_orfs";

my $cds_file = "$prefix.cds";

my $gff3_file = "$prefix.gff3";

my $pep_file = "$prefix.pep";

## Train a Markov model based on user-provided file or longest candidate CDS sequences, score all candidates, and select the final set.

my $top_cds_file = "$cds_file.top_${top_ORFs_train}_longest";

{

# to speed things up only check for redundancy up to x the number of entries we want

my $red_num = $top_ORFs_train * 10;

my $red_num_cds_longest_file = "$cds_file.top_longest_${red_num}";

$pipeliner->add_commands(new Command("$UTIL_DIR/get_top_longest_fasta_entries.pl $cds_file $red_num $max_prot_length_for_training > $red_num_cds_longest_file",

"get_longest_orfs.ok"));

$pipeliner->add_commands(new Command("$UTIL_DIR/exclude_similar_proteins.pl $red_num_cds_longest_file $genetic_code > $red_num_cds_longest_file.nr", "nr.ok"));

$pipeliner->add_commands(new Command("$UTIL_DIR/get_top_longest_fasta_entries.pl $red_num_cds_longest_file.nr $top_ORFs_train > $top_cds_file",

"top_train_select.ok"));

$pipeliner->run();

}

# get hexamer scores

my $hexamer_scores_file = "$workdir/hexamer.scores";

my $hexamer_checkpoint = "hexamer_scores_file.ok";

my $base_freqs_file = "$workdir/base_freqs.dat";

if ($RETAIN_LONG_ORFS_MODE eq 'dynamic') {

if (&genetic_code_has_standard_stops()) {

$RETAIN_LONG_ORFS_MIN_LENGTH = &get_dynamic_retain_long_orf_length($base_freqs_file, \@GC_TO_MIN_LONG_ORF_LENGTH);

}

else {

# the GC table was simulated with TAA/TAG/TGA; with fewer stops random ORFs are much

# longer (e.g. TGA-only ciliate code), so no ORF is retained on length alone

print STDERR "-genetic code stop codons differ from TAA/TAG/TGA: no length-based ORF retention\n";

$RETAIN_LONG_ORFS_MIN_LENGTH = 1000000000;

}

}

# score kmers for Markov model

$pipeliner->add_commands(new Command("$UTIL_DIR/seq_n_baseprobs_to_loglikelihood_vals.pl $top_cds_file $base_freqs_file $genetic_code > $hexamer_scores_file",

"hexamer_scores.ok") );

# score all cds entries

my $cds_scores_file = "$cds_file.scores";

my $cds_scores_checkpoint = "$cds_scores_file.ok";

$pipeliner->add_commands(new Command("$UTIL_DIR/score_CDS_likelihood_all_6_frames.pl $cds_file $hexamer_scores_file $genetic_code > $cds_scores_file",

"score_cdss.ok"));

$pipeliner->run();

##############

## select orfs

## Steps below may be rerun if input parameters are set differently (such as including pfam or blast results in a subsequent run).

my @checkpoint_inputs = ($transcripts_file, "short_orf:$short_orf_length:$short_orf_min_contrast");

if ($retain_pfam_hits_file) {

push (@checkpoint_inputs, $retain_pfam_hits_file);

}

if ($retain_blastp_hits_file) {

push (@checkpoint_inputs, $retain_blastp_hits_file);

}

my $big_checkpoint_string = join("^", @checkpoint_inputs);

my $checkpoint_hashcode = Digest::MD5::md5_hex($big_checkpoint_string);

my $select_cmd = "$UTIL_DIR/select_best_ORFs_per_transcript.pl --gff3_file $gff3_file --cds_scores $cds_scores_file "

. " --min_length_auto_accept $RETAIN_LONG_ORFS_MIN_LENGTH "

. " --training_orfs $top_cds_file --hexamer_scores $hexamer_scores_file $genetic_code "

. " --short_orf_length $short_orf_length --short_orf_min_contrast $short_orf_min_contrast ";

if ($retain_pfam_hits_file) {

$select_cmd .= " --pfam_hits $retain_pfam_hits_file ";

}

if ($retain_blastp_hits_file) {

$select_cmd .= " --blast_hits $retain_blastp_hits_file ";

}

unless ($retain_multiple_orfs_flag) {

$select_cmd .= " --single_best_orf ";

}

$select_cmd .= " > $cds_file.best_candidates.gff3 ";

$pipeliner->add_commands(new Command($select_cmd, "select_best_orfs.$checkpoint_hashcode.ok", "Selecting best orfs"));

my $final_output_prefix = &Pipeliner::ensure_full_path("$output_dir/" . basename($transcripts_file) . ".transdecoder");

my $target_final_file = "$cds_file.best_candidates.gff3";

## start codon refinement

unless ($NO_REFINE_START_CODONS_FLAG) {

my $cmd = "$UTIL_DIR/train_kozak_model.pl --transcripts $transcripts_file --training_orfs $top_cds_file "

. " --out_prefix $workdir/start_refinement";

$pipeliner->add_commands(new Command($cmd, "train_kozak_model.$checkpoint_hashcode.ok", "Training start codon (Kozak) context model"));

my $revised_orf_gff3_file = "$target_final_file.revised_starts.gff3";

$cmd = "$UTIL_DIR/start_codon_refinement.pl --transcripts $transcripts_file --gff3_file $target_final_file --workdir $workdir "

. " --training_orfs $top_cds_file --hexamer_scores $hexamer_scores_file $genetic_code "

. " --kozak_model $workdir/start_refinement.kozak.tsv "

. " --contrast_alpha $start_contrast_alpha --kozak_weight $start_kozak_weight > $revised_orf_gff3_file";

$pipeliner->add_commands(new Command($cmd, "revise_starts.$checkpoint_hashcode.ok", "Refining start codon selections."));

$target_final_file = $revised_orf_gff3_file;

}

$pipeliner->add_commands(new Command("cp $target_final_file $final_output_prefix.gff3",

"copy_to_final_outputfile.$checkpoint_hashcode.ok",

"copying output to final output file: $final_output_prefix.gff3") );

## write final outputs:

{

## make a BED file for viewing in IGV

my $gff3_file = "$final_output_prefix.gff3";

my $bed_file = $gff3_file;

$bed_file =~ s/\.gff3$/\.bed/;

my $cmd = "$UTIL_DIR/gff3_file_to_bed.pl $gff3_file > $bed_file";

$pipeliner->add_commands(new Command($cmd, "make_final_bed.$checkpoint_hashcode.ok", "Making bed file: $bed_file"));

# make a peptide file:

my $best_pep_file = $gff3_file;

$best_pep_file =~ s/\.gff3$/\.pep/;

$cmd = "$UTIL_DIR/gff3_file_to_proteins.pl --gff3 $gff3_file --fasta $transcripts_file $genetic_code > $best_pep_file";

$pipeliner->add_commands(new Command($cmd, "make_final_pep.$checkpoint_hashcode.ok", "Making pep file: $best_pep_file"));

# make a CDS file:

my $best_cds_file = $best_pep_file;

$best_cds_file =~ s/\.pep$/\.cds/;

$cmd = "$UTIL_DIR/gff3_file_to_proteins.pl --gff3 $gff3_file --fasta $transcripts_file --seqType CDS $genetic_code > $best_cds_file";

$pipeliner->add_commands(new Command($cmd, "make_final_cds.$checkpoint_hashcode.ok", "Making cds file: $best_cds_file"));

}

$pipeliner->run();

print STDERR "transdecoder is finished. See output files $final_output_prefix.\*\n\n\n";

exit(0);

}

sub process_cmd {

my ($cmd) = @_;

print "CMD: $cmd\n";

my $ret = system($cmd);

if ($ret) {

die "Error, cmd: $cmd died with ret $ret";

}

return;

}

sub parse_pfam_hits {

my ($pfam_hits_file) = @_;

my %has_pfam_hit;

if (! -e $pfam_hits_file) {

die "Error, cannot find pfam hits file: $pfam_hits_file";

}

print "PFAM output found and processing...\n";

# capture those proteins having pfam hits

open (my $fh, $pfam_hits_file) or die "Error, cannot open file: $pfam_hits_file";

while (my $ln=<$fh>) {

next if $ln=~/^\#/;

my @x = split(/\s+/,$ln);

next unless $x[3]; # domtbl

my $orf_acc = $x[3];

$has_pfam_hit{$orf_acc} = 1;

}

close $fh;

return(%has_pfam_hit);

}

sub parse_blastp_hits_file {

my ($blastp_file) = @_;

unless (-e $blastp_file) {

die "Error, cannot find file $blastp_file";

}

my %blastp_hits;

open (my $fh, $blastp_file) or die "Error, cannot open file $blastp_file";

while (<$fh>) {

chomp;

my @x = split(/\t/);

my $id = $x[0];

$blastp_hits{$id} = 1;

}

close $fh;

return(%blastp_hits);

}

sub check_program() {

my @paths;

foreach my $prog (@_) {

my $path = `which $prog`;

unless ($path =~ /\w/) {

die "Error, path to a required program ($prog) cannot be found\n\n"

}

chomp($path);

$path = readlink($path) if -l $path;

push( @paths, $path );

}

return @paths;

}

sub genetic_code_has_standard_stops {

# $genetic_code holds the option string passed to the util scripts

my ($code) = $genetic_code =~ /--genetic_code\s+(\S+)/;

&Nuc_translator::use_specified_genetic_code($code || "Universal");

my $stops = join(",", sort grep { !/U/ } &Nuc_translator::get_stop_codons()); # DNA spellings

&Nuc_translator::use_specified_genetic_code("Universal");

return($stops eq "TAA,TAG,TGA");

}

sub get_dynamic_retain_long_orf_length {

my ($base_freqs_file, $GC_to_min_orf_length_aref) = @_;

my $pct_GC = -1;

open(my $fh, $base_freqs_file) or confess "Error, cannot open file: $base_freqs_file";

while (<$fh>) {

chomp;

my ($base, $count, $ratio) = split(/\t/);

if ($base eq 'C') {

$pct_GC = 2 * $ratio * 100;

last;

}

}

close $fh;

print STDERR "PCT_GC: $pct_GC\n";

unless ($pct_GC > 0) {

confess "Error, didn't parse percent C from file: $base_freqs_file";

}

foreach my $bin (@$GC_to_min_orf_length_aref) {

my ($gc, $min_len) = @$bin;

if ($pct_GC <= $gc) {

return($min_len);

}

}

# if we got here, then GC content must be crazy high!

# in this case, rely entirely on Markov model for coding determination.

return(1000000); #effectively infinity here.

}