#!/usr/bin/env

#### MDMA002 CODE ####

# 1) We will import (combine the fastq and barcode sequencing files) and demlutiplex (assigning which reads belong to which sample) using the Colombia manifest file.
# Need to run the code in the woring directory of the manifest file

qiime tools import \
--type "SampleData[SequencesWithQuality]" \
--input-format SingleEndFastqManifestPhred33V2 \
--input-path /mnt/datasets/project_2/colombia/colombia_manifest.txt \
--output-path /data/colombian/demux.qza


# 2) Create a visualization of demultiplexed samples.
qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv

# 3) Denoising (Finding the appropriate trimming parameters and correcting sequencing errors) and clustering (finding the number of unique sequences we have based on ASVs) using the demux.qza file.    
# Occurs all in one-step based on the following code: 

qiime dada2 denoise-single --i-demultiplexed-seqs demux.qza \
--p-trim-left 0 \ 
--p-trunc-len 225 \
--o-representative-sequences rep-seqs.qza \
--o-table table.qza \
--o-denoising-stats stats.qza

# 4) Visualization of the representative sequences, table and stats files:

# a) DADA2 Stats: 
qiime metadata tabulate \
--m-input-file stats.qza \
--o-visualization stats.qzv

# b) ASV stats: 
qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt

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

# 5) Aligning our ASV reads with known reference sequences of microbes using a database.

# a) First, we will need to train a classifier. The original study had used the V4 hyper variable region of the 16S rRNA gene with the following primers:    

qiime feature-classifier extract-reads \
--i-sequences /mnt/datasets/silva_ref_files/silva-138-99-seqs.qza \
--p-f-primer GTGCCAGCMGCCGCGGTAA \
--p-r-primer GGACTACHVGGGTWTCTAAT \
--p-trunc-len 225 \
--o-reads ref-seqs-trimmed.qza

qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads ref-seqs-trimmed.qza \
--i-reference-taxonomy /mnt/datasets/silva_ref_files/silva-138-99-tax.qza \
--o-classifier classifier.qza

# b) Now, we will use the trained classifier to assign taxonomy to our reads: 

qiime feature-classifier classify-sklearn \
--i-classifier classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza

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

# 6) Creating the alpha rarefaction curve 
# a) First, we will generate a tree for phylogenetic 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

# b) Next, we will create our alpha-rarefaction curve based on the table.qza and rooted-tree.qza files that we generated.

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

# 7) We will be filtering our metadata based on taxonomy and our research objective. 

# a) Taxonomy-based filtering: We will be removing mitochondrial and chloroplast sequences as an quality control step:

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

qiime feature-table summarize \
--i-table table-no-mitochondria-no-chloroplast.qza \
--o-visualization table-no-mitochondria-no-chloroplast.qzv \
--m-sample-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt 

# b) Metadata-based filtering: We want to only keep the columns which we will be analyzing in downstream.  
# Those columns are: "#SampleID", "Cardiometabolic_status", age_range", "city", "Total_Cholesterol", "sex", "smoker" 

# We will be generating 4 tables of cardiometabolically healthy and abnormal samples as well as smokers and nonsmokers which will be used for downstream analysis

# Cardio Healthy: 
qiime feature-table filter-samples \
--i-table table-no-mitochondria-no-chloroplast.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--p-where "[Cardiometabolic_status]='Healthy'" \ 
--o-filtered-table cardio_healthy_filtered_table.qza

# Cardio Abnormal: 
qiime feature-table filter-samples \
--i-table table-no-mitochondria-no-chloroplast.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--p-where "[Cardiometabolic_status]='Abnormal'" \ 
--o-filtered-table CA_filtered_table.qza

# Smoker: 
qiime feature-table filter-samples \
--i-table table-no-mitochondria-no-chloroplast.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--p-where "[Cardiometabolic_status]='Healthy'" \ 
--o-filtered-table smoker_filtered_table.qza

# Nonsmoker: 
qiime feature-table filter-samples \
--i-table table-no-mitochondria-no-chloroplast.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--p-where "[Cardiometabolic_status]='Abnormal'" \ 
--o-filtered-table nonsmoker_filtered_table.qza


#### MDMA003 Code ####

# 1) Perfomring alpha and beta diversity metrics on the no mitochondria and chloroplast table:

qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table-no-mitochondria-no-chloroplast.qza \
--p-sampling-depth 20291 \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--output-dir main_core_metrics_results

# We will convert the files into .qzv to be able to visualize them. 

# Alpha Diversity Metrics 

qiime diversity alpha-group-significance \
--i-alpha-diversity core_metrics_results/faith_pd_vector.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--o-visualization core-metrics-results/faith-pd-group-significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity core_metrics_results/evenness_vector.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--o-visualization core-metrics-results/evenness-group-significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity core_metrics_results/observed_features_vector.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--o-visualization core-metrics-results/observed_features-group-significance.qzv

qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/shannon_vector.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--o-visualization core-metrics-results/shannon-group-significance.qzv

# Beta Diversity Metrics based on Cardiometabolic Status 

qiime diversity beta-group-significance \
--i-distance-matrix core_metrics_results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--m-metadata-column Cardiometabolic_status \ 
--o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
--p-pairwise

qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--m-metadata-column Cardiometabolic_status \ 
--o-visualization core-metrics-results/weighted_unifrac_distance_matrix.qzv \
--p-pairwise

qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/bray_curtis_distance_matrix.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--m-metadata-column Cardiometabolic_status \ 
--o-visualization core-metrics-results/bray_curtis_distance_matrix.qzv \
--p-pairwise

qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/jaccard_distance_matrix.qza \
--m-metadata-file /mnt/datasets/project_2/colombia/colombia_metadata.txt \
--m-metadata-column Cardiometabolic_status \ 
--o-visualization core-metrics-results/jaccard_distance_matrix.qzv \
--p-pairwise


#### MDMA004 Export Code ####

#To create Phyloseq files, we need the following files from QIIME:
#Metada, featuretable (table), taxonomy, rootedtree.

#CODE USED ON QIIME
#export table.qza
qiime tools export \
--input-path table.qza \
--output-path table_EXPORT

#export table-no-mitochondria-no-chloroplast.qza
qiime tools export \
--input-path table-no-mitochondria-no-chloroplast.qza \
--output-path table-noMC_EXPORT

#export taxonomy.qza
qiime tools export \
--input-path taxonomy.qza \
--output-path taxonomy_EXPORT

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







