#!/bin/sh
#
#   script: Here we run BLASTN on Set 5: 85% id. bacteria & archaea (7659 rRNA);
#           The same script can be launched for all database sets and
#           their respective reads. 
#
#			We parse the BLASTN results using blastn-parser.pl
#
#           Commented lines are for the two real metatranscriptomic read sets,
#           uncomment them for results. Make sure to create a 16S indexed library
#           for SRR106861, and a 23S indexed library for SRR013513.
#
#			Step (1) - format database
#			Step (2) - 300,000 Roche rRNA reads
#			Step (3) - 1,000,000 Illumina rRNA reads
#			Step (4) - 1,000,000 Roche non-rRNA reads
#			Step (5) - 1,000,000 Illumina non-rRNA
#			Step (6) - SRR106861 filtered reads (for 16S testing)
# 			Step (7) - SRR013513 filtered reads (for 23S testing)
#
#   note:   all fasta files should have .fasta suffix
#
#           BLASTN source ftp://ftp.ncbi.nlm.nih.gov/blast/executables/release/2.2.25/
#
#   contact: evguenia.kopylova@lifl.fr
#
#   date: July 2012
#

# path to 300,000 Roche 454 rRNA reads (produced by clean-reads.sh)
ROCHERRNA="set5-roche-rrna-reads.fasta"

# path to 1,000,000 Illumina rRNA reads (produced by clean-reads.sh)
ILLUMINARRNA="set5-illumina-rrna-reads.fasta"

# path to 1,000,000 Roche 454 non-rRNA reads
ROCHENONRRNA="roche-non-rrna-reads.fasta"

# path to 1,000,000 Illumina non-rRNA reads
ILLUMINANONRRNA="illumina-non-rrna-reads.fasta"

# path to metatranscriptomic filtered set SRR106861
METAT1="SRR106861-filtered.fasta"

# path to metatranscriptomic filtered set SRR013513
METAT2="SRR013513-filtered.fasta"

# path to rRNA database file
DBSPATH="set5-database.fasta"

# indexed representative database name created by formatdb (files .nhr, .nin, .nsq)
DBS="bac16S-85id"


# --------------------------------------
# Step (1) Index database

formatdb -i $DBSPATH -p F -n $DBS

# --------------------------------------
# Step (2) 300,000 Roche rRNA reads

echo "300,000 Roche rRNA reads"
time blastall -p blastn -m 8 -a 1 -d $DBS -i $ROCHERRNA -o ${ROCHERRNA%.fasta}-blastn.tabular

# Parse the tabular format
echo "Parsing ..."
time perl blastparser.pl ${ROCHERRNA%.fasta}-blastn.tabular > ${ROCHERRNA%.fasta}-blastn.tags

# Calculate unique reads
echo "Calculate unique reads ..."
awk '{print $1; }' ${ROCHERRNA%.fasta}-blastn.tags | sort - | uniq -c > ${ROCHERRNA%.fasta}-blastn-unique.tags
echo `wc -l ${ROCHERRNA%.fasta}-blastn-unique.tags`

# -------------------------------------- 
# Step (3) 1,000,000 Illumina rRNA reads

echo "1,000,000 Illumina rRNA reads"
time blastall -p blastn -m 8 -a 1 -d $DBS -i $ILLUMINARRNA -o ${ILLUMINARRNA%.fasta}-blastn.tabular

# Parse the tabular format
echo "Parsing ..."
time perl blastparser.pl ${ILLUMINARRNA%.fasta}-blastn.tabular > ${ILLUMINARRNA%.fasta}-blastn.tags

# Calculate unique reads
echo "Calculate unique reads ..."
awk '{print $1; }' ${ILLUMINARRNA%.fasta}-blastn.tags | sort - | uniq -c > ${ILLUMINARRNA%.fasta}-blastn-unique.tags
echo `wc -l ${ILLUMINARRNA%.fasta}-blastn-unique.tags`


# ---------------------------------------
# Step (4) 1,000,000 Roche non-rRNA reads

echo "1,000,000 Roche non-rRNA reads"
time blastall -p blastn -m 8 -a 1 -d $DBS -i $ROCHENONRRNA -o ${ROCHENONRRNA%.fasta}-blastn.tabular

# Parse the tabular format
echo "Parsing ..."
time perl blastparser.pl ${ROCHENONRRNA%.fasta}-blastn.tabular > ${ROCHENONRRNA%.fasta}-blastn.tags

# Calculate unique reads
echo "Calculate unique reads ..."
awk '{print $1; }' ${ROCHENONRRNA%.fasta}-blastn.tags | sort - | uniq -c > ${ROCHENONRRNA%.fasta}-blastn-unique.tags
echo `wc -l ${ROCHENONRRNA%.fasta}-blastn-unique.tags`


# ---------------------------------------
# Step (5) 1,000,000 Illumina non-rRNA

echo "1,000,000 Illumina non-rRNA reads"
time blastall -p blastn -m 8 -a 1 -d $DBS -i $ILLUMINANONRRNA -o ${ILLUMINANONRRNA%.fasta}-blastn.tabular

# Parse the tabular format
echo "Parsing ..."
time perl blastparser.pl ${ILLUMINANONRRNA%.fasta}-blastn.tabular > ${ILLUMINANONRRNA%.fasta}-blastn.tags

# Calculate unique reads
echo "Calculate unique reads ..."
awk '{print $1; }' ${ILLUMINANONRRNA%.fasta}-blastn.tags | sort - | uniq -c > ${ILLUMINANONRRNA%.fasta}-blastn-unique.tags
echo `wc -l ${ILLUMINANONRRNA%.fasta}-blastn-unique.tags`

# ---------------------------------------------------
# Step (6) SRR106861 filtered reads (for 16S testing)

#echo "SRR106861 filtered reads"
#time blastall -p blastn -m 8 -a 1 -d $DBS -i $METAT1 -o ${METAT1%.fasta}-blastn.tabular

# Parse the tabular format
#echo "Parsing ..."
#time perl blastparser.pl ${METAT1%.fasta}-blastn.tabular > ${METAT1%.fasta}-blastn.tags

# Calculate unique reads
#echo "Calculate unique reads ..."
#awk '{print $1; }' ${METAT1%.fasta}-blastn.tags | sort - | uniq -c > ${METAT1%.fasta}-blastn-unique.tags
#echo `wc -l ${METAT1%.fasta}-blastn-unique.tags`

# ---------------------------------------------------
# Step (7) SRR013513 filtered reads (for 23S testing)

#echo "SRR013513 filtered reads"
#time blastall -p blastn -m 8 -a 1 -d $DBS -i $METAT2 -o ${METAT2%.fasta}-blastn.tabular

# Parse the tabular format
#echo "Parsing ..."
#time perl blastparser.pl ${METAT2%.fasta}-blastn.tabular > ${METAT2%.fasta}-blastn.tags

# Calculate unique reads
#echo "Calculate unique reads ..."
#awk '{print $1; }' ${METAT2%.fasta}-blastn.tags | sort - | uniq -c > ${METAT2%.fasta}-blastn-unique.tags
#echo `wc -l ${METAT2%.fasta}-blastn-unique.tags`

echo "done."










