BASH script for SelfBLAST a genome
By Jit 2643 days ago
#!/bin/bash
#self BLAST a genome -- Expecting you have blast and samtools installed in your system
#Author: Jitendra Narayan
#USAGE: ./selfBlast.sh extract <chrName>
#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 :)"
