#!/bin/bash
#
# Name:         blastn_parallel.sh
#
#
# Purpose:      This program uses blastn to match the query sequences
#		against the database and generates the blast output file  
#		in the current working directory. 
#               The filename of blast output file is the original fasta filename 
#		plus "_blast.out".
#
#		In this program, I have used GNU parallel http://www.gnu.org/software/parallel/
#		which is a shell tool for executing jobs in parallel.
#		The script divides the fasta file into $totalChunks and then run blastn in parallel
#		on each chunk by piping it to GNU parallel program. Please modify $totalChunks to 
#		reflect the number of cores available for processing. 
#		
# Author(s):    Umer Zeeshan Ijaz
#
# Address:      University of Glasgow
#               School of Engineering   
#               Glasgow
#               G12 8LT
#               Umer.Ijaz@glasgow.ac.uk
#
# Version:      0.1 

#store the input arguments
args=("$@")

# = Error checks =================== #
if [ $# -ne 2 ]; then
	echo "Usage: `basename $0` <fastafile> <blastdatabasename>" 
	exit 0
fi

if [ ! -f ${args[0]} ]; then
	echo "File doesn't exist. Please try again!"
	exit 0
fi

if [ `echo "${args[0]}"| sed 's/.*\.//'` != "fasta" ]; then
	echo "File doesn't have .fasta extension. Please try again!"
	exit 0
fi
# ================================== #

# = Enable FP support ============== #
# By default, there is limited capability in bash to handle floating point
# operations. In this script bc is used to calculate the floating point operations.
# $float_scale parameter specifies the precision of the floating point.
# Reference: http://www.linuxjournal.com/content/floating-point-math-bash

float_scale=5
# Evaluate a floating point number expression.

function float_eval()
{
    local stat=0
    local result=0.0
    if [[ $# -gt 0 ]]; then
        result=$(echo "scale=$float_scale; $*" | bc -q 2>/dev/null)
        stat=$?
        if [[ $stat -eq 0  &&  -z "$result" ]]; then stat=1; fi
    fi
    echo $result
    return $stat
}


# Evaluate a floating point number conditional expression.

function float_cond()
{
    local cond=0
    if [[ $# -gt 0 ]]; then
        cond=$(echo "$*" | bc -q 2>/dev/null)
        if [[ -z "$cond" ]]; then cond=0; fi
        if [[ "$cond" != 0  &&  "$cond" != 1 ]]; then cond=0; fi
    fi
    local stat=$((cond == 0))
    return $stat
}

function ceil () {
  echo "define ceil (x) {if (x<0) {return x/1} \
        else {if (scale(x)==0) {return x} \
        else {return x/1 + 1 }}} ; ceil($1)" | bc;
 }

# =/Enable FP support ============== #

readFile=${args[0]}
blastDatabaseName=${args[1]}
blastFile=`echo "$(basename $readFile)" | cut -d'.' -f1`'_blast.out'

# =Get the file size in KB========== #
sizeFileBytes=$(du -b ${readFile} | sed 's/\([0-9]*\)\(.*\)/\1/')
totalChunks=45
sizeChunks=$(ceil $(float_eval "$sizeFileBytes / ($totalChunks * 1024)"))
sizeChunksString="${sizeChunks}k"
# =/Get the file size in KB========= #

startTime=`date +%s`

cat $readFile | parallel --block $sizeChunksString --recstart '>' --pipe blastn -perc_identity 90 -evalue 0.00001 -dust no -num_threads 16 -outfmt 6 -max_target_seqs 1 -db $blastDatabaseName -query - > $blastFile

echo blastn took $(expr `date +%s` - $startTime) seconds for $readFile ...