#!/bin/bash
# Apr 3, 2018
# use sequence data after pear (merged reads) generated by Cascabel (our
# NIOZ amplicon data analysis pipeline, see config_Cascabel.yaml,
# https://github.com/AlejandroAb/CASCABEL and https://doi.org/10.1101/809384),
# pick OTUs with trie (representative: most abundant),
# assign taxonomy with blast using customized query coverage and identity cutoffs.
#######################

# (path to) fasta file containing the paired sequence reads (forward and reverse merged using the overlap)
fasta=seqs_fw_rev_filtered.fasta
# mapping file with sequence identifiers and taxonomy path of sequences in the reference database
taxDB=12Sfish_20180320_id_to_tax.map

# align reads to ref db w/ mothur v.1.34.4. mothur does not tolerate dashes in file path, 
# therefore, link the file to the current dir, if it is not already there
ln -s $fasta pear_seqs_fw_rev_filtered.fasta

# mothur also requires the template (reference database) to be aligned
# link it to current working dir if not already there
ln -s Databases/12Sfish_20180320/12Sfish_20180320_aligned.fa .

# but blast does not like gaps in the sequences, so also generate a degapped reference file
# using 'degapseq' fr/ EMBOSS:6.6.0.0. 
degapseq 12Sfish_20180320_aligned.fa 12Sfish_20180320_aligned_degapped.fa
# convert from multi-line to one-line fasta
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n }' 12Sfish_20180320_aligned_degapped.fa > 12Sfish_20180320_aligned_degapped_oneLine.fa


################################################
# align sequence reads to the reference database
# run mothur in non-interactive mode: "# cmd "
mothur '#align.seqs(fasta=pear_seqs_fw_rev_filtered.fasta, reference=12Sfish_20180320_aligned.fa)'

# hard cut 19nt (the primer) from the front with cutadapt
# we need to cut first 19 nt because the alignment has 1 insertion.
# The file needs to end in 'fasta', therefore rename. 
# we used cutadapt version 1.15
mv pear_seqs_fw_rev_filtered.align pear_seqs_fw_rev_filtered.fasta
cutadapt -u 19 -o trimmed.fasta pear_seqs_fw_rev_filtered.fasta

# remove last 18 nt (remove reverse primer)
cutadapt -u -18 -o pear_seqs_fw_rev_filtered_noPrimers.fasta trimmed.fasta

# degap the sequences, then remove too short sequences:
# using 'degapseq' fr/ EMBOSS:6.6.0.0. 
degapseq  pear_seqs_fw_rev_filtered_noPrimers.fasta  pear_seqs_fw_rev_filtered_noPrimers_degapped.fasta

# The output is a multi-line, not a one-line fasta, so change to one-line fasta
awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n }' pear_seqs_fw_rev_filtered_noPrimers_degapped.fasta > pear_seqs_fw_rev_filtered_noPrimers_degapped_oneLine.fasta

# count the number of sequences
grep '>' pear_seqs_fw_rev_filtered_noPrimers_degapped_oneLine.fasta | wc -l

# remove very short and very long sequences. Actually, none are longer than 200,
# so we might as well remove only the short ones (<= 90 nt).
awk '!/^>/ { next } { getline seq } length(seq) >= 90 && length(seq) <=200 { print $0 "\n" seq }'  pear_seqs_fw_rev_filtered_noPrimers_degapped_oneLine.fasta > seqs_filtered.fasta

# pick OTUs using the 'trie' method from qiime (qiime v1)
echo 'now picking OTUs ...'
pick_otus.py -m trie -i seqs_filtered.fasta -o OTUs
echo 'Number of OTUs: '
wc -l OTUs/seqs_filtered_otus.txt

# pick representative sequences. 
pick_rep_set.py -m most_abundant -i OTUs/seqs_filtered_otus.txt -f seqs_filtered.fasta -o OTUs/pear_rep_set.fasta

# count the number of representative sequences
count_seqs.py -i OTUs/pear_rep_set.fasta >> seqCounts.txt

# run blast against our degapped file with reference sequences
blastn -query OTUs/pear_rep_set.fasta -task megablast -subject 12Sfish_20180320_aligned_degapped_oneLine.fa -out OTUs/pear_rep_set_noPrimers_blast_qcov98.out -evalue 0.01 -outfmt '6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs' -num_alignments 10 -perc_identity 98 -qcov_hsp_perc 98

# reorder blast output to get the hit with highest sequence identity first (did not change anything, but blast output 
# needs to be reformatted anyways)
# using an R script. Tested with R version 3.4.3
echo "start R"
R
source("utils.R")
res <- reorderByIdentity(taxMapFile="12Sfish_20180320_id_to_tax.map", blastOutFile="OTUs/rep_set_noPrimers_blast_qcov98.out",returnObjects=F, includeDate=F)
#the script removes blast***.out and substitutes with "_tax_ass_identity.txt"
q('no')
echo "stopped with R,"
echo "continue with bash script"

# generate an OTU table
make_otu_table.py -i OTUs/seqs_filtered_otus.txt -t OTUs/pear_rep_set_noPrimers_tax_ass_identity.txt -o OTUs/pear_rep_set_noPrimers_tax_ass_identity_otuTable.biom

#filter out overall singletons
filter_otus_from_otu_table.py -i OTUs/pear_rep_set_noPrimers_tax_ass_identity_otuTable.biom -o OTUs/pear_rep_set_noPrimers_tax_ass_identity_otuTable_noSingletons.biom -n 2

#convert otu table from biom to tsv
biom convert -i OTUs/pear_rep_set_noPrimers_tax_ass_identity_otuTable_noSingletons.biom -o OTUs/pear_rep_set_noPrimers_tax_ass_identity_otuTable_noSingletons.tsv --table-type 'OTU table' --header-key taxonomy --to-tsv