My study is comparing the fungal diversity between soil that experienced a wildfire (NoRxF) between soil that was treated with prescribed fire prior to experiencing the wildfire (RxF). It's a paired study design where 6 sites from each "treatment" are paired for a total of 12 sites. I also sampled at two different time-points, 6-months and 12-months post-wildfire, as well as took 9 soil samples from each site using a 3x3m grid for a total of 216 samples (108 per time point). Here is a picture that should help communicate the study design:
My hypothesis is that the soil treated with prescribed fire will have a higher fungal diversity than the soil that only experienced the wildfire. To test that hypothesis, I'm looking to compare each of the paired sites. I also want to compare the sites between time points to see how the communities, and their differences, have changed.
I've read the ANCOM-BC2 paper that shows that ANCOM-BC2 is appropriate for categorical treatments, and I understand that Qiime2 uses the conservative sensitivity score (SS) version. I see that folks use the formula parameters to do paired comparisons, but I have a column called "Sample-Time" (where samples from site 5 collected in January are '5-1' while samples collected in July are '5-7') that I'm using to group the data:
## Calculate differential abundance
qiime composition ancombc2 \
--i-table table_decont_pvlnc_unclst_noctl.qza \
--m-metadata-file metadata_noctl.txt \
--p-group 'Sample-Time' \
--p-fixed-effects-formula 'Treatment' \
--o-ancombc2-output ancom_pdecon_unclst.qza
As for the actual comparison between time points, a form moderator recommended that this user instead consider using Gemelli and construct an additive log ratio (ALR) using the factor loadings in order to account for the effect of each individual subject (for me, site).
One thing I should add is that my sampling depth varies a lot, and while the alpha rarefaction plotting showed a leveling out of Shannon diversity at 1556 and observed features at 6222, sampling at a depth of 6222 means I retain only 4.03% observations in 75.80% samples. I believe this means I'm losing A LOT of data, which makes metrics that don't need rarefaction (like differential abundance) more attractive.
Here is my data after filtering out sequences with <500 reads, extracting the ITS1 region with ITSxpress, denoinsing with Dada2, assigning taxonomy using Colin's 99% pretrained UNITE classifier, and decontaminating with Decontam:
asvs_decont_pvlnc_unclst_sgl.qzv (1.2 MB)