set -ue

# Load the conda commands.
source /home/www/miniconda3/etc/profile.d/conda.sh

# Input file that contains the accession numbers.
ACC=/export/www/biostar-central/export/media/projects/cookbook/data-275/Fish_marker_sequences

# How many reads to simulate.
N=1000

# Length of reads.
L=100

# Folder with results.
OUT=output

# Create directories to store files.
mkdir -p input code $OUT

# The final result file
ACCURACY=classification-qiime2-accuracy.txt

# The FASTA file storing the sequence for accession numbers.
FASTA=input/reference.fa

# The FASTQ file storing the reads simulated from the reference.
READS=input/reads.fq

# Fasta sequence in QIIME2 compatible format.
FASTA_QZA=$OUT/reference.qza

# Sequence reads in QIIME2 compatible format.
READS_QZA=$OUT/reads.qza

# Reads manifest file (This is like a sample sheet).
MANIFEST=input/reads_manifest.txt

# Qiime2 classifier (database).
CLASSIFIER=$OUT/trained_reference.qza

# Table with de-replicated information.
TABLE_QZA=$OUT/feature_table.qza

# File with representative sequences after dereplication.
REPS_QZA=$OUT/represetative_seqs.qza

# File with taxonomic classification in qiime format.
CLASSIFIED_QZA=$OUT/classified.qza

# The accesion taxonomy file.
TAXONOMY=input/taxonomy.txt

# File with taxa lineage.
LINEAGE=input/taxa_lineage.tsv

# QIIME2 compatable taxa lineage file.
LINEAGE_QZA=$OUT/taxa_lineage.qza

# File with acccessions without version.
ACCESSIONS=input/accessions.txt


# Remove version from accessions if present and copy to input folder.
cat $ACC | sed -E 's/\.[0-9]+//' >$ACCESSIONS

# Use the taxonomy specific files to build the custom database.
TAXDIR=/export/refs/taxonomy

#
# Download taxonomy specific files and create taxonomy lineage database using the script below.
# This operation only needs to be done only once for the entire website.
# It will create 'taxon_db' database in $TAXDIR.
# python -m recipes.code.taxon_lineage --taxdir $TAXDIR
#
DBPATH=$TAXDIR/lineage_db

# Create taxonomy file for the input accessions.
python -m recipes.code.get_taxa_lineage --dbpath $DBPATH --accessions $ACCESSIONS --outfile $LINEAGE

# Identify the taxonomy for each accession number.
epost -db nuccore -input $ACC  | esummary | xtract -pattern DocumentSummary -element Caption,TaxId,Title > $TAXONOMY

# Get sequences from NCBI.
epost -db nuccore -input $ACC | efetch -format fasta > $FASTA

# Remove version number from fasta header so that it matches with the taxonomy lineage file.
cat $FASTA | sed -E 's/\.[0-9]+//' >tmp.fa
mv tmp.fa $FASTA

# Install the plac command parser.
pip install plac -q

# The location of the code that simulates the reads.
URL1=https://raw.githubusercontent.com/biostars/biocode/master/scripts/fasta/simulate.py

# Get the code.
curl $URL1 > code/simulate.py

# Generate the simulated reads.
python code/simulate.py --fname $FASTA --count $N > $READS

# Switch to the qiime enviorment.
set +u
conda activate qiime2
set -ue

# Convert reference sequence to qiime compatible format.
qiime tools import --type 'FeatureData[Sequence]' --input-path $FASTA --output-path $FASTA_QZA 

# Convert taxonomy lineage to qiime compatible format.
qiime tools import --type 'FeatureData[Taxonomy]' --input-path $LINEAGE --output-path $LINEAGE_QZA

# Train classifer.
qiime feature-classifier fit-classifier-naive-bayes --i-reference-reads $FASTA_QZA --i-reference-taxonomy $LINEAGE_QZA --o-classifier $CLASSIFIER

# Create a manifest file for the input reads.
echo -e "sample-id\tabsolute-filepath" >$MANIFEST
echo -e "sample1\t$PWD/$READS" >>$MANIFEST

# Convert sequence reads to qiime compatible format using reads manifest file.
qiime tools import --type 'SampleData[SequencesWithQuality]' --input-path $MANIFEST --input-format SingleEndFastqManifestPhred33V2 --output-path $READS_QZA

# Cluster and dereplicate input reads to create representative sequences.
qiime vsearch dereplicate-sequences --i-sequences $READS_QZA --o-dereplicated-table $TABLE_QZA --o-dereplicated-sequences $REPS_QZA

# Create represetative sequence fasta.
qiime tools export --input-path $REPS_QZA --output-path $OUT

# Give meanigful name to representative sequences.fasta
mv $OUT/dna-sequences.fasta $OUT/representative_seqs.fa

# Classify input data with previously trained classifier.
qiime feature-classifier classify-sklearn --i-classifier $CLASSIFIER --i-reads $REPS_QZA --o-classification $CLASSIFIED_QZA --p-n-jobs 2

# Export feature table to biom format.
qiime tools export --input-path $TABLE_QZA --output-path $OUT

# Convert feature table from biom.
biom convert -i $OUT/feature-table.biom  -o $OUT/feature_table.tsv --to-tsv 

# Export classified qza file to tsv format. 
qiime tools export --input-path $CLASSIFIED_QZA --output-path $OUT

# TEST Rename file
mv $OUT/taxonomy.tsv $OUT/classified.tsv

# Reformat the header to make it work with biom.
sed -i $'1s/^.*$/#OTUID\ttaxonomy\tconfidence/' $OUT/classified.tsv

# Add metadata to biom table.
biom add-metadata -i $OUT/feature-table.biom -o $OUT/feature-table-taxonomy.biom --observation-metadata-fp $OUT/classified.tsv --sc-separated taxonomy

# Create the classification file.
biom convert -i $OUT/feature-table-taxonomy.biom -o $OUT/feature-table-taxonomy.tsv --to-tsv --header-key taxonomy

# Collapse at species level.
qiime taxa collapse  --i-table $TABLE_QZA --i-taxonomy $CLASSIFIED_QZA --p-level 7 --o-collapsed-table $OUT/species_table.qza

# Export  species table to biom format.
qiime tools export --output-path ./ --input-path $OUT/species_table.qza

# Move species biom table to TMP folder.
mv feature-table.biom  $OUT/feature-table-species.biom 

# Convert species biom table to tsv formt.
biom convert -i $OUT/feature-table-species.biom -o classification_qiime2.txt --to-tsv

# Install a required librar
pip install plac --upgrade

# Generate the 
# The location of the code that generates the validation.
URL2=https://raw.githubusercontent.com/biostars/biocode/master/scripts/classify/validate.py

# Get the code.
curl $URL2 > code/validate.py

# Run the validator
python code/validate.py -f classification_qiime2.txt -t $TAXONOMY -c $N > $ACCURACY