##QIIME2 Script for creating the taxonomy.tsv, feature-table.txt, and tree.nwk files needed to start the analysis on R##
#Create a directory where the outputs from the QIIME2 pipeline will be saved to

cd /data
mkdir parkinsons
cd parkinsons

#Create the demux_seqs.qza file to be able to determine the truncation and trim length for the denoising step
#use the parkinsons_manifest.txt file

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

#Convert the output to a .qzv and transfer it to your personal computer
#in order to visualize it in the qiime2 view website

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

scp <server address>:/data/parkinsons/demux_seqs.qzv <the directory pathway you want your file to be transferred to>

#Denoise the demux_seqs.qza file using DADA2 to attain the table.qza, rep-seqs.qza, and stats.qza files
#using the screen -S <name of the scree> command might be helpful as this takes a long time to complete.
#Based on our result, we decided to only trim the first 8 bases out of the total 251 bases

qiime dada2 denoise-single \
  --i-demultiplexed-seqs demux_seqs.qza \
  --p-trim-left 8 \
  --p-trunc-len 251 \
  --o-representative-sequences rep-seqs.qza \
  --o-table table.qza \
  --o-denoising-stats stats.qza

#Run a taxonomic analysis by using the silva-138-99-515-806-nb-classifier.qza classifier

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

#Then create the phylogenetic tree by using the rep-seqs.qza file

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

#Create a parkinsons_export folder, then use the export function to export the taxonomy.qza,
#rooted-tree.qza, and table.qza files into this folder

cd /data
mkdir parkinsons_export
cd parkinsons

qiime tools export \
   --input-path taxonomy.qza \
   --output-path /data/parkinsons_export/taxonomy_export

qiime tools export \
   --input-path rooted-tree.qza \
   --output-path /data/parkinsons_export/rooted_tree_export

qiime tools export \
   --input-path table.qza \
   --output-path /data/parkinsons_export/table_export

#Before transferring this folder to your personal computer to use the files in R studios
#convert the exported table.qza which is a .biom file into a .txt

cd /data/parkinsons_export/table_export

biom convert \
-i feature-table.biom \
--to-tsv \
-o feature-table.txt

#Transfer the whole folder onto your computer

scp -r <server address>:/data/parkinsons_export <the directory pathway you want your folder to be transferred to>

#Then transfer the parkinsons_metadata.txt file onto your personal computer

scp <server address>:/mnt/datasets/project_2/parkinsons/parkinsons_metadata.txt <the directory pathway you want your file to be transferred to>



#Rarefaction process script done in R studios using the transferred taxonomy.tsv, tree.nwk, feature-table.txt, and 
#parkinsons_metadata.txt van be found in the R script called Rarefaction_parkinsons.R
#Rest of the analyses done in R can be found in the Diversity_Metrics+Statistical_Testing+Binning.R and Microbiome_Analysis.R files



###Our group did not create an alpha-rarefaction curve for determining the###   
###sampling depth for rarefaction because we did this process in R, but here is the code for both parts of it###

#To set your max depth, visualize the table.qza file that is converted to .qzv, and base the depth 
#by using the maximum frequency value
#X is your determined max depth value, which for our data would have been around 23000

qiime diversity alpha-rarefaction \
  --i-table table.qza \
  --i-phylogeny rooted-tree.qza \
  --p-max-depth X \
  --m-metadata-file /mnt/datasets/project_2/parkinsons/parkinsons_metadata.txt \
  --o-visualization alpha-rarefaction.qzv

#Visualize the alpha-rarefaction.qzv file in order to set you sampling depth for rarefaction
#process that creates both alpha- and beta-diversity metrics
#Y is your determined sampling depth which for our data would have been around 5000 
#(based on the sequecning depth where most of the reads plateau)

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table.qza \
  --p-sampling-depth Y \
  --m-metadata-file /mnt/datasets/project_2/parkinsons/parkinsons_metadata.txt \
  --output-dir core-metrics-results
