#
# Codes for the gut microbiome analysis
#
# [Analysis on Qiime2]
# 
# Sequencing was carried out on two separate sequences, so each was denoise and then merged.
#

# Data import1

qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path ManifestFile1.csv --output-path demux1.qza --input-format PairedEndFastqManifestPhred33

qiime demux summarize --i-data demux1.qza --o-visualization demux1.qzv


# Data denoising1

qiime dada2 denoise-paired \
--i-demultiplexed-seqs demux1.qza \
--p-trunc-len-f 270 --p-trunc-len-r 210 \
--p-trim-left-f 0 --p-trim-left-r 0 \
--p-trunc-q 2 \
--p-n-reads-learn 1000000 \
--p-max-ee-f Inf --p-max-ee-r Inf \
--p-n-threads 4 --o-table table1.qza \
--o-denoising-stats stats1.qza \
--o-representative-sequences rep-seqs1.qza \
--verbose


# Data import2

qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path ManifestFile2.csv --output-path demux2.qza --input-format PairedEndFastqManifestPhred33

qiime demux summarize --i-data demux2.qza --o-visualization demux2.qzv


# Data denoising2

qiime dada2 denoise-paired \
--i-demultiplexed-seqs demux2.qza \
--p-trunc-len-f 260 --p-trunc-len-r 260 \
--p-trim-left-f 0 --p-trim-left-r 0 \
--p-trunc-q 2 \
--p-n-reads-learn 1000000 \
--p-max-ee-f Inf --p-max-ee-r Inf \
--p-n-threads 4 --o-table table.2qza \
--o-denoising-stats stats2.qza \
--o-representative-sequences rep-seqs2.qza \
--verbose


# Combining tables

qiime feature-table merge \
    --i-tables table1.qza \
    --i-tables table2.qza \
    --o-merged-table merged-table.qza


# Combining rep-seqs

qiime feature-table merge-seqs \
        --i-data rep-seqs1.qza \
        --i-data rep-seqs2.qza \
        --o-merged-data merged-rep-seqs.qza


# Taxonomy analysis

qiime feature-classifier classify-sklearn  \
--i-classifier silva-138-99-515-806-nb-classifier.qza  \
--i-reads merged-rep-seqs.qza  \
--o-classification taxonomy.qza


# creating a barplot

qiime taxa barplot   \
--i-table merged-table.qza   \
--i-taxonomy taxonomy.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization taxa-bar-plots.qzv


# Creating phylogenetic trees

qiime phylogeny align-to-tree-mafft-fasttree   \
--i-sequences merged-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


# alpha rarefaction with goods_coverage 

qiime diversity alpha-rarefaction   \
--i-table merged-table.qza   \
--p-metrics 'goods_coverage' \
--i-phylogeny rooted-tree.qza   \
--p-iterations 1 \
--p-max-depth 100000   \
--p-steps 1000   \
--m-metadata-file Metadata_202306.tsv \
--o-visualization alpha-rarefaction_good.qzv


# Calculation of diversity index.　

qiime diversity core-metrics-phylogenetic   \
--i-phylogeny rooted-tree.qza   \
--i-table merged-table.qza   \
--p-sampling-depth 4305   \
--m-metadata-file Metadata_202306.tsv   \
--output-dir core-metrics-results


# Calculation of diversity index (chao1, Simpson, Simpson_e)

qiime diversity alpha   \
--i-table rarefied_table.qza  \
--p-metric 'chao1' \
--o-alpha-diversity chao1_vector.qza

qiime diversity alpha   \
--i-table rarefied_table.qza  \
--p-metric 'simpson' \
--o-alpha-diversity simpson_vector.qza

qiime diversity alpha   \
--i-tablerarefied_table.qza   \
--p-metric 'simpson_e' \
--o-alpha-diversity simpson_e_vector.qza


# Comparison of alpha diversity  

qiime diversity alpha-group-significance   \
--i-alpha-diversity shannon_vector.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualizationshannon_significance.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity observed_features_vector.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization observed_features_significance.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity chao1_vector.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization chao1_significance.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity evenness_vector.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization evenness_significance.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity faith_pd_vector.qza   \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization faith_pd_significance.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity simpson_e_vector.qza  \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization simpson_e_vector.qzv

qiime diversity alpha-group-significance   \
--i-alpha-diversity simpson_vector.qza  \
--m-metadata-file Metadata_202306.tsv   \
--o-visualization simpson_vector.qzv


# Quantification of the alpha diversity index.

qiime metadata tabulate --m-input-file shannon_vector.qza --o-visualization shannon.qzv

qiime metadata tabulate --m-input-file chao1_vector.qza --o-visualization chao1.qzv

qiime metadata tabulate --m-input-file simpson_e_vector.qza --o-visualization simpson_e.qzv

qiime metadata tabulate --m-input-file simpson_vector.qza --o-visualization simpson.qzv

qiime metadata tabulate --m-input-file observed_features_vector.qza --o-visualization observed_features.qzv

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

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


# Calculation of β-diversity (comparing similarity of gut microbiota between mitochondrial lineages).

qiime diversity beta-group-significance   \
--i-distance-matrix jaccard_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column lineage   \
--o-visualization lineage_jaccard_distance_significance.qzv   \
--p-pairwise

qiime diversity beta-group-significance   \
--i-distance-matrix bray_curtis_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column lineage   \
--o-visualization lineage_bray_curtis_distance_significance.qzv   \
--p-pairwise	

qiime diversity beta-group-significance   \
--i-distance-matrix unweighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column lineage   \
--o-visualization lineage_unweighted_unifrac_significance.qzv   \
--p-pairwise

qiime diversity beta-group-significance   \
--i-distance-matrix weighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column lineage   \
--o-visualization lineage_weighted_unifrac_significance.qzv   \
--p-pairwise


# Calculation of β-diversity (comparing similarity of gut microbiota between management groups).

qiime diversity beta-group-significance   \
--i-distance-matrix jaccard_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column data5   \
--o-visualization origin_jaccard_distance_significance.qzv   \
--p-pairwise

qiime diversity beta-group-significance   \
--i-distance-matrix bray_curtis_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column data5   \
--o-visualization origin_bray_curtis_distance_significance.qzv   \
--p-pairwise	

qiime diversity beta-group-significance   \
--i-distance-matrix unweighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column data5   \
--o-visualization origin_unweighted_unifrac_significance.qzv   \
--p-pairwise

qiime diversity beta-group-significance   \
--i-distance-matrix weighted_unifrac_distance_matrix.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column data5   \
--o-visualization origin_weighted_unifrac_significance.qzv   \
--p-pairwise


#
# [ANCOM (Investigating differences in the gut bacteria between mitochondrial lineages)]
# 

# sample extraction

qiime feature-table filter-samples   \
--i-table rarefied_table.qza   \
--m-metadata-file Metadata_202306.tsv   \
--p-where "[lineage]='C' OR [lineage]='N2' OR [lineage]='S'"  \
--o-filtered-table lineage_table.qza


# Compared with level6 (genus level)

qiime taxa collapse   \
--i-table lineage_table.qza  \
--i-taxonomy taxonomy.qza   \
--p-level 6   \
--o-collapsed-table lineage_L6_table.qza


# pseudocount+1

qiime composition add-pseudocount   \
--i-table lineage_L6_table.qza   \
--o-composition-table lineage_L6_comp_table.qza


# ANCOM

qiime composition ancom   \
--i-table lineage_L6_comp_table.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column lineage   \
--o-visualization linegae_L6_ancom.qzv


#
# [ANCOM (Investigating differences in the gut bacteria between management groups)]
#

# sample extraction

qiime feature-table filter-samples   \
--i-table rarefied_table.qza   \
--m-metadata-file Metadata_202306.tsv   \
--p-where "[data5]='E' OR [data5]='O'"  \
--o-filtered-table origin_table.qza


# Compared with level6 (genus level)

qiime taxa collapse   \
--i-table origin_table.qza  \
--i-taxonomy taxonomy.qza   \
--p-level 6   \
--o-collapsed-table origin_L6_table.qza


# pseudocount+1

qiime composition add-pseudocount   \
--i-table origin_L6_table.qza   \
--o-composition-table origin_L6_comp_table.qza


# ANCOM

qiime composition ancom   \
--i-table origin_L6_comp_table.qza   \
--m-metadata-file Metadata_202306.tsv   \
--m-metadata-column data5   \
--o-visualization origin_L6_ancom.qzv


# Decomposition of the feature table by classification hierarchy (levels 2-7)

for i in {2..7}; do \
qiime taxa collapse \
--i-table rarefied_table.qza \
--i-taxonomy taxonomy.qza \
--p-level $i \
--o-collapsed-table table-l${i}.qza \
;done


# Conversion to occupancy.

for i in {2..7}; do \
qiime feature-table relative-frequency \
--i-table table-l${i}.qza \
--o-relative-frequency-table table-l${i}._relabn.qza \
;done


# Visualization

for i in {2..7}; do \
qiime metadata tabulate \
--m-input-file table-l${i}._relabn.qza \
--o-visualization table-l${i}._relabn.qzv \
;done