#Script developped by group 9 of the Bioinformatics course at Universitat Pompeu Fabra. Autumn-Winter of 2012.
#You are allowed to use and modify this sript.

#Script to blast all proteins in a directory against a genome



#!/usr/bin/perl -w
use strict;
use warnings;
use File::Spec;

#The programs called through this script (blast, formatdb and fastaindex) should be able to being run from the directory where this script will be executed


my $genome_path = "./Pelodiscus_sinensis/genome.fa";

my $proteins_path = "./selenoproteins/";


(my $volume, my $genome_dir, my $genome_name) = File::Spec->splitpath($genome_path);


###Format the genome if not previously done(formatdb)###
#Output: formated database of the genome stored inside the genome directory provided by the standard input

	#Put all the files in the genome path into a list

opendir my $dir_genome, "$genome_dir" or die "Cannot open genome directory: $!";
my @genome_files = readdir $dir_genome; 
closedir $dir_genome;

my $b;
my $a;
my $database_name;
	
	#Check if the genome is already formated

foreach my $file (@genome_files) {  
	if ($file =~ m/\.nhr/) {
		$file =~ m/([a-zA-Z0-9\.\s_]+)\.nhr/;
		$database_name = $1;
		$b = 1;
	}
	elsif ($file =~ m/\.fa/) {
		$file =~ m/([a-zA-Z0-9\.\s_]+)\.fa/;
		$database_name = $1;
	}
	elsif ($file =~ m/\.index/) {
		$a = 1;
	}
	
}


	#if not formated use formatdb to format it

if (!$b) {
	my $formatdb_command = "formatdb -i $genome_path -p F -n $database_name.fa";	
	system($formatdb_command);
	system("mv $database_name.* $genome_dir");
	print "The genome has been formated\n";
}

else {
	print "The genome is already formated\n";
}

if (!$a) {
	system("fastaindex $genome_path $database_name.index");  #creates an index of the genome so as to be able to run fastafetch on it later
	system("mv $database_name.index $genome_dir");
}

###BLAST the proteins against the genome(tblastn)###
#Output: blast_out.txt

	#Put all the files in the proteins directory into a list

opendir my $dir_proteins, "$proteins_path" or die "Cannot open proteins directory: $!";
my @protein_files = readdir $dir_proteins; 
closedir $dir_proteins;


my $prot_name;
my $database = "$genome_dir/$database_name\.fa";
my $blast_out;
my $blast_command; 

system("mkdir ./BLAST_OUTPUT"); #All the output files from the blast will be stored in this directory

foreach my $file (@protein_files) {
	if ($file =~ m/^[^\.]/) {
		$prot_name = $file =~ m/([a-zA-Z0-9\.\s_]+)\.aa/;
		$blast_out = "./BLAST_OUTPUT/" . $1 . "blast_out.txt";
		$blast_command = "blastall -p tblastn -i " . $proteins_path . $file . " -d $database -o $blast_out -m8";
		system($blast_command);
	}
}




