Install and test qiime2-docker

#Cannot run under QIIME1, switch to QIIME2: pick_open_reference_otus.py -r/home/jhuang/REFs/SILVA_132_QIIME_release/rep_set/rep_set_16S_only/99/silva_132_99_16S.fna -i test.fna -o clustering_test/ -p clustering_params.txt --parallel --verbose

docker pull quay.io/qiime2/core:2023.9

docker run -it --rm \
-v /mnt/md1/DATA/Data_Marius_16S_2025:/data \
-v /home/jhuang/REFs:/home/jhuang/REFs \
quay.io/qiime2/core:2023.9 bash
cd /data

Import the fastq-files to paired-end-demux.qza

#https://docs.qiime2.org/2018.8/tutorials/importing/

qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path pe-33-manifest --output-path paired-end-demux.qza --input-format PairedEndFastqManifestPhred33
#--> 1095204304 Mai 27 15:11 paired-end-demux.qza

qiime demux summarize \
--i-data paired-end-demux.qza \
--o-visualization demux_pe.qzv

#https://view.qiime2.org
#qiime tools view demux_pe.qzv

Optimizing the parameters trunc-len-f and trunc-len-r and denoising with DADA2: optimized parameters is f240_r240

#Your amplicon (V3–V4 region) is ~464 bp, so you need ≥20–30 bp overlap
#464-38=426; 440 is the longst +12 nt for overlapping=we need 452 nt!

#Optimize the parameters --p-trunc-len-f and --p-trunc-len-r
./dada2_batch_test.sh

        #!/bin/bash

        # Set your base inputs
        INPUT=paired-end-demux.qza
        TRIM_LEFT_F=17
        TRIM_LEFT_R=21

        # Output base
        OUTPUT_DIR=dada2_tests
        mkdir -p $OUTPUT_DIR

        # Loop over trunc-len-f and trunc-len-r combinations
        # Forward: from 220 to 240
        # Reverse: from 210 to 230
        i=1
        for f in 240 235 230 225; do
        for r in 225 220 215; do
            OUT=test_${i}_f${f}_r${r}
            echo "Running: $OUT"
            mkdir -p $OUTPUT_DIR/$OUT

            qiime dada2 denoise-paired \
            --i-demultiplexed-seqs $INPUT \
            --p-trim-left-f $TRIM_LEFT_F \
            --p-trim-left-r $TRIM_LEFT_R \
            --p-max-ee-f 3 --p-max-ee-r 5 \
            --p-trunc-len-f $f \
            --p-trunc-len-r $r \
            --p-n-threads 32 \
            --o-table $OUTPUT_DIR/$OUT/table.qza \
            --o-representative-sequences $OUTPUT_DIR/$OUT/rep-seqs.qza \
            --o-denoising-stats $OUTPUT_DIR/$OUT/denoising-stats.qza \
            --verbose > $OUTPUT_DIR/$OUT/log.txt 2>&1

            ((i++))
        done
        done

for f in dada2_tests2/test_*/denoising-stats.qza; do
qiime metadata tabulate \
    --m-input-file $f \
    --o-visualization ${f%.qza}.qzv
done

#pandaseq.out: grep ">" A1_R1.fastq.gz_merged.fasta | wc -l #8229;  grep ">" A10_R1.fastq.gz_merged.fasta | wc -l #9165

# The best choice is f251_r251, since the first 17 and 21 bases with bad quality are anyway removed!
#f251_r251: sample-A1   18299   10989   60.05   10609   7129    38.96   6535    35.71;  sample-A10  18736   12249   65.38   11778   7444    39.73   6978    37.24
#f251_r250: sample-A1   18299   11855   64.78   11435   7431    40.61   6823    37.29;  sample-A10  18736   13092   69.88   12590   7714    41.17   7206    38.46
#f251_r245: sample-A1   18299   12642   69.09   12180   7860    42.95   7193    39.31;  sample-A10  18736   13981   74.62   13457   8218    43.86   7649    40.83
#f251_r240: sample-A1   18299   12678   69.28   12214   8060    44.05   7387    40.37;  sample-A10  18736   14018   74.82   13498   8412    44.9    7758    41.41

#f250_r240: sample-A1   18299   13705   74.89   13191   8796    48.07   7984    43.63
#f250_r235: sample-A1   18299   13716   74.95   13198   8793    48.05   7969    43.55
#f250_r230: sample-A1   18299   13739   75.08   13217   9023    49.31   8113    44.34;  sample-A10  18736   14838   79.2    14159   8895    47.48   8101    43.24
#f245_r240: sample-A1   18299   14513   79.31   14019   9472    51.76   8739    47.76;  sample-A10  18736   15609   83.31   15102   9605    51.26   8880    47.4
#f245_r235: sample-A1   18299   14524   79.37   14026   9485    51.83   8746    47.79;  sample-A10  18736   15634   83.44   15127   9685    51.69   8869    47.34
#f245_r230: sample-A1   18299   14547   79.5    14045   8845    48.34   8058    44.04;  sample-A10  18736   15664   83.6    15156   8812    47.03   7999    42.69
#f240_r240: sample-A1   18299   14647   80.04   14164   9869    53.93   8932    48.81;  sample-A10  18736   15728   83.95   15213   10488   55.98   9081    48.47 *
#f240_r235: sample-A1   18299   14661   80.12   14172   9125    49.87   8194    44.78;  sample-A10  18736   15755   84.09   15242   9579    51.13   8105    43.26
#f240_r230: sample-A1   18299   14686   80.26   14191   4952    27.06   4666    25.5;   sample-A10  18736   15785   84.25   15267   3489    18.62   3341    17.83

#f240_r225: sample-A1   18299   14701   80.34   14206   4575    25  4297    23.48
#f240_r220: sample-A1   18299   14720   80.44   14223   4588    25.07   4310    23.55
#f240_r225: sample-A1   18299   14747   80.59   14250   3976    21.73   3758    20.54
#f230_r220: sample-A1   18299   14972   81.82   14514   3   0.02    3   0.02

#The output of the optimized denoising: (*) dada2_tests2/test_7_f240_r240/table.qza and dada2_tests2/test_7_f240_r240/rep-seqs.qza.

Visualize outputs

qiime feature-table summarize \
--i-table dada2_tests2/test_7_f240_r240/table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file qiime2_metadata.tsv

#Table summary
#Metric Sample
#Number of samples  137
#Number of features 3,039
#Total frequency    1,641,484
#
#Frequency per sample
#Minimum frequency  413.0
#1st quartile   10,319.0
#Median frequency   11,530.0
#3rd quartile   13,146.0
#Maximum frequency  40,022.0
#Mean frequency 11,981.635036496351
#
#Frequency per feature
#Minimum frequency  1.0
#1st quartile   3.0
#Median frequency   8.0
#3rd quartile   95.5
#Maximum frequency  56,472.0
#Mean frequency 540.1395195788089

#qiime tools peek table.qza
#qiime tools peek qiime2_metadata.tsv

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

qiime metadata tabulate \
--m-input-file dada2_tests2/test_7_f240_r240/denoising-stats.qza \
--o-visualization denoising-stats.qzv

Import reference sequences and taxonomy (SILVA 132)

qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path /home/jhuang/REFs/SILVA_132_QIIME_release/rep_set/rep_set_16S_only/99/silva_132_99_16S.fna \
--output-path silva_132_99_otus.qza \
--input-format DNAFASTAFormat

qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path /home/jhuang/REFs/SILVA_132_QIIME_release/taxonomy/16S_only/99/consensus_taxonomy_7_levels.txt \
--output-path silva_132_99_taxonomy.qza

Assign taxonomy

qiime feature-classifier classify-consensus-vsearch \
--i-query dada2_tests2/test_7_f240_r240/rep-seqs.qza \
--i-reference-reads silva_132_99_otus.qza \
--i-reference-taxonomy silva_132_99_taxonomy.qza \
--p-perc-identity 0.97 \
--p-threads 64 \
--o-classification taxonomy.qza \
--o-search-results search-results.qza

Visualize taxonomy

qiime taxa barplot \
--i-table dada2_tests2/test_7_f240_r240/table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file qiime2_metadata.tsv \
--o-visualization taxa-bar-plots.qzv

Build phylogenetic tree

qiime alignment mafft \
--i-sequences dada2_tests2/test_7_f240_r240/rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza

qiime alignment mask \
--i-alignment aligned-rep-seqs.qza \
--o-masked-alignment masked-aligned-rep-seqs.qza

qiime phylogeny fasttree \
--i-alignment masked-aligned-rep-seqs.qza \
--o-tree unrooted-tree.qza

(*) qiime phylogeny midpoint-root \
--i-tree unrooted-tree.qza \
--o-rooted-tree rooted-tree.qza

Core diversity analysis

#The -e 6389 flag sets the even sampling depth (rarefaction depth) to 6,389 reads for diversity analyses.
#All samples will be rarefied to 4,753 reads.
#Samples with fewer reads are excluded.
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table dada2_tests2/test_7_f240_r240/table.qza \
--p-sampling-depth 6389 \
--m-metadata-file qiime2_metadata.tsv \
--output-dir core_metrics_results

qiime diversity alpha \
--i-table table.qza \
--p-metric chao1 \
--o-alpha-diversity core_metrics_results/chao1_vector.qza

qiime tools export --input-path core_metrics_results/shannon_vector.qza --output-path exported_alpha/shannon
qiime tools export --input-path core_metrics_results/faith_pd_vector.qza --output-path exported_alpha/faith_pd
qiime tools export --input-path core_metrics_results/observed_features_vector.qza --output-path exported_alpha/observed_features
qiime tools export --input-path core_metrics_results/chao1_vector.qza --output-path exported_alpha/chao1

qiime tools export \
--input-path core_metrics_results/unweighted_unifrac_distance_matrix.qza \
--output-path exported_unweighted_unifrac
qiime tools export \
--input-path core_metrics_results/weighted_unifrac_distance_matrix.qza \
--output-path exported_weighted_unifrac

qiime diversity beta-group-significance \
--i-distance-matrix core_metrics_results/weighted_unifrac_distance_matrix.qza \
--m-metadata-file qiime2_metadata.tsv \
--m-metadata-column Group \
--p-pairwise \
--p-method permanova \
--o-visualization beta_group_significance.qzv

qiime tools export \
--input-path beta_group_significance.qzv \
--output-path exported_beta_group

#✅ Group 1 / Group 2 — The pairwise comparisons.
#✅ Sample size — Number of samples used in the test.
#✅ Permutations — Number of permutations in the PERMANOVA test.
#✅ pseudo-F — The test statistic from PERMANOVA.
#✅ p-value — The unadjusted p-value.
#✅ q-value — The adjusted p-value (Bonferroni in QIIME2).
#The q-value column (also sometimes called p-adj) is the multiple-testing corrected p-value.
#👉 q < 0.05 means the difference is statistically significant between those two groups, even after correction.
#“There is a significant difference in community composition between Group 1 and Group 2 (p=0.002).”

#The --p-sampling-depth 6389 parameter is directly equivalent to QIIME 1’s -e 6389!
#The QIIME 2 command will compute:
#✅ Alpha diversity metrics (Observed OTUs, Shannon, Faith PD, Evenness)
#✅ Beta diversity distance matrices (UniFrac, Bray-Curtis)
#✅ PCoA plots
#✅ Excludes samples with fewer than 6,389 reads.

#📦 Output Folders and Files
#Output Description
#table.qzv  Visual of feature table (sample counts)
#rep-seqs.qzv   Sequences per ASV
#denoising-stats.qzv    DADA2 read tracking
#taxonomy.qza/.qzv  Taxonomic classification of ASVs
#taxa-bar-plots.qzv Interactive bar plots
#core_metrics_results/  Alpha/beta diversity metrics + PCoA plots

Prepare three files feeding to Phyloseq.Rmd: table.qza (see above with ), rooted-tree.qza (see above with ), qiime2_metadata_for_qza_to_phyloseq.tsv edited from qiime2_metadata.tsv.

# Rarefying can be performed here, or in Phyloseq.Rmd (default), therefore, we don't need this step any more.
qiime feature-table summarize \
--i-table core_metrics_results/rarefied_table.qza \
--o-visualization rarefied_table.qzv \
--m-sample-metadata-file qiime2_metadata.tsv

#Table summary
#Metric Sample
#Number of samples  136
#Number of features 2,781
#Total frequency    868,904

# In QIIME2, we need table.qza, not biom-file, therefore, we don't need this step any more.
qiime tools export \
--input-path core_metrics_results/rarefied_table.qza \
--output-path exported_rarefied_table
#-->
exported_rarefied_table/feature-table.biom

biom convert \
-i exported_rarefied_table/feature-table.biom \
-o exported_rarefied_table/feature-table.tsv \
--to-tsv

#✅ Old QIIME 1 table with GenBank IDs (like EF603722.1.1487) as feature labels.
#✅ QIIME 2 table where feature IDs are hashes (like 0b438323a296b5f2ce2c8bbe3949ee8d).

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

#Feature ID    Taxon                               Confidence
#0b4383...     k__Bacteria; p__Proteobacteria...   0.98
#dfa833...     k__Bacteria; p__Firmicutes...       0.87
#...

# ---- I used the following to generate two file for feeding in the Phyloseq.Rmd ----

#-1- exported_table/feature-table.biom corresesponds to table_even6389.biom in QIIME1, but in QIIME2, we don't need biom-file, instead of table.qza.

qiime tools export \
--input-path dada2_tests2/test_7_f240_r240/table.qza \
--output-path exported_table
#--> exported_table/feature-table.biom

#-2- exported-tree/tree.nwk corresesponds to rep_set.tre in QIIME1

qiime tools export \
--input-path rooted-tree.qza \
--output-path exported-tree
#--> exported-tree/tree.nwk

# ---- The code in Phyloseq.Rmd ----

#install.packages("remotes")
#remotes::install_github("jbisanz/qiime2R")
#"core_metrics_results/rarefied_table.qza", rarefying performed in the code, therefore import the raw table.
library(qiime2R)
ps.ng.tax <- qza_to_phyloseq(
    features =  "dada2_tests2/test_7_f240_r240/table.qza",
    tree = "rooted-tree.qza",
    metadata = "qiime2_metadata_for_qza_to_phyloseq.tsv"
)
# or
#biom convert \
#      -i ./exported_table/feature-table.biom \
#      -o ./exported_table/feature-table-v1.biom \
#      --to-json
#ps.ng.tax <- import_biom("./exported_table/feature-table-v1.biom", treefilename="./exported-tree/tree.nwk")

#Note that the alpha- and beta-diversity-files needed in Phyloseq.Rmd has been prepared in the step 9.

Figures generated by Phyloseq.Rmd and MicrobiotaProcess_*.R

The following files can be found under server.

./Phyloseq.Rmd (Result Phyloseq.html)
./MicrobiotaProcess_cluster1_Group9-11_vs_cluster2_Group12-14_orig.R
./MicrobiotaProcess_Group1_vs_Group2.R
./MicrobiotaProcess_Group3_vs_Group4.R
./MicrobiotaProcess_PCA_Group1-4.R
./MicrobiotaProcess_PCA_Group9-14.R