#!/bin/bash
# RNA-Seq 
# 2021-01

#SBATCH --nodes=1			
#SBATCH --ntasks-per-node=2
#SBATCH --mem=64GB
#SBATCH --job-name=BLAST_Small_OFS


# Choose between long or short: A short job will start quickly after submitting, but will be terminated after 30 minutes

#SBATCH --time=24:00:00
##SBATCH -p short

export PERL5LIB=/data/pg-molgen/software/molgentools_mirror/lib
module load BioPerl/1.7.7-GCCcore-9.3.0
module load BLAST+/2.7.1-foss-2018a
module load Python/3.6.4-foss-2018a


# parameters: Modify for your folders
	cpu=2
	SCRATCHDIR=/scratch/p276244
	molgentools=/data/pg-molgen/software/molgentools_mirror
	
	DIR_small_orfs=/data/p127804/PROJECTS/JanMaarten/Elisa/Hybrid_Assembly_Stau_fasta
	DIR_RAST_genomes=/data/p127804/PROJECTS/JanMaarten/Elisa/RAST_Annotation
	DIR_results=/data/p127804/PROJECTS/JanMaarten/Elisa/blast_small_orfs

# =====================================================================================================================================================
# ==                                                          Database setup                                                                         ==
# =====================================================================================================================================================

# only once after update;  makeblastdb -in uniprot_sprot_clean.fasta -dbtype prot
# For uniprot update see GSEA-Pro and FACOP

# =====================================================================================================================================================
# ==                                                          MAIN                                                                                   ==
# =====================================================================================================================================================

function my_blast {
	# 1. Build blast database
	makeblastdb -in $DIR_RAST_genomes/$2.faa -dbtype prot
	# 2. Blast all small orfs to the RAST genome annotation
	blastp -outfmt 6 -db $DIR_RAST_genomes/$2.faa -query $DIR_small_orfs/$1.faa -max_target_seqs 1 -num_threads $cpu -evalue 10E-09 -out $DIR_results/$1.blast
	# 3. Filter for best blast hit only
	$molgentools/blast/parse_blast_uniref_table.pl -i $DIR_results/$1.blast -o $DIR_results/$1.blast.filtered
	# 4. Filter for 100% perfect match over full length of the protein
	python3 $molgentools/blast/Blast_Table_Parser.py -query $DIR_results/$1.blast.filtered -out $DIR_results/$1.blast
	# 5. Make new fasta file for small orfs  ==> $1.imperfect.match.faa
		# a. make tab file of all small orfs
		$molgentools/genomics/g2d_conversions.pl -s $DIR_small_orfs -l $1.faa -r Cvrt_fasta2tab -o $1.tab	
		# b. algin $1.blast.imperfect.match   to  small orfs
		$molgentools/genomics/g2d_align_tables.pl -s $DIR_results -t $DIR_results/$1.blast.imperfect.match -k $DIR_small_orfs/$1.tab -r Align_firstcolumn -o $1.imperfect.match.tab
		# c. get columns 1 and 4
		awk '{print $1"|"$2,$4}' OFS="\t" $DIR_results/$1.imperfect.match.tab > $DIR_results/$1.imperfect.match.14.tab
		# d. convert back to fasta
		$molgentools/genomics/g2d_conversions.pl -s $DIR_results -l $1.imperfect.match.14.tab -r Cvrt_tab2fasta -o $1.imperfect.match.faa	
}


mkdir $DIR_results

# call the BLAST
my_blast 'hybrid_171_131.30.101.sOFS'    'Staphylococcus_aureus_hybrid_171_131'
my_blast 'hybrid_LL58.30.101.sOFS'       'Staphylococcus_aureus_hybrid_LL58'
my_blast 'hybrid_LL059.30.101.sOFS'      'Staphylococcus_aureus_hybrid_LL059'
my_blast 'hybrid_LL327_S11.30.101.sOFS'  'Staphylococcus_aureus_hybrid_LL327_S11'
my_blast 'hybrid_LL071.30.101.sOFS'      'Staphylococcus_aureus_hybrid_LL071'
my_blast 'hybrid_173_102.30.101.sOFS'    'Staphylococcus_aureus_hybrid_173_102'