Here, we analyse abundances with three different methods: Wilcoxon test (CLR), DESeq2,
and ANCOM-BC. All of these test statistical differences between groups.
We will analyse Genus level abundances.
We might want to first perform prevalence filtering to reduce the amount of multiple tests. In this particular dataset, all genera pass a prevalence threshold of 10%, therefore, we do not perform filtering.
Wilcoxon test
A Wilcoxon test estimates the difference in an outcome between two groups. It is a
non-parametric alternative to a t-test, which means that the Wilcoxon test
does not make any assumptions about the data.
Let’s first combine the data for the testing purpose.
Now we can start with the Wilcoxon test. We test all the taxa by looping through columns,
and store individual p-values to a vector. Then we create a data frame from collected
data.
The code below does the Wilcoxon test only for columns that contain abundances,
not for columns that contain patient status.
Multiple tests were performed. These are not independent, so we need
to adjust p-values for multiple testing. Otherwise, we would increase
the chance of a type I error drastically depending on our p-value
threshold. By applying a p-value adjustment, we can keep the false
positive rate at a level that is acceptable. What is acceptable
depends on our research goals. Here we use the fdr method, but there
are several other methods as well.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
DESeq2
Our second analysis method is DESeq2. This method performs the data
normalization automatically. It also takes care of the p-value
adjustment, so we don’t have to worry about that.
DESeq2 utilizes a negative binomial distribution to detect differences in
read counts between groups. Its normalization takes care of the
differences between library sizes and compositions. DESeq2 analysis
includes multiple steps, but they are done automatically. More
information can be found, e.g., from Harvard Chan Bioinformatic Core’s
tutorial Introduction to DGE –
ARCHIVED
Now let us show how to do this. First, run the DESeq2 analysis.
## converting counts to integer mode
## Warning in DESeqDataSet(tse_genus, ~patient_status): 2 duplicate rownames were
## renamed by adding numbers
## Warning in DESeqDataSet(tse_genus, ~patient_status): some variables in design
## formula are characters, converting to factors
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 11 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | taxon | |
---|---|---|---|---|---|---|---|
Genus:Ruminococcaceae_UCG-014 | 22.548297 | -24.891268 | 2.460684 | -10.115589 | 0.0000000 | 0.0000000 | Genus:Ruminococcaceae_UCG-014 |
Order:Bacteroidales | 40.353733 | -9.241798 | 2.136205 | -4.326270 | 0.0000152 | 0.0002730 | Order:Bacteroidales |
Genus:Faecalibacterium | 231.079502 | -7.074433 | 1.745612 | -4.052694 | 0.0000506 | 0.0006835 | Genus:Faecalibacterium |
Genus:Catabacter | 18.045614 | -6.615454 | 1.716150 | -3.854823 | 0.0001158 | 0.0012508 | Genus:Catabacter |
Genus:Butyricicoccus | 2.392885 | -5.179608 | 2.948055 | -1.756957 | 0.0789251 | 0.3278426 | Genus:Butyricicoccus |
Order:Gastranaerophilales | 2.067972 | -3.054975 | 2.938641 | -1.039588 | 0.2985315 | 0.7269742 | Order:Gastranaerophilales |
ANCOM-BC
The analysis of composition of microbiomes with bias correction (ANCOM-BC)
is a recently developed method for differential abundance testing. It is based on an
earlier published approach.
The former version of this method could be recommended as part of several approaches:
A recent study
compared several mainstream methods and found that among another method, ANCOM produced the most consistent results and is probably a conservative approach. Please note that based on this and other comparisons, no single method can be recommended across all datasets. Rather, it could be recommended to apply several methods and look at the overlap/differences.
As the only method, ANCOM-BC incorporates the so called sampling fraction into the model. The latter term could be empirically estimated by the ratio of the library size to the microbial load. Variations in this sampling fraction would bias differential abundance analyses if ignored. Furthermore, this method provides p-values, and confidence intervals for each taxon.
It also controls the FDR and it is computationally simple to implement.
As we will see below, to obtain results, all that is needed is to pass
a phyloseq object to the ancombc()
function. Therefore, below we first convert
our tse
object to a phyloseq
object. Then, we specify the formula. In this formula, other covariates could potentially be included to adjust for confounding.
Please check the function documentation
to learn about the additional arguments that we specify below. Also, see here for another example for more than 1 group comparison.
The object out
contains all relevant information. Again, see the
documentation of the function
under Value for an explanation of all the output objects. Our question can be answered
by looking at the res
object, which now contains dataframes with the coefficients,
standard errors, p-values and q-values. Conveniently, there is a dataframe diff_abn
.
Here, we can find all differentially abundant taxa. Below we show the first 6 entries of this dataframe:
patient_statusControl | |
---|---|
172647198 | FALSE |
1726478 | FALSE |
172647201 | FALSE |
17264798 | FALSE |
172647195 | FALSE |
1726472 | FALSE |
In total, this method detects 14 differentially abundant taxa.
Comparison of the methods
Let’s compare results that we got from the methods.
As we can see from the scatter plot, DESeq2 gives lower p-values than Wilcoxon test.
Prints number of p-values under 0.05
## [1] "Wilcoxon test p-values under 0.05: 2/54"
## [1] "DESeq2 p-values under 0.05: 7/54"
## [1] "ANCOM p-values under 0.05: 14/49"
We can also look at the intersection of identified taxa
## [1] "Faecalibacterium" "[Ruminococcus]_gauvreauii_group"
Comparison of abundance
In previous steps, we got information which taxa vary between ADHD and control groups.
Let’s plot those taxa in the boxplot, and compare visually if abundances of those taxa
differ in ADHD and control samples. For comparison, let’s plot also taxa that do not
differ between ADHD and control groups.
Let’s first gather data about taxa that have highest p-values.
Next, let’s do the same but for taxa with lowest p-values.
Then we can plot these six different taxa. Let’s arrange them into the same picture.
We plotted those taxa that have the highest and lowest p values according to DESeq2. Can you create a plot that shows the difference in abundance in “[Ruminococcus]_gauvreauii_group”, which is the other taxon that was identified by all tools. Try for yourself! Below you find one way how to do it.
Read more here: Source link