#!/usr/bin/env ruby

require 'optimist'

ARGV.push("--help") if ARGV.empty?
opts = Optimist::options do
  banner File.basename($0)
  opt :input_dir, "dir with fastqs", :required=>true, :type=>:string
  opt :classification_db, "qiime classification db", :required=>true,
      :type=>:string
  opt :metadata, "metadata tsv", :required=>true, :type=>:string
  opt :threads, "number of threads to use", :default=>4
  opt :f_trim, "number of bases to trim on the forward read", :default=>19
  opt :r_trim, "number of bases to trim on the reverse read", :default=>20
  opt :f_trunc, "number of bases to truncate forward read", :default=>250
  opt :r_trunc, "number of bases to truncate forward read", :default=>250
end

if !File.exists? "paired.qza"
  STDERR << "Pairing sequences...\n"
  system("qiime tools import --type " +
         "'SampleData[PairedEndSequencesWithQuality]' --input-path " +
         opts.input_dir + " --output-path paired.qza " +
         "--source-format CasavaOneEightSingleLanePerSampleDirFmt")
end
	
if !File.exists? "paired.qzv"
  system "qiime demux summarize --i-data paired.qza " +
    "--o-visualization paired.qzv"
end

if !File.exists? "table.qza"
  STDERR << "Running DADA2...\n"
  system "qiime dada2 denoise-paired \
  --verbose \
  --p-n-threads #{opts.threads} \
  --i-demultiplexed-seqs paired.qza \
  --o-table table.qza \
  --o-representative-sequences rep-seqs.qza \
  --p-trim-left-f #{opts.f_trim} \
  --p-trim-left-r #{opts.r_trim} \
  --p-trunc-len-f #{opts.f_trunc} \
  --p-trunc-len-r #{opts.r_trunc}"
  system "qiime tools export rep-seqs.qza --output-dir ."
  system "qiime tools export table.qza --output-dir ."
end

if !File.exists? "aligned-rep-seqs.qza"
  STDERR << "Aligning Sequences...\n"
  system "qiime alignment mafft \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza"
end

if !File.exists? "masked-aligned-rep-seqs.qza"
  system "qiime alignment mask \
  --i-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza"
end

if !File.exists? "unrooted-tree.qza"
  STDERR << "Making Phylogeny...\n"
  system "qiime phylogeny fasttree \
  --i-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza"
end

if !File.exists? "rooted-tree.qza"
  system "qiime phylogeny midpoint-root \
  --i-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza"
  system "qiime tools export rooted-tree.qza --output-dir ."
end

if !File.exists? "taxonomy.qza"
  STDERR << "Classifying Sequences...\n"
  system "qiime feature-classifier classify-sklearn --i-reads rep-seqs.qza \
  --i-classifier #{opts.classification_db} --o-classification taxonomy.qza"
  system "qiime tools export taxonomy.qza --output-dir ."
end

if !File.exists? "out.biom"
  system "biom add-metadata -i feature-table.biom -o out.biom \
  --output-as-json --observation-metadata-fp taxonomy.tsv \
  --sample-metadata-fp #{opts.metadata} \
  --observation-header 'Id,Taxonomy,Confidence'"
end