vicsin

#!/usr/bin/env perl
$|++; # Forces stdout flush

# VICSIN Pipeline
# Copyright 2017 University of Illinois at Urbana-Champaign
# Author: Joe Leigh <jleigh@illinois.edu>
# Requirements:
#  Prodigal
#  PhiSpy
#  MUMmer
#  Spine/Agent
#  MCL
#  hmmer
#  Metagene Annotator
#  MUSCLE
#  BLAST
#  BLAST+
#  VirSorter
#  BioPerl

# Inputs:
#  1) input.txt file with genome file prefixes, like:
#    NC_004663
#    GENOMEA
#  2) config.txt:
#    Basic config:
#      input_path=
#        path to the directory containing the input gff/fasta files
#      output_path=
#        path to the base directory where all output should be saved
#    Program paths: (default for all: . or in the case of executables, PATH)
#      genbank_to_seed
#      genbank_to_fasta
#      gff_to_seed
#      prodigal
#      virsorter
#      phiSpy
#      blastn
#      spine
#      agent
#      mcl
#    Program params:
#      phispy_windowsize=
#        size of window to scan genes (set default = 40)
#      phispy_threshold=
#        number of consecutive genes required to call element (set default = 20)
#      spine_percent_input=
#        Number of genome % to include as core (default = 100)
#      spine_max_distance=
#        Max distance between elements (default = 10)
#      spine_agent_min_perc_id=
#        Minimum percent ID for matching regions (default = 85)
#      spine_agent_min_size_core=
#        Minimum core lenght (default = 10)
#      spine_core_file=
#        precomputed core file from separate/previous spine run
#      spacer_fasta_file=
#        Fasta file of spacers
#      known_viral_types=
#        Fasta of known viruses to include clusters
#      virsorter_database=
#        1 for RefseqABVir only, 2 for RefseqABVir + Viromes.
#      masking_file=
#        file of regions to ignore:
#          NC_004663	x	y
#          NC_004663	z	a
#          GENOMEA		m	n

use strict;
use Getopt::Long qw(GetOptions);
use YAML qw(LoadFile);
use File::Copy;
use File::Path qw(make_path rmtree);
use File::Spec;
use Bio::SeqIO;

#TODO only used for debugging; remove in production
use Data::Dumper;

# Constants
use define CONVERTED_INPUT_DIR => "Converted_Input_Files";
use constant OUTPUT_DIR => "Output_Files";
use constant INTERMEDIATE_OUTPUT_DIR => "Pre_Reblast_Output_Files";

use VICSIN;
use VH_VirSorter;
use VH_PhiSpy;
use VH_SpineAgent;
use VH_CRISPR;
use VH_Blast;
use VH_helpers;
use VH_ReBlast;
use VH_Cluster;
use VH_Database;

##### STEP 1: Parse arguments #####
print `python -c "import Bio"`;
# quit unless we have the correct number of args
my $num_args = $#ARGV+1;
if($num_args<1 || $ARGV[0] =~ /^-/){
	print "\nUsage: vicsin input.txt [config.txt] [options]\n";
	exit;
}

# We got at least one argument, the first one is the input text file
my $input_file = $ARGV[0];

# If we got a second argument, it's the config file
my $config_file = "";
if($num_args>1 && $ARGV[1] =~ /^[^-]/){
	$config_file = $ARGV[1];
}
VICSIN::parseParameters($config_file);

VH_helpers::log("Arguments:",1);
VH_helpers::log("\tInput file: $input_file",1);
VH_helpers::log("\tConfig file: $config_file",1);
VICSIN::logParameters();

if(VICSIN::param("stop") eq "parameters"){
	VH_helpers::log("Halting as requested after parsing parameters");
	exit;
}

##### STEP 2. READ input.txt #####
# Create output directory, if necessary
make_path(VICSIN::param("output_path"));

# Save each file prefix to array
my @prefixes;
open(my $fh, '<', $input_file)
	or die "Could not open input file '$input_file' $!";
while (my $row = <$fh>) {
	chomp $row;
	push @prefixes, $row;
}

# Test each file prefix for format
make_path(VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR);
VH_helpers::log("Processing input files from $input_file...");
my @valid_prefixes;
my %genomes;
my %contigs;
foreach my $prefix (@prefixes) {
	my $valid_prefix_found = 0;
	my $parse_length = 0;
	if (my $gbk_file_name = VH_helpers::file_of_type_exists($prefix,"gbk","gb")){
		$gbk_file_name = VICSIN::param("input_path")."/$gbk_file_name";
		# Skip if already done this step and running again
		if ( -f VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna" and 
			-d VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix" ){
			VH_helpers::log("\t$prefix already converted. Skipping.",1);
			$valid_prefix_found = 1;
			push @valid_prefixes, $prefix;
		} else {
			## Gbk. Generate fasta+seed
			VH_helpers::log("\t$gbk_file_name found. Converting... ",1);

			my $seed_file_name = VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix";
			my $fasta_file_name = VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna";

			# Convert gbk to seed
			VH_helpers::run_cmd("python \"".VICSIN::param('genbank_to_seed')."\" \"$gbk_file_name\" \"$seed_file_name\"");
			# Convert gbk to fasta
			VH_helpers::run_cmd(VICSIN::param('genbank_to_fasta')." -i \"$gbk_file_name\" -m genbank -s whole -a accessions -o \"$prefix.fna\"");
			# Move fasta file into new directory
			move VICSIN::param("input_path")."/$prefix.fna", $fasta_file_name;
			# TODO check formats of fasta definition lines
			$valid_prefix_found = 1;
			push @valid_prefixes, $prefix;
		}
		# TODO gather gbk metadata
		my $gbkobj = Bio::SeqIO->new(-file => $gbk_file_name,
									 -format => 'GenBank');
		my $gbkseq = $gbkobj->next_seq;
		my $gbkanno = $gbkseq->annotation;
		my @dblinks = $gbkanno->get_Annotations('dblink');
		my $gene_count = 0;

		foreach my $feat ($gbkseq->get_SeqFeatures){
			if($feat->primary_tag eq "source"){
				for my $tag ($feat->get_all_tags){
					if($tag eq "organism"){
						my @vals = $feat->get_tag_values($tag);
						$genomes{$prefix}{'organism'} = $vals[0];
					}
					if($tag eq "strain"){
						my @vals = $feat->get_tag_values($tag);
						$genomes{$prefix}{'strain'} = $vals[0];
					}
				}
			}
			if($feat->primary_tag eq "gene"){
				$gene_count++;
			}
		}

		$genomes{$prefix}{'name'} = $prefix;
		$genomes{$prefix}{'version'} = $gbkseq->version();
		$genomes{$prefix}{'format'} = 'genbank';
		$genomes{$prefix}{'definition'} = $gbkseq->desc();
		$genomes{$prefix}{'accession'} = $gbkseq->accession_number();
		$genomes{$prefix}{'keywords'} = $gbkseq->keywords();
		$genomes{$prefix}{'length'} = $gbkseq->length();
		# TODO what to do with dblink?
		$genomes{$prefix}{'genes'} = $gene_count;
	} elsif ( (my $fasta_file_name = VH_helpers::file_of_type_exists($prefix,"fna","fasta","fa")) && (my $gff_file_name = VH_helpers::file_of_type_exists($prefix,"gff")) ) {
		# Skip if already done this step and running again
		if ( -f VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna" and 
			-d VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix" ){
			VH_helpers::log("\t$prefix already converted. Skipping.",1);
			$parse_length = 1;
			$valid_prefix_found = 1;
			push @valid_prefixes, $prefix;
		} else {
			## Fasta+gff. Generate seed
			VH_helpers::log("\t$fasta_file_name and $gff_file_name found. Converting... ",1);

			$fasta_file_name = VICSIN::param("input_path")."/$fasta_file_name";
			$gff_file_name = VICSIN::param("input_path")."/$gff_file_name";
			my $seed_file_name = VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix";

			# Run gff_to_seed to generate seed file
			VH_helpers::run_cmd(VICSIN::param('gff_to_seed')." $gff_file_name $fasta_file_name $prefix");
			# Move seed output into new directory
			rmtree($seed_file_name);
			move "_SEED_$prefix", $seed_file_name or die "Could not copy _SEED_/$prefix";
			# Copy fasta file into new directory
			copy $fasta_file_name, VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna" or die "Could not copy $prefix.fna";
			# TODO check formats of fasta definition lines
			$valid_prefix_found = 1;
			$parse_length = 1;
			push @valid_prefixes, $prefix;
		}

		$genomes{$prefix}{'name'} = $prefix;
		$genomes{$prefix}{'version'} = "";
		$genomes{$prefix}{'format'} = 'fasta/gff';
		$genomes{$prefix}{'genes'} = 0;
	} elsif (my $fasta_file_name = VH_helpers::file_of_type_exists($prefix,"fna","fasta","fa")){
		# Skip if already done this step and running again
		if ( -f VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna" and 
			-d VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix" ){
			VH_helpers::log("\t$prefix already converted. Skipping.",1);
			$parse_length = 1;
			$valid_prefix_found = 1;
			push @valid_prefixes, $prefix;
		} else {
			## Fasta only. Generate gff, then generate seed
			VH_helpers::log("\t$fasta_file_name found. Converting... ",1);

			$fasta_file_name = VICSIN::param("input_path")."/$fasta_file_name";
			my $gff_file_name = VICSIN::param("input_path")."/$prefix.gff";
			my $seed_file_name = VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/_SEED_$prefix";
			
			# Run prodigal to generate gff file
			my $prod_out = VH_helpers::run_cmd(VICSIN::param('prodigal')." -f gff -c -m -i $fasta_file_name -o $gff_file_name 2>&1");
			if ($? == 0) {
				# Run gff_to_seed to generate seed file
				VH_helpers::run_cmd(VICSIN::param('gff_to_seed')." $gff_file_name $fasta_file_name $prefix");
				# Move seed output into new directory
				rmtree $seed_file_name;
				move "_SEED_$prefix", $seed_file_name or die "Could not copy _SEED_/$prefix";
				# Copy fasta file into new directory
				copy $fasta_file_name, VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna" or die "Could not copy $prefix.fna";

				# TODO check formats of fasta definition lines
				push @valid_prefixes, $prefix;
				$valid_prefix_found = 1;
				$parse_length = 1;
			} else {
				VH_helpers::log("\tProdigal returned an error: $?. Skipping $prefix.");
				# Prodigal creates the gff file *before* checking to see if the fasta file is valid.
				#  So if it errors out, delete the empty gff file it creates.
				unlink $gff_file_name;
				# TODO save prodigal output to file
			}
			# TODO if cancelled halfway through prodigal run, empty/invalid gff file results
		}
		$genomes{$prefix}{'name'} = $prefix;
		$genomes{$prefix}{'version'} = "";
		$genomes{$prefix}{'format'} = 'fasta';
		$genomes{$prefix}{'genes'} = 0;
		
	} else {
		# No workable files found.
		VH_helpers::log("\tNo genebank or fasta files found for $prefix. Skipping.");
	}
	if($valid_prefix_found == 1){
		# read fasta file to get sequence borders
		my $fasta_file_name = VICSIN::param("output_path")."/".CONVERTED_INPUT_DIR."/$prefix.fna";
		open(my $fastafh, '<', $fasta_file_name) or die "Could not open fasta file: ".$fasta_file_name."\n";
		my $sequence = "";
		my $seqlength = 0;
		my $seqstart = 1;
		while (my $row = <$fastafh>) {
			chomp $row;
			if(substr($row,0,1) eq '>'){
				if($sequence ne ""){
					# Save previous sequence
					$contigs{$prefix}{$sequence} = {'start'=>$seqstart,'end'=>$seqstart+$seqlength-1,'length'=>$seqlength};
					$seqstart = $seqstart + $seqlength;
					$seqlength = 0;
				}
				$sequence = substr($row,1);
			} else {
				$seqlength += length($row);
			}
		}
		# Save last sequence
		$contigs{$prefix}{$sequence} = {'start'=>$seqstart,'end'=>$seqstart+$seqlength-1,'length'=>$seqlength};
		if($parse_length == 1){
			$genomes{$prefix}{'length'} = $seqstart+$seqlength;
		}
		$genomes{$prefix}{'scaffolds'} = scalar(keys %{$contigs{$prefix}});
	}
}

### STEP 2B. Read Masking File ###
my %masks;
if(VICSIN::param('masking_file') ne '' and -f VICSIN::param('masking_file')){
	VH_helpers::log("Masking File Found. Parsing... ",1);
	open(my $masking_fh, '<', VICSIN::param('masking_file'));
	while(my $row = <$masking_fh>){
		chomp $row;
		if($row =~ m/^(.+?)\s+(\d+?)\s+(\d+)/m){
			my $prefix = $1;
			my $start = $2;
			my $end = $3;
			if(exists $contigs{$prefix}){
				foreach my $sequence ( keys %{ $contigs{$prefix} } ){
					if($end > $contigs{$prefix}{$sequence}{'start'} and $start < $contigs{$prefix}{$sequence}{'end'}){
						my $maskstart = $start;
						my $maskend = $end;
						if ( $start < $contigs{$prefix}{$sequence}{'start'} ){
							$maskstart = $contigs{$prefix}{$sequence}{'start'};
						}
						if ( $end > $contigs{$prefix}{$sequence}{'end'} ){
							$maskend = $contigs{$prefix}{$sequence}{'end'};
						}
						push @{$masks{$prefix}{$sequence}}, {'start'=>$maskstart, 'end'=>$maskend};
					}
				}
			}
		}
	}
}

if(VICSIN::param("stop") eq "input"){
	VH_helpers::log("Halting as requested after converting input files");
	exit;
}

##### STEP 3. Run Subprograms #####
### STEP 3A. Run VirSorter
if(index(VICSIN::param("skip"), 'virsorter') == -1){
	VH_VirSorter::run(\@valid_prefixes);
} else {
	VH_helpers::log("Skipping VirSorter as requested");
}

### STEP 3B. Run PhiSpy
if(index(VICSIN::param("skip"), 'phispy') == -1){
	VH_PhiSpy::run(\@valid_prefixes);
} else {
	VH_helpers::log("Skipping PhiSpy as requested");
}

### STEP 3C. Run CRISPR (optional)
my $ran_crispr = 0;
if(index(VICSIN::param("skip"), 'crispr') == -1){
	if (VICSIN::param("spacer_fasta_file") eq ""){
		VH_helpers::log("No spacer fasta file given. Skipping CRISPR match.",1);
	} elsif (! -f VICSIN::param("spacer_fasta_file")) {
		VH_helpers::log("Spacer file not found. Skipping CRISPR match.");
	} else {
		$ran_crispr = 1;
		VH_CRISPR::run(\@valid_prefixes);
	}
} else {
	VH_helpers::log("Skipping CRISPR as requested");
}

### STEP 3D. Run Spine/Agent
if(index(VICSIN::param("skip"), 'agent') == -1){
	VH_SpineAgent::run(\@valid_prefixes);
} else {
	VH_helpers::log("Skipping Spine/AGEnt as requested");
}

### STEP 3E. Run Blastn (optional)
my $ran_known_types = 0;
if(index(VICSIN::param("skip"), 'blast') == -1){
	if (VICSIN::param("known_viral_types") eq ""){
		VH_helpers::log("Known viral types not given. Skipping.",1);
	} else {
		$ran_known_types = 1;
		VH_Blast::run(\@valid_prefixes);
	}
} else {
	VH_helpers::log("Skipping BLAST as requested");
}

if(VICSIN::param("stop") eq "prediction"){
	VH_helpers::log("Halting as requested after running subprograms");
	exit;
}

### STEP 4. INITIAL PROCESS FOR GENOME COORDINATES OF OVERLAPPING ELEMENTS
VH_helpers::log("Processing output");
my %predictions;
my %merged_predictions;
my %binned_predictions;

foreach my $prefix (@valid_prefixes){
	VH_helpers::log("\tProcessing $prefix...",1);
	#Process VirSorter output
	VH_helpers::log("\t\tParsing VirSorter output... ",2);
	$predictions{$prefix}{'virsorter'} = VH_VirSorter::get_predictions($prefix);
	#Process PhiSpy output
	VH_helpers::log("\t\tParsing PhiSpy output... ",2);
	$predictions{$prefix}{'phispy'} = VH_PhiSpy::get_predictions($prefix);
	#Process CRISPR output
	if($ran_crispr == 1){
		VH_helpers::log("\t\tParsing CRISPR output... ",2);
		$predictions{$prefix}{'crispr'} = VH_CRISPR::get_predictions($prefix);
	} else {
		$predictions{$prefix}{'crispr'} = {};
	}
	#Process AGEnt output
	VH_helpers::log("\t\tParsing Agent output... ",2);
	$predictions{$prefix}{'agent'} = VH_SpineAgent::get_predictions($prefix);
	#Process Known Types Blast output
	if($ran_known_types == 1){
		VH_helpers::log("\t\tParsing Known Types output... ",2);
		$predictions{$prefix}{'blast'} = VH_Blast::get_predictions($prefix);
	} else {
		$predictions{$prefix}{'blast'} = {};
	}

	#Apply masking file
	VICSIN::apply_mask(\%predictions,$prefix,\%masks);

	#Merge results
	VH_helpers::log("\t\tMerging predictions... ",2);
	my %mergeable_predictions;
	# Dump mergeable predictions into one array
	for(my $i=0; $i<scalar(@VICSIN::methods); $i++){
		if($VICSIN::methods[$i]{'extend'} == 1){ # We only want to merge Agent, Virsorter, and Blast hits
			foreach my $sequence (keys %{$predictions{$prefix}{$VICSIN::methods[$i]{'key'}}}){
				my @predictions_to_add;
				foreach my $prediction (@{$predictions{$prefix}{$VICSIN::methods[$i]{'key'}}{$sequence}}){
					$prediction->{'methods'} = [$VICSIN::methods[$i]{'abbr'}];
					push @predictions_to_add, $prediction;
				}
				push @{$mergeable_predictions{$sequence}}, @predictions_to_add;
			}
		}
	}
	
	my %merged_predictions = VICSIN::merge_predictions(\%mergeable_predictions);

	#Compare results
	VH_helpers::log("\t\tCross-referencing predictions... ",2);
	$binned_predictions{$prefix} = VICSIN::bin_predictions(\%merged_predictions,\%predictions,$prefix);
}
print "\n";

# Output intermediate results
VH_helpers::log("Saving Consensus Output...");
foreach my $prefix (@valid_prefixes){
	#Output to file
	VH_helpers::log("\tSaving $prefix... ",1);
	sub output_bin {
		my ($fh, $bin) = @_;
		foreach my $prediction (@$bin){
			if(not exists $prediction->{'masked'}){
				print $fh $prediction->{'name'}."\t".$prediction->{'sequence'}."\t".$prediction->{'methods'}."\t".$prediction->{'start'}."\t".$prediction->{'end'}."\n";
			}
		}
	}
	make_path(VICSIN::param('output_path')."/".INTERMEDIATE_OUTPUT_DIR);
	my $output_file_name = VICSIN::param('output_path')."/".INTERMEDIATE_OUTPUT_DIR."/$prefix.txt";
	open my $output_fh, '>', $output_file_name;
	# print $output_fh "# Type 1: Predicted by >2 methods\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[0]);
	# print $output_fh "# Type 2: Predicted by 2 methods\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[1]);
	# print $output_fh "# Type 3: Predicted by 1 1° method\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[2]);
	# print $output_fh "# Type 4: Predicted by 1 2° method\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[3]);
	# print $output_fh "# Type 5: Protospacer match\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[4]);
	close $output_fh;
}

if(VICSIN::param('stop') eq "consensus"){
	VH_helpers::log('Halting as requested after consensus step');
	exit;
}

### STEP 5. RE-SCREEN PREDICTIONS AGAINST OTHER GENOMES FOR MISSED ELEMENTS ###
%binned_predictions = %{VH_ReBlast::run(\@valid_prefixes,\%binned_predictions,\%masks,\%contigs)};

VH_helpers::log("Saving Output...");
foreach my $prefix (@valid_prefixes){
	#Output to file
	VH_helpers::log("\tSaving $prefix... ",1);
	sub output_bin {
		my ($fh, $bin) = @_;
		foreach my $prediction (@$bin){
			if(not exists $prediction->{'masked'}){
				print $fh $prediction->{'name'}."\t".$prediction->{'sequence'}."\t".$prediction->{'methods'}."\t".$prediction->{'start'}."\t".$prediction->{'end'}."\n";
			}
		}
	}
	make_path(VICSIN::param('output_path')."/".OUTPUT_DIR);
	my $output_file_name = VICSIN::param('output_path')."/".OUTPUT_DIR."/$prefix.txt";
	open my $output_fh, '>', $output_file_name;
	# print $output_fh "# Type 1: Predicted by >2 methods\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[0]);
	# print $output_fh "# Type 2: Predicted by 2 methods\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[1]);
	# print $output_fh "# Type 3: Predicted by 1 1° method\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[2]);
	# print $output_fh "# Type 4: Predicted by 1 2° method\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[3]);
	# print $output_fh "# Type 5: Protospacer match\n# Prediction\tSequence\tmethods\tstart\tend\n";
	output_bin($output_fh,$binned_predictions{$prefix}[4]);
	close $output_fh;
}
print "\n";

if(VICSIN::param('stop') eq "reblast"){
	VH_helpers::log('Halting as requested after ReBLAST');
	exit;
}

### STEP 6. COMPARE/CLUSTER PREDICTIONS ###
my $clusters = VH_Cluster::run(\@valid_prefixes,\%binned_predictions);

### Finally, do database insertions ###
if(VICSIN::param('use_database') eq 'true'){
	VH_Database::insert(\@valid_prefixes, \%genomes, \%contigs, \%predictions, \%binned_predictions, $clusters);
}

VH_helpers::log("Done.");