BASH script for SelfBLAST a genome

Jit — Mon, 30 Jan 2017 09:31:33 -0600

#!/bin/bash

#self BLAST a genome -- Expecting you have blast and samtools installed in your system
#Author: Jitendra Narayan
#USAGE: ./selfBlast.sh extract 
#USAGE: ./selfBlast.sh all

#Common settings 
FASTAFILE=MergedContigs.fasta
MYDB=myDB
OUTFILE=seeRES
THREAD=20
SEQ=""

echo "User $USER provided $# arguments, Detail of the arguments: $@"

if [ -f $MYDB.nhr ]
then
  echo "BLAST database for MergedContigs.fasta genome exists"
else
  echo "Thanks for testing this script $USER; Me creating creating blastDB named $MYDB for you";
  makeblastdb -in $FASTAFILE -parse_seqids -dbtype nucl -out $MYDB
fi

if [ $1 = "extract" ]
then
  echo "Extracting the sequence $2 for you from $FASTAFILE -- MAKE SURE U HAVE ADDED CORRECT NAME"
  samtools faidx MergedContigs.fasta
  samtools faidx MergedContigs.fasta $2 > $2.fa
  SEQ=$2.fa
elif [ $1 = "all" ]
then
  echo "You want entire sequence to blast"
  SEQ=$FASTAFILE
else
  echo "Something went wrong $USER - Contact jitendra"
fi

echo "Doing alignments -- BLASting";
blastn -task megablast -query $SEQ -db $MYDB -evalue 1e-5 -num_threads $THREAD -max_target_seqs 1 -outfmt '6 qseqid staxid qstart qend sseqid sstart send evalue length frames qcovs' -out $OUTFILE;

echo "DONE successfully :)"

BOL: BASH script for SelfBLAST a genome

BASH script for SelfBLAST a genome