Author: Ilze Brila
This is the code used in data processing and analysis in QIIME2 for manuscript Idiosyncratic effects of coinfection on the association between systemic pathogens and the gut microbiota of a wild rodent, the bank vole (Myodes glareolus)".
All analysis and processing up to differential abundance analysis done exactly as for a previous manuscript on the same data.
Raw data are deposited in GenBank under accession number PRJNA702897.
All analysis was done in Qiime2 v.2019.10.
Table of contents:
Raw reads were demultiplexed at the sequencing facility Institute for Molecular Medicine Finland, HiLIFE (FIMM). Illumina adapter sequences were removed using cutadapt v.1.10.
# import data
qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path manifest_gut.csv\
--input-format PairedEndFastqManifestPhred33\
--output-path gut-demux.qza
# look at the summary
qiime demux summarize \
--i-data gut-demux.qza \
--o-visualization gut-demux_sum.qzv
DADA2 plugin was used to:
# denoise data and trim
qiime dada2 denoise-paired \
--i-demultiplexed-seqs gut16S-demux.qza \
--p-trim-left-f 19 \
--p-trunc-len-f 250 \
--p-trim-left-r 20 \
--p-trunc-len-r 176 \
--p-n-threads $SLURM_CPUS_PER_TASK \
--o-table gut_feat_table_dada2.qza \
--o-representative-sequences gut_rep_seqs_dada2.qza \
--o-denoising-stats gut_dada2_stats.qza
# summarize and examine the output
qiime feature-table summarize \
--i-table gut_feat_table_dada2.qza \
--m-sample-metadata-file metadata_qiime_guts.txt \
--o-visualization gut_feat_table_dada2.qzv
qiime feature-table tabulate-seqs \
--i-data gut_rep_seqs_dada2.qza \
--o-visualization gut_rep_seqs_dada2.qzv
qiime metadata tabulate \
--m-input-file gut_dada2_stats.qza \
--o-visualization gut_dada2_stats.qzv
Rare ASVs (ASVs with <10 reads across all samples) were removed, as these might be potential errors or artifacts.
# filter the features
qiime feature-table filter-features \
--i-table gut_feat_table_dada2.qza \
--p-min-frequency 10 \
--o-filtered-table gut_feat_table_filt10.qza
# visualize and summarize new filtered feature table
qiime feature-table summarize \
--i-table gut_feat_table_filt10.qza \
--m-sample-metadata-file metadata_qiime_gut.txt \
--o-visualization gut_feat_table_filt10.qzv
# filter representative sequences with the same parameters
qiime feature-table filter-seqs \
--i-data gut_rep_seqs_dada2.qza \
--i-table gut_feat_table_filt10.qza \
--o-filtered-data gut_rep_seqs_filt10.qza
Taxonomy to the representative sequences was assigned using the SILVA 132 database.
# create taxonomy
qiime feature-classifier classify-sklearn \
--i-reads gut_rep_seqs_filt10.qza \
--i-classifier silva-132-99-515-806-nb-classifier.qza \
--o-classification gut_taxonomy_filt10_silva.qza
# summarize the taxonomy
qiime metadata tabulate \
--m-input-file gut_taxonomy_filt10_silva.qza \
--o-visualization gut_taxonomy_filt10_silva.qzv
Removal of ASVs that could not be assigned to a bacterial phylum, and the ASVs classified as mitochondria or chloroplasts.
# filter the feature table
qiime taxa filter-table \
--i-table gut_feat_table_filt10.qza \
--i-taxonomy gut_taxonomy_filt10_silva.qza \
--p-include "D_0__Bacteria;D_1__" \
--p-exclude mitochondria,chloroplast \
--o-filtered-table gut_feat_table_filt10_taxafilt.qza
# summarize the feature table
qiime feature-table summarize \
--i-table gut_feat_table_filt10_taxafilt.qza \
--m-sample-metadata-file metadata_qiime_gut.txt \
--o-visualization gut_feat_table_filt10_taxafilt.qzv
# filter representative sequences
qiime taxa filter-seqs \
--i-sequences gut_rep_seqs_filt10.qza \
--i-taxonomy gut_taxonomy_filt10_silva.qza \
--p-include "D_0__Bacteria;D_1__" \
--p-exclude mitochondria,chloroplast \
--o-filtered-sequences gut_rep_seqs_filt10_taxafilt.qza
# summarize filtered representative sequences
qiime feature-table tabulate-seqs \
--i-data gut_rep_seqs_filt10_taxafilt.qza \
--o-visualization gut_rep_seqs_filt10_taxafilt.qzv
# make a barplot
qiime taxa barplot \
--i-table gut_feat_table_filt10_taxafilt.qza \
--i-taxonomy gut_taxonomy_filt10_silva.qza \
--m-metadata-file metadata_qiime_gut.txt \
--o-visualization gut-barplot-filt10-taxafilt.qzv
Creation of a mid-point rooted phylogenetic tree using FastTree.
# align representative sequences
qiime alignment mafft \
--i-sequences gut_rep_seqs_filt10_taxafilt.qza \
--o-alignment aligned_gut_rep_seqs_filt10_taxafilt.qza
# mask ambiguous sequences
qiime alignment mask \
--i-alignment aligned_gut_rep_seqs_filt10_taxafilt.qza \
--o-masked-alignment aligned_masked_gut_rep_seqs_filt10_taxafilt.qza
# create a phylogenetic tree
qiime phylogeny fasttree \
--i-alignment aligned_masked_gut_rep_seqs_filt10_taxafilt.qza \
--o-tree unrooted_gut_tree_filt10_taxafilt.qza
# root the tree
qiime phylogeny midpoint-root \
--i-tree unrooted_gut_tree_filt10_taxafilt.qza \
--o-rooted-tree rooted_gut_tree_filt10_taxafilt.qza
This was done to adjust for differences in sequencing depth across samples.
# look at raefaction curves using median feature frequency - with metadata
qiime diversity alpha-rarefaction \
--i-table gut_feat_table_filt10_taxafilt.qza \
--i-phylogeny rooted_gut_tree_filt10_taxafilt.qza \
--p-max-depth 54119 \
--p-metrics faith_pd \
--p-metrics shannon \
--p-metrics observed_otus \
--m-metadata-file metadata_qiime_gut.txt \
--o-visualization gut_a_rare_med_meta.qzv
# look at raefaction curves using median feature frequency - without metadata
qiime diversity alpha-rarefaction \
--i-table gut_feat_table_filt10_taxafilt.qza \
--i-phylogeny rooted_gut_tree_filt10_taxafilt.qza \
--p-max-depth 54119 \
--p-metrics faith_pd \
--p-metrics shannon \
--p-metrics observed_otus \
--o-visualization gut_a_rare_med.qzv
Most rarefaction curves level off at depth of 27'000 reads, with loss of 1 sample.
# rarefy
qiime feature-table rarefy \
--i-table gut_feat_table_filt10_taxafilt.qza \
--p-sampling-depth 27000 \
--o-rarefied-table gut_feat_table_rarefied.qza
# look at the rarefied feature table
qiime feature-table summarize \
--i-table gut_feat_table_rarefied.qza \
--m-sample-metadata-file metadata_qiime_gut.txt \
--o-visualization gut_feat_table_rarefied.qzv
# make a new barplot with rarefied data
qiime taxa barplot \
--i-table gut_feat_table_rarefied.qza \
--i-taxonomy gut_taxonomy_filt10_silva.qza \
--m-metadata-file metadata_qiime_gut.txt \
--o-visualization gut_barplot_rarefied.qzv
Differential abundnace analysis was done using Songbird v.1.0.4 plugin. Songbird is a compositionally aware method, as such unrarefied data were used. The model optimization was done according to Songbird tutorial on github . For brevity only final model code included.
As the Q2 values indicated low predictive power of Songbird models, we compared each model to three models where the respective explanatory variables were shuffled among the animals, to verify that model fit is better than random.
Pruning the feature table to include only non-gravid animals with full metadata available (n=179):
# prune the non-rarefied feature table
qiime feature-table filter-samples \
--i-table gut_feat_table_filt10_taxafilt.qza \
--m-metadata-file meta_179_jul2020.txt \
--o-filtered-table gut_feat_table_filt10_taxafilt_179.qza
# summarize the new feature table
qiime feature-table summarize \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--o-visualization gut_feat_table_filt10_taxafilt_179.qzv \
--m-sample-metadata-file meta_179_jul2020.txt
Running the full, null and baseline(covariates only) models on a dataset with 173 animals.
# full model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173.txt \
--p-formula "C(infected, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1 \
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials results_inf/inf_full_18_diff.qza \
--o-regression-stats results_inf/inf_full_18_stats.qza \
--o-regression-biplot results_inf/inf_full_18_biplot.qza
# look at model statistics
qiime songbird summarize-single --i-regression-stats inf_full_18_stats.qza --o-visualization inf_full_18_stats.qzv
# null model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173.txt \
--p-formula "1" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1 \
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials inf_null_18_diff.qza \
--o-regression-stats inf_null_18_stats.qza \
--o-regression-biplot inf_null_18_biplot.qza
# baseline model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173.txt \
--p-formula "poll_group + sex + BCI + head_w + location" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1 \
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials inf_base_18_diff.qza \
--o-regression-stats inf_base_18_stats.qza \
--o-regression-biplot inf_base_18_biplot.qza
# compare models
qiime songbird summarize-paired \
--i-regression-stats results_inf/inf_full_18_stats.qza \
--i-baseline-stats inf_null_18_stats.qza \
--o-visualization paired_stats_inf_18_full-null.qzv
qiime songbird summarize-paired \
--i-regression-stats results_inf/inf_full_18_stats.qza \
--i-baseline-stats inf_base_18_stats.qza \
--o-visualization paired_stats_inf_18_full-base.qzv
Compare to three random models:
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(inf_rand1, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1 \
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials results_inf/inf_rand_18_1_diff.qza \
--o-regression-stats results_inf/inf_rand_18_1_stats.qza \
--o-regression-biplot results_inf/inf_rand_18_1_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(inf_rand2, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1 \
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials results_inf/inf_rand_18_2_diff.qza \
--o-regression-stats results_inf/inf_rand_18_2_stats.qza \
--o-regression-biplot results_inf/inf_rand_18_2_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(inf_rand3, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 2 \
--p-epochs 20000 \
--p-summary-interval 1
--p-num-random-test-examples 17 \
--p-learning-rate 0.001 \
--o-differentials results_inf/inf_rand_18_3_diff.qza \
--o-regression-stats results_inf/inf_rand_18_3_stats.qza \
--o-regression-biplot results_inf/inf_rand_18_3_biplot.qza
# compare models to the full model
qiime songbird summarize-paired \
--i-regression-stats results_inf/inf_rand_18_1_stats.qza \
--i-baseline-stats results_inf/inf_full_18_stats.qza \
--o-visualization results_inf/paired_stats_inf_18_full-rand1.qzv
qiime songbird summarize-paired \
--i-regression-stats results_inf/inf_rand_18_2_stats.qza \
--i-baseline-stats results_inf/inf_full_18_stats.qza \
--o-visualization results_inf/paired_stats_inf_18_full-rand2.qzv
qiime songbird summarize-paired \
--i-regression-stats results_inf/inf_rand_18_3_stats.qza \
--i-baseline-stats results_inf/inf_full_18_stats.qza \
--o-visualization results_inf/paired_stats_inf_18_full-rand3.qzv
Running the full, null and baseline(covariates only) models on a dataset with 173 animals.
# full model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "C(coinf_stat, Treatment(0)) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_full_11_diff.qza \
--o-regression-stats results_coinf/coinf_full_11_stats.qza \
--o-regression-biplot results_coinf/coinf_full_11_biplot.qza \
# baseline model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "poll_group + sex + BCI + head_w + location" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_base_11_diff.qza \
--o-regression-stats results_coinf/coinf_base_11_stats.qza \
--o-regression-biplot results_coinf/coinf_base_11_biplot.qza
# null model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "1" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_null_11_diff.qza \
--o-regression-stats results_coinf/coinf_null_11_stats.qza \
--o-regression-biplot results_coinf/coinf_null_11_biplot.qza
# compare models
qiime songbird summarize-paired \
--i-regression-stats coinf_full_11_stats.qza \
--i-baseline-stats coinf_base_11_stats.qza \
--o-visualization paired_stats_coinf_11_base-full.qzv
qiime songbird summarize-paired \
--i-regression-stats coinf_full_11_stats.qza \
--i-baseline-stats coinf_null_11_stats.qza \
--o-visualization paired_stats_coinf_11_null-full.qzv
Compare to three random models:
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(coinf_rand1, Treatment(0)) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_rand_11_r1_diff.qza \
--o-regression-stats results_coinf/coinf_rand_11_r1_stats.qza \
--o-regression-biplot results_coinf/coinf_rand_11_r1_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(coinf_rand2, Treatment(0)) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_rand_11_r2_diff.qza \
--o-regression-stats results_coinf/coinf_rand_11_r2_stats.qza\
--o-regression-biplot results_coinf/coinf_rand_11_r2_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(coinf_rand3, Treatment(0)) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 1 \
--p-epochs 15000 \
--p-summary-interval 1 \
--p-training-column test_coinf_2 \
--p-learning-rate 0.001 \
--o-differentials results_coinf/coinf_rand_11_r3_diff.qza \
--o-regression-stats results_coinf/coinf_rand_11_r3_stats.qza \
--o-regression-biplot results_coinf/coinf_rand_11_r3_biplot.qza
# compare to the full model
qiime songbird summarize-paired \
--i-regression-stats results_coinf/coinf_full_11_stats.qza \
--i-baseline-stats results_coinf/coinf_rand_11_r1_stats.qza \
--o-visualization results_coinf/paired_stats_coinf_11_full-rand1.qzv
qiime songbird summarize-paired \
--i-regression-stats results_coinf/coinf_full_11_stats.qza \
--i-baseline-stats results_coinf/coinf_rand_11_r2_stats.qza \
--o-visualization results_coinf/paired_stats_coinf_11_full-rand2.qzv
qiime songbird summarize-paired \
--i-regression-stats results_coinf/coinf_full_11_stats.qza \
--i-baseline-stats results_coinf/coinf_rand_11_r3_stats.qza \
--o-visualization results_coinf/paired_stats_coinf_11_full-rand3.qzv
# full model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "C(ana, Treatment('no')) + C(bab, Treatment('no')) + C(bor, Treatment('no')) + C(puu, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1\
--p-learning-rate 0.001 \
--o-differentials results_path/path_full_1_diff.qza \
--o-regression-stats results_path/path_full_1_stats.qza \
--o-regression-biplot results_path/path_full_1_biplot.qza
# baseline model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "C(poll_group + sex + BCI + head_w + location" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1 \
--p-learning-rate 0.001 \
--o-differentials results_path/path_base_1_diff.qza \
--o-regression-stats results_path/path_base_1_stats.qza \
--o-regression-biplot results_path/path_base_1_biplot.qza
# null model
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird.txt \
--p-formula "1" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1 \
--p-learning-rate 0.001 \
--o-differentials results_path/path_null_1_diff.qza \
--o-regression-stats results_path/path_null_1_stats.qza \
--o-regression-biplot results_path/path_null_1_biplot.qza
# compare models
qiime songbird summarize-paired \
--i-regression-stats path_full_1_stats.qza \
--i-baseline-stats path_null_1_stats.qza \
--o-visualization paired_stats_path_1_full-null.qzv
qiime songbird summarize-paired \
--i-regression-stats path_full_1_stats.qza \
--i-baseline-stats path_base_1_stats.qza \
--o-visualization paired_stats_path_1_full-base.qzv
Compare to three random models:
qiime songbird multinomial
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(ana_rand1, Treatment('no')) + C(bab_rand1, Treatment('no')) + C(bor_rand1, Treatment('no')) + C(puu_rand1, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1 \
--p-learning-rate 0.001 \
--o-differentials results_path/path_rand_1_r1_diff.qza \
--o-regression-stats results_path/path_rand_1_r1_stats.qza \
--o-regression-biplot results_path/path_rand_1_r1_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(ana_rand2, Treatment('no')) + C(bab_rand2, Treatment('no')) + C(bor_rand2, Treatment('no')) + C(puu_rand2, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1 \
--p-learning-rate 0.001 \
--o-differentials results_path/path_rand_1_r2_diff.qza \
--o-regression-stats results_path/path_rand_1_r2_stats.qza \
--o-regression-biplot results_path/path_rand_1_r2_biplot.qza
qiime songbird multinomial \
--i-table gut_feat_table_filt10_taxafilt_179.qza \
--m-metadata-file meta_173_songbird_randomR.txt \
--p-formula "C(ana_rand3, Treatment('no')) + C(bab_rand3, Treatment('no')) + C(bor_rand3, Treatment('no')) + C(puu_rand3, Treatment('no')) + poll_group + sex + BCI + head_w + location" \
--p-differential-prior 0.5 \
--p-epochs 10000 \
--p-summary-interval 1 \
--p-training-column test_path_1 \
--p-learning-rate 0.001 \
--o-differentials results_path/path_rand_1_r3_diff.qza \
--o-regression-stats results_path/path_rand_1_r3_stats.qza \
--o-regression-biplot results_path/path_rand_1_r3_biplot.qza
# compare to the full model
qiime songbird summarize-paired \
--i-regression-stats results_path/path_full_1_stats.qza \
--i-baseline-stats results_path/path_rand_1_r1_stats.qza \
--o-visualization results_path/paired_stats_path_1_full-rand1.qzv
qiime songbird summarize-paired \
--i-regression-stats results_path/path_full_1_stats.qza \
--i-baseline-stats results_path/path_rand_1_r2_stats.qza \
--o-visualization results_path/paired_stats_path_1_full-rand2.qzv
qiime songbird summarize-paired \
--i-regression-stats results_path/path_full_1_stats.qza \
--i-baseline-stats results_path/path_rand_1_r3_stats.qza \
--o-visualization results_path/paired_stats_path_1_full-rand3.qzv