#!/bin/sh
#
#   script: (1) Remove the circular reads created by MetaSim, these
#           are defined to have more than one KEY in the SOURCES
#           brackets given in the tag of the generated read.
#           MetaSim source http://ab.inf.uni-tuebingen.de/software/metasim/
#
#			(2) Further filter the linear reads to remove sequences
#           with <200nt and >1% of ambiguous N's. 
#           PRINSEQ source http://sourceforge.net/projects/prinseq/files/standalone/prinseq-lite-0.17.4.tar.gz
#
#           (3) Randomly choose 300,000 Roche 454 reads
#
#   note: For Illumina reads, remove "-min_len 200" in Step (2) and change "head -n 300000" to
#         "head -n 1000000" in Step (3).
#
#   contact: evguenia.kopylova@lifl.fr
#
#   date: July 2012
#

# path to 454 reads generated by MetaSim (suffix must be .fna)
METASIMREADS="bac16s-85id-454.3357f545.fna"
# path to PRINSEQ
PRINSEQ="prinseq-lite.pl"


# Step (1) Remove the circular reads
awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' \
		< $METASIMREADS \
		| egrep -v '^$' \
		| awk -F "KEY=" '{i=NF-1; if(i==1) print $0}' \
		| tr "\t" "\n" > "${METASIMREADS%.fna}-nocircular-reads.fasta"

# Step (2) PRINSEQ on linear reads to remove sequences with <200nt and >1% ambiguous N's
perl $PRINSEQ -fasta "${METASIMREADS%.fna}-nocircular-reads.fasta" -min_len 200 \
    ns_max_p 1 -out_good "${METASIMREADS%.fna}-prinseq" -out_bad null -verbose

rm "${METASIMREADS%.fna}-nocircular-reads.fasta"

# Step (3) Randomly choose 300,000 Roche 454 linear-PRINSEQ filtered reads
awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' \
		< "${METASIMREADS%.fna}-prinseq.fasta" | egrep -v '^$' | shuf \
        | head -n 300000 \
        | tr "\t" "\n" > "${METASIMREADS%.fna}-final.fna"

rm "${METASIMREADS%.fna}-prinseq.fasta"

echo "Done."













