#make a directory and move into it

mkdir ./HiSeas
cd ./ HiSeas

#tabulate the metadata to ensure it is sorted correctly by Qiime2

qiime metadata tabulate \
  --m-input-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization hiseasmetadata.qzv

#import demultiplexed sequences

qiime tools import \
  --type "SampleData[SequencesWithQuality]" \
  --input-format SingleEndFastqManifestPhred33V2 \
  --input-path /mnt/datasets/project_2/hiseas/hiseas_manifest.txt \
  --output-path ./demux_seqs.qza

#summarize demultiplexed sequences

qiime demux summarize \
  --i-data ./demux_seqs.qza \
  --o-visualization ./demux_seqs.qzv

#denoise demultiplexed sequences
#set truncation length to 292, because at this value the bottom 9th and 25th percentiles are still at a quality of 24 and 34 #respectively, while at the 50th percentile the quality is high, at 37
#There is an incredibly steep drop-off in quality immediately after that, starting at 293 the quality at the bottom 9th and #25th percentiles are 10 and 23 respectively. The quality at the 50th percentile is 34, so low quality reads have started to #influence the quality of the majority - at position 294, for comparison the 50th percentile has dropped to 18, unnacceptable

qiime dada2 denoise-single \
  --i-demultiplexed-seqs ./demux_seqs.qza \
  --p-trunc-len 292 \
  --o-table ./dada2_table.qza \
  --o-representative-sequences ./dada2_rep_set.qza \
  --o-denoising-stats ./dada2_stats.qza

#tabulate denoising stats

qiime metadata tabulate \
  --m-input-file ./dada2_stats.qza  \
  --o-visualization ./dada2_stats.qzv

#output feature table summary

qiime feature-table summarize \
  --i-table ./dada2_table.qza \
  --m-sample-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./dada2_table.qzv

#import the silva sequence database (version 138, 515F/806R)

wget \
  -O "sepp-refs-silva-138.qza" \
  "https://data.qiime2.org/2021.8/common/sepp-refs-silva-128.qza"

#use the database to generate a phylogenetic tree for diversity analysis, set threads to 8 for this server

qiime fragment-insertion sepp \
  --i-representative-sequences ./dada2_rep_set.qza \
  --i-reference-database sepp-refs-silva-138.qza \
  --o-tree ./tree.qza \
  --o-placements ./tree_placements.qza \
  --p-threads 8 

#Use a pre-trained classifier already on the server *(trained on the same silva version used for phylogeny) to assign #taxonomy

qiime feature-classifier classify-sklearn \
  --i-classifier /mnt/datasets/classifiers/silva-138-99-515-806-nb-classifier.qza \
  --i-reads ./dada2_rep_set.qza \
  --o-classification taxonomy.qza

#visualize taxonomy

qiime metadata tabulate \
  --m-input-file ./taxonomy.qza \
  --o-visualization ./taxonomy.qzv

#generate interactive BLAST table for all taxonomic assignments

qiime feature-table tabulate-seqs \
  --i-data ./dada2_rep_set.qza \
  --o-visualization ./dada2_rep_set.qzv

#filter out mitochondria and chloroplasts from the feature table

qiime taxa filter-table \
  --i-table ./dada2_table.qza \
  --i-taxonomy taxonomy.qza \
  --p-exclude mitochondria,chloroplast \
  --o-filtered-table table-no-mitochondria-no-chloroplast.qza

#filter out all data except plastic and wood from the feature table

qiime feature-table filter-samples \
  --i-table table-no-mitochondria-no-chloroplast.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --p-where "[orig_env_material] IN ('wood', 'plastic')" \
  --o-filtered-table wood-plastic-filtered-table.qza

#create alpha rarefaction curves 

qiime diversity alpha-rarefaction \
  --i-table ./wood-plastic-filtered-table.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./alpha_rarefaction_curves.qzv \
  --p-min-depth 10 \
  --p-max-depth 47391

#generate taxonomy barchart for rarefaction depth 20000 and 36000, based on the 
#alpha rarefaction curves

qiime feature-table filter-samples \
  --i-table ./wood-plastic-filtered-table.qza \
  --p-min-frequency 20000 \
  --o-filtered-table ./table_20k.qza

qiime taxa barplot \
  --i-table ./table_20k.qza \
  --i-taxonomy ./taxonomy.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./taxa_barplot20k.qzv

qiime feature-table filter-samples \
  --i-table ./wood-plastic-filtered-table.qza \
  --p-min-frequency 36000 \
  --o-filtered-table ./table_36k.qza

qiime taxa barplot \
  --i-table ./table_36k.qza \
  --i-taxonomy ./taxonomy.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./taxa_barplot36k.qzv

#calculate core diversity metrics for rarefaction depth 20000 and 36000, based on the 
#alpha rarefaction curves

qiime diversity core-metrics-phylogenetic \
  --i-table ./wood-plastic-filtered-table.qza \
  --i-phylogeny ./tree.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --p-sampling-depth 20000 \
  --output-dir ./core-metrics-results-20000

qiime diversity core-metrics-phylogenetic \
  --i-table ./wood-plastic-filtered-table.qza \
  --i-phylogeny ./tree.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --p-sampling-depth 36000 \
  --output-dir ./core-metrics-results-36000

#Use core diversity metrics to determine significance of differences between sample surface material #and usage

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column orig_env_material \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-material-significance.qzv

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column description \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-location-significance.qzv

#Use core diversity metrics for pairwise comparisons of weighted unifrac distances
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column orig_env_material \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-material-pairwise-significance.qzv \
  --p-pairwise

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column description \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-location-pairwise-significance.qzv \
  --p-pairwise

#repeat using permdisp method to determine if differences are due to within group variation
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column orig_env_material \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-material-permdisp-significance.qzv \
  --p-method permdisp

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-20000-filtered/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --m-metadata-column description \
  --o-visualization core-metrics-results-20000-filtered/weighted-unifrac-location-permdisp-significance.qzv \
  --p-method permdisp

#Alpha diversity
#Observed features
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-20000-filtered/observed_features_vector.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./core-metrics-results-20000-filtered/observed_features_statistics.qzv

#Shannon diversity
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-20000-filtered/shannon_vector.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./core-metrics-results-20000-filtered/shannon_statistics.qzv

#Faiths 
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-20000-filtered/faith_pd_vector.qza \
  --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
  --o-visualization ./core-metrics-results-20000-filtered/faiths_pd_statistics.qzv

#Piellou's
qiime diversity alpha-group-significance \
 --i-alpha-diversity core-metrics-results-20000-filtered/evenness_vector.qza \
 --m-metadata-file /mnt/datasets/project_2/hiseas/hiseas_metadata.txt \
 --o-visualization ./core-metrics-results-20000-filtered/evenness_statistics.qzv






