#Team3_script#2 447_21WT1.txt

###The following script outlines the code ran for the analyses performed in QIIME2###

#Create symbolic link while in the working directory
ln -S /mnt/datasets/project_2/infant

#Import filtered metadata and manifest file into server using the scp command

# Import demultiplexed data into data directory with the rearranged manifest file
qiime tools import \
  --type "SampleData[SequencesWithQuality]" \
  --input-format SingleEndFastqManifestPhred33V2 \
  --input-path ./rearranged_manifest.tsv \
  --output-path ./demux_seqs.qza

#Visualize demultiplexed samples 
qiime demux summarize \
  --i-data demux_seqs.qza \
  --o-visualization demux_seqs.qzv

#Perform sequence quality control and identify ASVs using DADA2
qiime dada2 denoise-single \
  --i-demultiplexed-seqs ./demux_seqs.qza \
  --p-trunc-len 150 \
  --o-table ./dada2_table.qza \
  --o-representative-sequences ./dada2_rep_set.qza \
  --o-denoising-stats ./dada2_stats.qza

#Visualize DADA2 stats
qiime metadata tabulate \
  --m-input-file ./dada2_stats.qza  \
  --o-visualization ./dada2_stats.qzv

#Generate and visualize ASV feature table
qiime feature-table summarize \
  --i-table ./dada2_table.qza \
  --m-sample-metadata-file ./combined_filtered_infant_metadata.tsv \
  --o-visualization ./dada2_table.qzv

#Rename files
mv dada2_rep_set.qza rep-seqs.qza
mv dada2_table.qza table.qza

#Generate phylogenetic tree for diversity analyses
qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza

#Generate alpha rarefaction curve to determine an appropriate rarefaction depth
qiime diversity alpha-rarefaction \
  --i-table table.qza \
  --i-phylogeny rooted-tree.qza \
  --p-max-depth 48000 \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --o-visualization alpha-rarefaction.qzv

#Assign taxonomy to sequences
qiime feature-classifier classify-sklearn \
  --i-classifier /mnt/datasets/classifiers/silva-138-99-515-806-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza

#Generate Taxonomy barplots
qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --o-visualization taxa-bar-plots.qzv

#Filter feature table to remove Chloroplast and Mitochondrial sequences
qiime taxa filter-table \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --p-exclude mitochondria,chloroplast \
  --o-filtered-table table-no-mitochondria-no-chloroplast.qza

#Calculate alpha and beta diversity metrics with a sampling depth of 15,000
qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table-no-mitochondria-no-chloroplast.qza \
  --p-sampling-depth 15000 \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --output-dir core-metrics-results-15000

#Calculate alpha group significance and generate boxplots for the following alpha diversity metrics:
#Observed features
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-15000/observed_features_vector.qza \
  --m-metadata-file  /data/combined_filtered_infant_metadata.tsv \
  --o-visualization core-metrics-results-15000/observed-features.qzv\
  --p-pairwise

#Shannon Diversity
qiime diversity alpha-group-significance \
  --i-alpha-diversity /data/core-metrics-results-15000/shannon_vector.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --o-visualization /data/core-metrics-results-15000/shannon-significance.qzv
  --p-pairwise

#Faith's Phylogenetic Distance
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-15000/faith_pd_vector.qza \
  --m-metadata-file  /data/combined_filtered_infant_metadata.tsv \
  --o-visualization core-metrics-results-15000/faith-pd-group-significance.qzv\
  --p-pairwise

#Pielou's Evenness
qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-15000/evenness_vector.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --o-visualization core-metrics-results-15000/evenness-group-significance.qzv
  --p-pairwise

#Calculate beta group significance and generate boxplots for the following beta diversity metrics:
#unweighted unifrac distance
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-15000/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --m-metadata-column feed\
  --o-visualization core-metrics-results-15000/unweighted-unifrac-feed-significance.qzv\
  --p-pairwise

#weighted unifrac distance
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-15000/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --m-metadata-column feed \
  --o-visualization core-metrics-results-15000/weighted-unifrac-feed-significance.qzv\
  --p-pairwise

#Jaccard
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-15000/jaccard_distance_matrix.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --m-metadata-column feed \
  --o-visualization core-metrics-results-15000/jaccard-distance-feed-significance.qzv\
  --p-pairwise

#Bray-Curtis
qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-15000/bray_curtis_distance_matrix.qza \
  --m-metadata-file /data/combined_filtered_infant_metadata.tsv \
  --m-metadata-column feed \
  --o-visualization core-metrics-results-15000/bray-curtis-feed-significance.qzv\
  --p-pairwise

# Export BIOM data from QIIME
qiime tools export \
 --input-path table.qza \
 --output-path exported

qiime tools export \
 --input-path taxonomy.qza \
 --output-path exported

qiime tools export \
 --input-path rooted-tree.qza \
 --output-path exported

# Edit the column names exported/taxonomy.tsv
sed \
  '1 s/Feature ID/#OTUID/;1 s/Taxon/taxonomy/;1 s/Confidence/confidence/' \
  exported/taxonomy.tsv \
  > exported/biom-taxonomy.tsv

#Combine Taxonomy with BIOM data 
biom add-metadata \
  -i exported/feature-table.biom \
  -o exported/table-with-taxonomy.biom \
  --observation-metadata-fp exported/biom-taxonomy.tsv \
  --sc-separated taxonomy