#!/bin/bash

if [ x$usearch == x ] ; then
	echo Must set \$usearch >> /dev/stderr
	exit 1
fi

rm -rf ../out
mkdir -p ../out
cd ../out

# Assign reads to samples (demultiplex)
$u64 -fastx_demux ../data/reads.fq -barcodes ../data/barcodes.fa -sample_delim . \
  -fastqout reads_demux.fq

# Strip primer (CCGTCAATTCMTTTRAGT, len 18)
$u64 -fastx_truncate reads_demux.fq -stripleft 18 -fastqout reads_stripped.fq

# Quality vs. length report to choose trim length
$usearch -fastq_eestats2 reads_stripped.fq -output qual.txt

# Quality filter & trim to length 350 (chosen after review of qual report)
$usearch -fastq_filter reads_stripped.fq -fastq_maxee 1.0 \
  -fastq_trunclen 350 -fastaout filtered.fa -relabel Filt

# Find unique read sequences and abundances
$usearch -fastx_uniques filtered.fa -sizeout -relabel Uniq -fastaout uniques.fa

# Make 97% OTUs and filter chimeras
$usearch -cluster_otus uniques.fa -otus otus.fa -relabel Otu

# Make OTU table
$usearch -otutab reads_demux.fq -otus otus.fa -otutabout otutab_raw.txt

# Normalize to 1k reads / sample
$usearch -otutab_norm otutab_raw.txt -sample_size 1000 -output otutab.txt

# Alpha diversity
$usearch -alpha_div otutab.txt -output alpha.txt

# Make OTU tree
$usearch -cluster_agg otus.fa -treeout otus.tree

# Beta diversity
mkdir beta/
$usearch -beta_div otutab.txt -tree otus.tree -filename_prefix beta/

# Rarefaction
$usearch -alpha_div_rare otutab.txt -output rare.txt

# Predict taxonomy
$usearch -sintax otus.fa -db ../data/rdp_16s_v16.fa -strand both \
  -tabbedout sintax.txt -sintax_cutoff 0.8

# Taxonomy summary reports
$usearch -sintax_summary sintax.txt -otutabin otutab.txt -rank g -output genus_summary.txt
$usearch -sintax_summary sintax.txt -otutabin otutab.txt -rank p -output phylum_summary.txt
