skbio.stats.composition.ancombc2#
- skbio.stats.composition.ancombc2(table, metadata, formula, grouping=None, pseudocount=0, aggregator=None, var_quantile=0.05, max_iter=100, tol=1e-05, p_adjust='holm', alpha=0.05)[source]#
Perform differential abundance test using ANCOM-BC2.
Analysis of Compositions of Microbiomes with Bias Correction 2 (ANCOM-BC2) [1] extends ANCOM-BC by explicitly estimating sample-specific sampling fractions, correcting the transformed abundance data by these estimates, and refitting the model before inference.
Added in version 0.7.4.
- Parameters:
- tabletable_like of shape (n_samples, n_features)
Matrix containing count or proportional abundance data of the samples. See supported formats.
- metadatapd.DataFrame or 2-D array_like
Metadata of the samples. Rows correspond to samples and columns correspond to covariates (attributes). Must be a pandas DataFrame or convertible to a pandas DataFrame.
- formulastr or generic Formula object
Formula defining the model using factors included in the metadata columns. Refer to Patsy’s documentation on how to specify a formula.
- groupingstr, optional
Metadata column defining sample groups for post-hoc analyses. Must be a term in
formulaand contain at least three groups. Other terms are treated as adjustment covariates. It has no impact on the main result. If None (default), post-hoc analyses will be unavailable.- pseudocountint or float, optional
Pseudocount to add to all abundance data. Default is 0.
- aggregatorcallable, mapping, or 1-D array_like, optional
Rule for aggregating features before final regression. Can be a function or a dictionary that maps each feature ID to an aggregate ID, or a plain list or array of aggregate ID per feature in table order. By default, no aggregation is performed.
- var_quantilefloat, optional
Quantile of coefficient variances used as a regularization factor to reduce spurious significance from extremely small standard errors, particularly for rare features. Must be between 0 and 1. Set to 0 to disable. Default is 0.05.
- max_iterint, optional
Maximum number of iterations for the bias estimation process. Default is 100.
- tolfloat, optional
Absolute convergence tolerance for the bias estimation process. Default is 1e-5.
- alphafloat, optional
Significance level for the statistical tests. Must be in the range of (0, 1). Default is 0.05.
- p_adjuststr, optional
Method to correct p-values for multiple comparisons. Options are: Bonferroni (“bonf”/”bonferroni”), Holm-Bonferroni (“holm”/”holm-bonferroni”, default), Benjamini-Hochberg (“bh”/”benjamini-hochberg”), and Benjamini-Yekutieli (“by”/”benjamini-yekutieli”), or any method supported by statsmodels’
multipletestsfunction. Case-insensitive. If None, no correction will be performed.
- Returns:
ANCOMBCResultResult object with primary results and post-hoc analysis methods.
See also
ancombcANCOM-BC without explicit sampling-fraction correction.
struc_zeroStandalone structural zero detection.
Notes
This function is a Python re-implementation of the ANCOM-BC2 method [1], which was originally implemented in the R package
ANCOMBC. This function provides an efficient and scalable algorithm, with a simple interface consistent with other scikit-bio components. The output of this function should match that of the R package.Comparing with the R command
ancombc2, which completes the entire workflow in one go with extensive parameter settings, this function only executes the core ANCOM-BC2 algorithm and reports the primary results, while:Filtering of samples and features should be completed by the user prior to the function call.
Post-hoc analyses (global test, pairwise test, Dunnett’s test and trend test) are available as methods of the result object.
The structural zero test is available as a standalone function
struc_zero.Pseudocount sensitivity analysis is performed by repeating function calls with alternative
pseudocountsettings and merging the results.
A comprehensive tutorial on running the ANCOM-BC2 analysis and matching the R workflow is provided in the “Examples” section below.
References
Examples
>>> from skbio.stats.composition import ancombc2, struc_zero >>> import pandas as pd
Consider a dataset with 15 samples and eight features. The counts are sparse and skewed. Samples belong to three disease-status groups, while age is included as a potential confounder.
>>> samples = [f"S{i}" for i in range(1, 16)] >>> features = [f"F{i}" for i in range(1, 9)] >>> data = [ ... [1, 10, 16, 0, 0, 12, 0, 7], ... [0, 10, 45, 2, 1, 25, 1, 9], ... [4, 8, 22, 1, 5, 14, 1, 5], ... [1, 11, 22, 1, 4, 14, 1, 3], ... [4, 7, 36, 0, 0, 19, 0, 9], ... [2, 30, 26, 0, 0, 12, 1, 0], ... [1, 36, 10, 1, 5, 15, 2, 9], ... [1, 25, 11, 2, 2, 8, 0, 5], ... [1, 36, 23, 0, 2, 11, 1, 2], ... [3, 60, 33, 3, 4, 14, 1, 6], ... [3, 57, 30, 3, 1, 3, 0, 4], ... [0, 65, 23, 0, 2, 2, 0, 8], ... [2, 46, 19, 0, 1, 2, 0, 1], ... [0, 41, 26, 0, 2, 0, 0, 3], ... [3, 46, 20, 0, 0, 3, 0, 0], ... ] >>> table = pd.DataFrame(data, index=samples, columns=features) >>> status = ["mild"] * 5 + ["moderate"] * 5 + ["severe"] * 5 >>> age = [25, 48, 23, 35, 52, 51, 18, 29, 40, 29, 44, 39, 26, 37, 46] >>> metadata = pd.DataFrame({"status": status, "age": age}, index=samples)
Pre-processing
Before start, refer to the documentation of
ancombcon data pre-processing, such as filtering samples by total abundance, filtering features by prevalence, and defining the reference level of a categorical column in the metadata.There are two differences from
ancombc: First, apseudocountparameter is present inancombc2, defaulting to zero. There is no need to add a pseudocount to the data table prior to the analysis.Second, an
aggregatorparameter automates data aggregation within the function call. Aggregating data manually before the analysis will produce different output. This is becauseancombc2estimates initial model parameters on raw data before it aggregates them for final parameter estimation.Primary analysis
Fit disease status while adjusting for age. Specifying
grouping="status"also enables post-hoc analyses for this factor. The alphabetically first levelmildserves as the reference group. When done, display the primary result table.>>> res = ancombc2(table, metadata, formula="status + age", grouping="status") >>> res_main = res.result >>> res_main.round(3) Log(FC) SE W pvalue qvalue Signif FeatureID Covariate F1 Intercept 0.329 0.788 0.417 0.687 1.000 False status[T.moderate] -0.130 0.428 -0.305 0.768 1.000 False status[T.severe] 0.450 0.446 1.010 0.342 1.000 False age -0.004 0.021 -0.193 0.852 1.000 False F2 Intercept -0.097 0.514 -0.189 0.854 1.000 False status[T.moderate] 1.578 0.287 5.503 0.000 0.001 True status[T.severe] 1.880 0.300 6.272 0.000 0.000 True age -0.024 0.014 -1.794 0.100 0.703 False F3 Intercept 0.068 0.494 0.137 0.894 1.000 False status[T.moderate] -0.054 0.271 -0.200 0.845 1.000 False status[T.severe] -0.004 0.297 -0.014 0.989 1.000 False age 0.004 0.014 0.248 0.809 1.000 False F4 Intercept -0.990 0.781 -1.268 0.294 1.000 False status[T.moderate] 0.933 0.427 2.186 0.117 0.817 False status[T.severe] 0.600 0.516 1.164 0.329 1.000 False age 0.021 0.026 0.814 0.475 1.000 False F5 Intercept 2.376 0.663 3.586 0.009 0.071 False status[T.moderate] 0.069 0.388 0.178 0.864 1.000 False status[T.severe] -0.448 0.467 -0.960 0.369 1.000 False age -0.060 0.022 -2.798 0.027 0.213 False F6 Intercept 1.117 0.522 2.141 0.058 0.406 False status[T.moderate] -0.075 0.283 -0.264 0.797 1.000 False status[T.severe] -1.729 0.343 -5.042 0.001 0.004 True age -0.012 0.016 -0.729 0.483 1.000 False F7 Intercept 0.834 0.787 1.060 0.349 1.000 False status[T.moderate] 0.348 0.440 0.791 0.473 1.000 False status[T.severe] 0.000 0.516 0.000 1.000 1.000 False age -0.025 0.026 -0.957 0.393 1.000 False F8 Intercept 0.775 0.714 1.084 0.306 1.000 False status[T.moderate] 0.028 0.411 0.069 0.947 1.000 False status[T.severe] -0.530 0.444 -1.193 0.263 1.000 False age -0.013 0.021 -0.609 0.558 1.000 False
resis an instance ofANCOMBCResult. It behaves like its primary result tableres.resultfor display and column selection. Here we display only significant feature-covariate pairs:>>> res[res["Signif"]].round(3) Log(FC) SE W pvalue qvalue Signif FeatureID Covariate F2 status[T.moderate] 1.578 0.287 5.503 0.000 0.001 True status[T.severe] 1.880 0.300 6.272 0.000 0.000 True F6 status[T.severe] -1.729 0.343 -5.042 0.001 0.004 True
Log(FC)is the estimated coefficient on the natural-log scale. As shown, F2 is more abundant in both moderate and severe samples than in mild samples after adjusting for age, while F6 is less abundant in severe samples.Global test
The global test asks whether a feature differs between at least two of the three status groups, without identifying which groups differ.
>>> res_global = res.global_test() >>> res_global.round(3) W pvalue qvalue Signif FeatureID F1 1.786 0.457 1.000 False F2 60.287 0.000 0.000 True F3 0.040 0.078 0.469 False F4 4.780 0.233 1.000 False F5 1.385 0.623 1.000 False F6 25.906 0.000 0.002 True F7 0.865 0.975 1.000 False F8 1.589 0.513 1.000 False
This identifies F2 and F6 as globally differentially abundant.
Pairwise test
The pairwise test compares every pair of status groups. In addition to the two comparisons against the reference group (mild), it directly compares severe with moderate samples.
>>> res_pair = res.pairwise_test() >>> cols = ["Log(FC)", "qvalue", "Signif"] >>> res_pair.loc[res_pair["Signif"], cols].round(3) Log(FC) qvalue Signif FeatureID Comparison F2 status[T.moderate] 1.578 0.002 True status[T.severe] 1.880 0.001 True F6 status[T.severe] -1.729 0.006 True status[T.severe]_status[T.moderate] -1.654 0.023 True
Thus, F6 is not only depleted in severe samples relative to mild samples, but also relative to moderate samples.
Dunnett’s test
Dunnett’s test is useful when the scientific question specifically concerns each group versus a reference group (mild). It therefore omits the severe-vs-moderate comparison. A seed is supplied because the procedure uses bootstrapping.
>>> res_dunn = res.dunnett_test(seed=42) >>> res_dunn[res_dunn["Signif"]].round(3) Log(FC) SE W pvalue qvalue Signif FeatureID Comparison F2 status[T.moderate] 1.578 0.287 5.503 0.000 0.001 True status[T.severe] 1.880 0.300 6.272 0.000 0.000 True F6 status[T.severe] -1.729 0.343 -5.042 0.001 0.004 True
Trend test
Finally, the status column has a natural mild-to-moderate-to-severe ordering. The trend test evaluates ordered patterns in group effects. Here this ordering also matches the factor level order used by the fitted model.
>>> res_trend = res.trend_test(seed=42) >>> res_trend.round(3) W pvalue qvalue Signif FeatureID F1 0.513 0.28 1.00 False F2 1.880 0.00 0.00 True F3 0.031 1.00 1.00 False F4 0.833 0.13 0.78 False F5 0.485 0.32 1.00 False F6 1.729 0.00 0.00 True F7 0.232 0.79 1.00 False F8 0.539 0.23 1.00 False
The trend test again identifies F2 and F6, consistent with their increasing and decreasing abundance patterns across disease severity, respectively.
Structural zero test
The structural zero test supplied by the standalone function
struc_zeroidentifies features that are systematically absent from certain sample groups. As we may notice, F7 is absent from every severe sample but remains present in the mild and moderate groups. This information can be captured by the structural zero test:>>> struc_zero(table, metadata, "status").loc["F7"] mild False moderate False severe True Name: F7, dtype: bool
Refer to the documentation of
ancombcon updating test results with identified structural zeros.Pseudocount sensitivity analysis
The choice of pseudocount can affect differential abundance results. A sensitivity analysis repeats the complete analysis using multiple pseudocounts. A result is considered “pass” when its significance decision is identical across all analyses, and is marked as “robust” when it is consistently significant.
Note
This approach corresponds to the procedure used by ANCOMBC 2.10.1 and later. Earlier versions added alternative pseudocounts after bias correction. Since 2.10.1, each pseudocount is added before bias correction and the full analysis is repeated.
In addition to the default analysis (with the default
pseudocount=0), run the analysis on three pseudocounts: 0.1, 0.5 and 1.0, following the ANCOMBC package:>>> fits = [ancombc2(table, metadata, formula='status + age', grouping='status', ... pseudocount=p) for p in (0.1, 0.5, 1)] >>> signif = pd.concat([x.result["Signif"] for x in fits], axis=1)
Compare with the default analysis and append the output to the primary results:
>>> res_main["Pass"] = signif.eq(res_main["Signif"], axis=0).all(axis=1) >>> res_main["Robust"] = res_main["Signif"] & res_main["Pass"]
The result reveals that none of the three initially significant feature-covariate associations is robust to the tested pseudocounts. Thus, they should be interpreted cautiously rather than as supported conclusions.
>>> res_main.query("Robust == True").shape[0] 0
The same approach can be applied to global test, pairwise test and Dunnett’s test.
Note
For trend test, the ANCOMBC package uses the global test sensitivity decision rather than repeating the stochastic trend test at each pseudocount.
More broadly, the sensitivity analysis can also be used with
ancombcand other analyses involving pseudocounts.Sensitivity analysis can substantially reduce the number of significant results and thereby lower the risk of pseudocount-driven false positives, but this conservatism may reduce power and yield fewer significant results. The ANCOM-BC2 authors strongly recommend incorporating sensitivity analysis into the final assessment of taxon significance, unless maximizing power is the primary goal.
Variance regularization
Another mechanism for reducing the risk of false positives is the
var_quantileparameter (corresponding to the R package’ss0_percparameter), which adds the selected quantile (default: 0.05) of coefficient variances to the variance before inference. This helps reducing spurious significance caused by extremely small standard errors, especially for rare features. Larger values generally make the analysis more conservative and can further reduce false positives, but may also reduce statistical power.