#!/bin/bash

#PIPELINE FROM IMPORTED SEQUENCES TO ABUNDANCE TABLES USING QIIME2-2022.2
#3 SEPARATE RUNS: A, B AND C
#RUN-A: FARM NASAL SAMPLES
#RUN-B: FARM RECTAL SAMPLES
#RUN-C: FIVE PIGLETS DIFFERENT NASAL SAMPLES

##################################################
##################################################
#1: IMPORT
#EDIT THE THREE MANIFEST FILES TO ADD THE PATH WERE SEQUENCES ARE

runs='A B C'
for i in $runs; do

qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path manifest-${i} \
  --output-path paired-end-demux-${i}.qza \
  --input-format PairedEndFastqManifestPhred33 ;

qiime cutadapt trim-paired \
  --i-demultiplexed-sequences paired-end-demux-${i}.qza \
  --p-cores 24 \
  --p-front-f CCTACGGGNGGCWGCAG \
  --p-front-r NACTACHVGGGTATCTAATCC \
  --p-match-adapter-wildcards \
  --p-match-read-wildcards \
  --p-discard-untrimmed \
  --o-trimmed-sequences trimmed-demux-${i}.qza ;

qiime demux summarize \
  --i-data trimmed-demux-${i}.qza \
  --p-n 100000 \
  --o-visualization trimmed-demux-${i}.qzv ;

done ;

##################################################
##################################################
#2: DENOISING WITH DADA2

mkdir dada2-files ;

#RunA
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs trimmed-demux-A.qza \
  --p-trunc-len-f 267 \
  --p-trunc-len-r 182 \
  --o-table dada2-files/table-A.qza \
  --o-representative-sequences dada2-files/rep-seqs-A.qza \
  --o-denoising-stats dada2-files/denoising-stats-A.qza \
  --p-n-threads 12 ;

# Denoising stats visualization
qiime metadata tabulate \
  --m-input-file dada2-files/denoising-stats-A.qza \
  --o-visualization dada2-files/denoising-stats-A.qzv ;

#RunB

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs trimmed-demux-B.qza \
  --p-trunc-len-f 247 \
  --p-trunc-len-r 198 \
  --o-table dada2-files/table-B.qza \
  --o-representative-sequences dada2-files/rep-seqs-B.qza \
  --o-denoising-stats dada2-files/denoising-stats-B.qza \
  --p-n-threads 12 ;

# Denoising stats visualization
qiime metadata tabulate \
  --m-input-file dada2-files/denoising-stats-B.qza \
  --o-visualization dada2-files/denoising-stats-B.qzv ;

#RunC

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs trimmed-demux-C.qza \
  --p-trunc-len-f 271 \
  --p-trunc-len-r 182 \
  --o-table dada2-files/table-C.qza \
  --o-representative-sequences dada2-files/rep-seqs-C.qza \
  --o-denoising-stats dada2-files/denoising-stats-C.qza \
  --p-n-threads 12 ;

# Denoising stats visualization
qiime metadata tabulate \
  --m-input-file dada2-files/denoising-stats-C.qza \
  --o-visualization dada2-files/denoising-stats-C.qzv ;

##################################################
##################################################
#3: TABLE AND SEQUENCES MERGING


qiime feature-table merge \
  --i-tables dada2-files/table-A.qza \
  --i-tables dada2-files/table-B.qza \
  --i-tables dada2-files/table-C.qza \
  --o-merged-table dada2-files/table-merged.qza ;

qiime feature-table merge-seqs \
  --i-data dada2-files/rep-seqs-A.qza \
  --i-data dada2-files/rep-seqs-B.qza \
  --i-data dada2-files/rep-seqs-C.qza \
  --o-merged-data dada2-files/rep-seqs-merged.qza ;
  
qiime feature-table tabulate-seqs \
  --i-data dada2-files/rep-seqs-merged.qza \
  --o-visualization dada2-files/rep-seqs-merged.qzv ; 

##################################################
##################################################
#4: NON PROKARYOTIC SEQUENCES REMOVAL

qiime quality-control exclude-seqs \
  --i-query-sequences dada2-files/rep-seqs-merged.qza \
  --i-reference-sequences 88_otus.qza \
  --p-method vsearch \
  --p-perc-identity 0.65 \
  --p-perc-query-aligned 0.5 \
  --p-threads 12 \
  --o-sequence-hits dada2-files/filtered-seqs.qza \
  --o-sequence-misses dada2-files/misses.qza ;

qiime feature-table filter-features \
  --i-table dada2-files/table-merged.qza \
  --m-metadata-file dada2-files/filtered-seqs.qza \
  --o-filtered-table dada2-files/filtered-table.qza ;

##################################################
##################################################
#5: ASVs TAXONOMY

qiime feature-classifier classify-sklearn \
  --i-classifier gg_v3v4_classifier.qza \
  --i-reads dada2-files/filtered-seqs.qza \
  --o-classification taxonomy.qza ;

##################################################
##################################################
#6: OTHER 16S REMOVAL

qiime taxa filter-table \
  --i-table dada2-files/filtered-table.qza \
  --i-taxonomy taxonomy.qza \
  --p-exclude Archaea,Chloroplast,mitochondria \
  --o-filtered-table final-table.qza ;

qiime feature-table filter-seqs \
  --i-data dada2-files/filtered-seqs.qza \
  --i-table final-table.qza \
  --p-no-exclude-ids \
  --o-filtered-data dada2-files/final-rep-seqs.qza ;

##################################################
##################################################
#7: COMPUTATION OF PHYLOGENETIC TREE


# 2.1. Multiple sequence alignment with MAFFT
qiime alignment mafft \
  --i-sequences dada2-files/final-rep-seqs.qza  \
  --o-alignment dada2-files/final-aligned-rep-seqs.qza \
  --p-n-threads 8 ;

# 2.2. Masking the hypevariable positions in the alignement
qiime alignment mask \
  --i-alignment dada2-files/final-aligned-rep-seqs.qza \
  --o-masked-alignment dada2-files/final-masked-aligned-rep-seqs.qza ;

# 2.3. Build phylogenetic tree with FastTree
qiime phylogeny fasttree \
  --i-alignment dada2-files/final-masked-aligned-rep-seqs.qza \
  --o-tree dada2-files/unrooted-tree.qza \
  --p-n-threads 8 ;

# 2.4. Tree rooting
qiime phylogeny midpoint-root \
  --i-tree dada2-files/unrooted-tree.qza \
  --o-rooted-tree dada2-files/rooted-tree.qza ;

##################################################
##################################################
#8: NASAL DEEP SAMPLES DATASET

mkdir deep ;

#Table filtering (dna and rna)
qiime feature-table filter-samples \
  --i-table final-table.qza \
  --m-metadata-file metadata.tsv \
  --p-where "[run]='C' AND [dna]='dna'" \
  --o-filtered-table deep/deep-table-dna.qza

qiime feature-table filter-samples \
  --i-table final-table.qza \
  --m-metadata-file metadata.tsv \
  --p-where "[run]='C' AND [dna]='rna'" \
  --o-filtered-table deep/deep-table-rna.qza
  
#Negative control ASVs filtering

qiime feature-table filter-features \
  --i-table deep/deep-table-dna.qza \
  --m-metadata-file asv-to-exclude-dna.csv \
  --p-exclude-ids \
  --o-filtered-table deep/deep-table-dna-filtered.qza

qiime feature-table filter-features \
  --i-table deep/deep-table-rna.qza \
  --m-metadata-file asv-to-exclude-rna.csv \
  --p-exclude-ids \
  --o-filtered-table deep/deep-table-rna-filtered.qza
  
#DNA nasal samples diversity analysis

qiime feature-table filter-samples \
  --i-table deep/deep-table-dna-filtered.qza \
  --m-metadata-file metadata.tsv \
  --p-where "[bodysite]='nasal'" \
  --o-filtered-table deep/table-nasal-dna-filtered.qza ;

depth=2291

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny dada2-files/rooted-tree.qza \
  --i-table deep/table-nasal-dna-filtered.qza \
  --p-sampling-depth $depth \
  --m-metadata-file metadata.tsv \
  --output-dir deep/core-metrics \
  --p-n-jobs-or-threads 4 ;

qiime diversity alpha \
  --i-table deep/table-nasal-dna-filtered.qza \
  --p-metric chao1 \
  --o-alpha-diversity deep/core-metrics/chao1_vector.qza 

qiime diversity alpha-group-significance \
  --i-alpha-diversity deep/core-metrics/chao1_vector.qza \
  --m-metadata-file metadata.tsv \
  --o-visualization deep/core-metrics/chao1-group-significance.qzv

qiime diversity beta-group-significance \
  --i-distance-matrix deep/core-metrics/jaccard_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column depth \
  --o-visualization deep/core-metrics/jaccard_depth-significance.qzv \
  --p-pairwise ;

qiime diversity beta-group-significance \
  --i-distance-matrix deep/core-metrics/bray_curtis_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --m-metadata-column depth \
  --o-visualization deep/core-metrics/bray_curtis_depth-significance.qzv \
  --p-pairwise

qiime diversity adonis \
  --i-distance-matrix deep/core-metrics/jaccard_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --p-formula depth \
  --o-visualization deep/core-metrics/jaccard_ADONIS-depth.qzv ;

qiime diversity adonis \
  --i-distance-matrix deep/core-metrics/bray_curtis_distance_matrix.qza \
  --m-metadata-file metadata.tsv \
  --p-formula depth \
  --o-visualization deep/core-metrics/bray_curtis_ADONIS-depth.qzv ;

#Relative frequency table generation

mkdir deep/rel-freq ;

sequencing='dna rna'
levels="2 3 4 5 6 7"

for a in $sequencing; do for i in $levels; do

qiime taxa collapse \
  --i-table deep/deep-table-${a}-filtered.qza \
  --i-taxonomy taxonomy.qza \
  --p-level ${i} \
  --o-collapsed-table deep/rel-freq/table-${a}-l${i}.qza

qiime feature-table relative-frequency \
  --i-table deep/rel-freq/table-${a}-l${i}.qza \
  --o-relative-frequency-table deep/rel-freq/table-rel-freq-${a}-l${i}.qza 
 
qiime tools export \
  --input-path deep/rel-freq/table-rel-freq-${a}-l${i}.qza \
  --output-path deep/rel-freq/l${i}-rel-freq-table-${a} 

biom convert \
-i deep/rel-freq/l${i}-rel-freq-table-${a}/feature-table.biom \
-o deep/rel-freq/l${i}-rel-freq-table-${a}.tsv --to-tsv

#

qiime feature-table relative-frequency \
  --i-table deep/deep-table-${a}-filtered.qza \
  --o-relative-frequency-table deep/rel-freq/table-rel-freq-${a}-asv.qza 
 
qiime tools export \
  --input-path deep/rel-freq/table-rel-freq-${a}-asv.qza \
  --output-path deep/rel-freq/rel-freq-asv-table-${a} 

biom convert \
-i deep/rel-freq/rel-freq-asv-table-${a}/feature-table.biom \
-o deep/rel-freq/rel-freq-asv-table-${a}.tsv --to-tsv

done ;
done ;

##################################################
##################################################
#: GENERATION OF FARM-SAMPLES RELATIVE FREQUENCY TABLE


mkdir rel-freq-farms ;

qiime feature-table filter-samples \
  --i-table final-table.qza \
  --m-metadata-file metadata.tsv \
  --p-where "[study]='farms'" \
  --o-filtered-table table-farms.qza

levels="2 3 4 5 6 7"

for i in $levels; do

qiime taxa collapse \
  --i-table table-farms.qza \
  --i-taxonomy taxonomy.qza \
  --p-level ${i} \
  --o-collapsed-table rel-freq-farms/table-l${i}.qza

qiime feature-table relative-frequency \
  --i-table rel-freq-farms/table-l${i}.qza \
  --o-relative-frequency-table rel-freq-farms/table-rel-freq-l${i}.qza 
 
qiime tools export \
  --input-path rel-freq-farms/table-rel-freq-l${i}.qza \
  --output-path rel-freq-farms/l${i}-rel-freq-table 

biom convert \
-i rel-freq-farms/l${i}-rel-freq-table/feature-table.biom \
-o rel-freq-farms/l${i}-rel-freq-table-farms.tsv --to-tsv

done ;

#asv table

qiime feature-table relative-frequency \
  --i-table table-farms.qza \
  --o-relative-frequency-table rel-freq-farms/table-rel-freq-asv.qza 
 
qiime tools export \
  --input-path rel-freq-farms/table-rel-freq-asv.qza \
  --output-path rel-freq-farms/rel-freq-asv-table 

biom convert \
-i rel-freq-farms/rel-freq-asv-table/feature-table.biom \
-o rel-freq-farms/rel-freq-asv-table-farms.tsv --to-tsv


##################################################
##################################################

