#!/bin/bash -l

# Step 1: Importing the data into Qiime2

module load bioinfo-tools
module load qiime2
source activate qiime2-2018.11

qiime tools import \
   --type 'SampleData[PairedEndSequencesWithQuality]' \
   --input-path /proj/snic2018-8-341/P10156/P10156_1001/02-FASTQ/180420_M01548_0170_000000000-BRNFY \
   --input-format CasavaOneEightSingleLanePerSampleDirFmt \
   --output-path /proj/snic2018-8-341/18S/demux-paired-end-1.qza


# Step 1b: Summarizing the Imported Data

qiime demux summarize \
   --i-data demux-paired-end-1.qza \
   --o-visualization demux-paired-end-1.qzv


# Step 2: DADA2

qiime dada2 denoise-paired \
   --i-demultiplexed-seqs demux-paired-end-1.qza \
   --p-trim-left-f 9 \
   --p-trim-left-r 7 \
   --p-trunc-len-f 300 \
   --p-trunc-len-r 300 \
   --o-table table-1.qza \
   --o-representative-sequences rep-seqs-1.qza \
   --o-denoising-stats denoising-stats-1.qza


# Step 2b: Summarizing DADA2 

qiime feature-table summarize \
  --i-table table-1.qza \
  --o-visualization table-1.qzv
qiime feature-table tabulate-seqs \
   --i-data rep-seqs-1.qza \
   --o-visualization rep-seqs-1.qzv &
qiime metadata tabulate \
   --m-input-file denoising-stats-1.qza \
   --o-visualization denoising-stats-1.qzv


#Step 3: Rarefaction Curves

qiime feature-table merge \
  --i-tables table-1.qza \
  --i-tables table-2.qza \
  --i-tables table-3.qza \
  --i-tables table-16.qza \
  --i-tables table-17.qza \
  --i-tables table-18.qza \
  --o-merged-table table-Tjarno.qza
qiime feature-table merge-seqs \
  --i-data rep-seqs-1.qza \
  --i-data rep-seqs-2.qza \
  --i-data rep-seqs-3.qza \
  --i-data rep-seqs-16.qza \
  --i-data rep-seqs-17.qza \
  --i-data rep-seqs-18.qza \
  --o-merged-data rep-seqs-Tjarno.qza
qiime feature-table summarize \
  --i-table table-T12.qza \
  --o-visualization table-Tjarno.qzv
qiime feature-table tabulate-seqs \
  --i-data rep-seqs-Tjarno.qza \
  --o-visualization rep-seqs-Tjarno.qzv
qiime diversity alpha-rarefaction \
  --i-table table-Tjarno.qza \
  --p-max-depth 900000 \
  --o-visualization alpha-rarefaction-Tjarno.qzv


#Step 3b: Diversity Measurements

qiime feature-table rarefy \
  --i-table table-ALL-CO1.qza \
  --p-sampling-depth 269894 \
  --o-rarefied-table CO1_ALL_rare.qza
qiime diversity alpha-rarefaction \
  --i-table table-ALL-CO1.qza \
  --p-max-depth 269894 \
  --p-metrics chao1 \
  --m-metadata-file CO1_Metadata.txt \
  --o-visualization alpha_rare_chao1_CO1.qzv
qiime diversity alpha \
  --i-table CO1_ALL_rare.qza \
  --p-metric chao1 \
  --o-alpha-diversity CO1_ALL_chao1.qza
qiime diversity alpha-group-significance \
  --i-alpha-diversity CO1_All_chao1.qza \
  --m-metadata-file Metadata3.txt \
  --o-visualization CO1_GROUPED_chao1_sign.qzv
qiime diversity beta \
  --i-table CO1_ALL_rare.qza \
  --p-metric jaccard \
  --o-distance-matrix CO1_ALL_jaccard.qza 
qiime diversity pcoa \
  --i-distance-matrix CO1_ALL_jaccard.qza \
  --o-pcoa CO1_ALL_jaccard_pcoa.qza
qiime diversity beta-group-significance \
  --i-distance-matrix CO1_ALL_jaccard.qza \
  --m-metadata-file Metadata3.txt \
  --m-metadata-column site \
  --p-method anosim \
  --p-pairwise \q
  --o-visualization CO1_beta_sign.qzv
qiime emperor plot \
  --i-pcoa CO1_ALL_jaccard_pcoa.qza \
  --m-metadata-file Metadata3.txt \
  --o-visualization CO1_ALL_emp_plot.qzv


#Step 4: Initial Identification with SILVA128

qiime feature-classifier classify-sklearn \
  --i-classifier /proj/snic2018-8-341/TestQiime/SILVA128_18S_99_classifier.qza \
  --i-reads rep-seqs-1.qza \
  --o-classification SILVA128_18S_99_taxonomy-1.qza
qiime metadata tabulate \
  --m-input-file SILVA128_18S_99_taxonomy-1.qza \
  --o-visualization SILVA128_18S_99_taxonomy-1.qzv


#Step 5: Creating the reference tree for taxa of interest and identification of OTUs with PPlacer

module load raxml

raxmlHPC -m GTRGAMMA -s Acoela18S.fasta -n ref.tre -f d -p 12354 -# 2 -T 2
raxmlHPC -m GTRGAMMA -f I -t RAxML_bestTree.ref.tre -n root.ref.tre -# 2
raxmlHPC -m GTRGAMMA -f J -p 12354 -t RAxML_rootedTree.root.ref.tre -n conf.root.ref.tre -s Acoela18S.fasta -# 2

source ~/taxtastic/taxtastic-env/bin/activate

taxit create -l 18S_DNA -P Acoela18S.refpkg --aln-fasta Acoela18S.fasta --tree-stats RAxML_info.ref.tre --tree-file RAxML_fastTreeSH_Support.conf.root.ref.tre

module load pplacer

pplacer -o Acoela18S.jplace -p --keep-at-most 20 -c Acoela18S.refpkg Acoela18S_All.fasta
guppy to_csv --point-mass --pp -o Acoela18S.csv Acoela18S.jplace
guppy fat --node-numbers --point-mass --pp -o Acoela18S.phyloxml Acoela18S.jplace