#!/bin/sh
#
#   script: (1-2) This script takes as input an rRNA database in the format 
#           of a fasta file and an xml file, both exported from the 
#           ARB package for some complete rRNA database (ex. SILVA).
#           We use the xml description tags to filter out rRNA
#           which are not 16S or 23S rRNA. 
#           ARB source http://www.arb-home.de/
#
#			(3) The resulting (only 16S or 23S)
#           fasta file is further filtered by PRINSEQ to remove
#           long sequences and sequences with more than 1% occurrences
#           of the ambiguous letter 'N'. Long sequences are defined to be
#           >1600 nucleotides for 16S rRNA and >5000 nucleotides for 
#           23S rRNA. This aids to remove sequences labeled as `complete 
#           genome', which can include a complete rRNA molecule along
#           with other coding or non-coding nucleic acids appended to 
#           either end of the molecule (mRNA, tRNA ..). Since only rRNA
#           sequences are masked out from the NCBI complete bacterial 
#           genomes, the resulting non-rRNA reads can match to an rRNA
#           database if sequences labeled as `complete genome' are not 
#           filtered out. 
#           PRINSEQ source http://sourceforge.net/projects/prinseq/files/standalone/prinseq-lite-0.17.4.tar.gz
#
#			(4) Finally, we apply UCLUST to the filtered database to create
#           a representative rRNA database with x% id.
#	        UCLUST source http://www.drive5.com/usearch/
#
#   contact: evguenia.kopylova@lifl.fr
#
#   date: July 2012
#


# path to ARB-exported rRNA sequences xml file
SILVAxml=""
# path to ARB-exported rRNA sequences fasta file
SILVAfasta=""
# path to PRINSEQ
PRINSEQ="prinseq-lite.pl"
# path to UCLUST
UCLUST="usearch5.2.32_i86linux32"
# UCLUST -id value
ID=98

# Step (1) find non-23s sequences

cat $SILVAxml | grep "5S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna5s.txt
cat $SILVAxml | grep "5.8S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna5.8s.txt
cat $SILVAxml | grep "16S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna16s.txt
cat $SILVAxml | grep "18S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna18s.txt
cat $SILVAxml | grep "26S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna26s.txt
cat $SILVAxml | grep "28S" -B 45 | sed -n 's/<species name="\(.*\)">/\1/p' > rrna28s.txt

cat rrna*.txt > totalnon23srrna_temp.txt
sort totalnon23srrna_temp.txt | uniq -c | awk '{print $2;} ' > totalnon23srrna.txt
rm rrna*.txt
rm totalnon23srrna_temp.txt

echo "Found all non-23S rRNA."
echo `wc -l totalnon23srrna.txt`

# Step (2) remove non-23s sequences

# put fasta file into 1 line per entry
awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < $SILVAfasta \
	| egrep -v '^$' > ${SILVAfasta%.fasta}-1line.fasta

while read not23s
do
	sed -i "/$not23s /d" ${SILVAfasta%.fasta}-1line.fasta
done < "totalnon23srrna.txt"
rm totalnon23srrna.txt

# put back into 2 lines per entry
tr "\t" "\n" < ${SILVAfasta%.fasta}-1line.fasta > ${SILVAfasta%.fasta}-no23s.fasta

echo "Removed all non-23S rRNA."

# Step (3) PRINSEQ (maximum length 5000, maximum ambigious N's 1%)
$PRINSEQ -fasta ${SILVAfasta%.fasta}-no23s.fasta \
   -max_len 5000 -ns_max_p 1 -verbose -out_good ${SILVAfasta%.fasta}-no23s-prinseq -out_bad null

# Step (4) UCLUST
$UCLUST -sort ${SILVAfasta%.fasta}-no23s-prinseq.fasta -output ${SILVAfasta%.fasta}-no23s-prinseq-sorted.fasta
$UCLUST -cluster ${SILVAfasta%.fasta}-no23s-prinseq-sorted.fasta -seedsout ${SILVAfasta%.fasta}-database-id$ID.fasta \
	-id 0.$ID

echo "Results:"

echo `grep -c '>' ${SILVAfasta%.fasta}-database-id$ID.fasta`

echo "Done."




