MetaDeconfound()

The main (metadeconfoundR::MetaDeconfound()) analysis is a two step process:

Figure 1: metadeconfoundR pipeline overview.

First, significant associations between single omics features (like gut microbial OTUs) and metadata (like disease status, drug administration, BMI) are identified. Based on the data type of the respective metadata, either wilcox.test() (for binary), cor.test() (for continuous numerical) or kruskal.test() (for neither numerical nor binary) is used. All three tests are rank-based to minimize assumptions about data distribution (Fig. 1, left). In addition to collecting p-values for all computed tests, effect size is measured as Cliff’s Delta and Spearman’s Rho for binary and continuous data, respectively. Since there is no suitable effect size metric for categorical data with more than 2 levels, no value is reported here. It is recommended to introduce binary pseudo-variables for each level of the categorical metadata to partially circumvent this drawback.

In the second step, all hits are checked for confounding effects and a status is reported for each feature-metadata combination (Fig. 1, center and right; Fig2). A “hit” here is defined as a feature-metadata association with small enough fdr-corrected p-value and big enough effect size reported from the first, naive part of the pipeline. Thresholds for both parameters can be set via QCutoff and DCutoff when starting the analysis. Since confounding of signal can only happen with more than one metadata variable associated to a certain feature, all features with only one significant metadata association are trivially deconfounded and get status “No Covariates (OK_nc)”.

The actual confounder detection is done by performing a set of two likelihood ratio tests (LRTs) of nested linear models. For each possible combination of a feature and two of its associated metavariables, three models are fitted to the rank-transformed feature:

• lm(rank(feature) ~ covariate1 + covariate2), the full model
• lm(rank(feature) ~ covariate1), a model with only covariate1 as independent variable
• lm(rank(feature) ~ covariate2), a model with only covariate2 as independent variable.

LRTs reveal whether inclusion of covariate1 and/or covariate2 significantly improves the performance of the model. Random and/or fixed effects can be added to all models by the user. These additional effects will, however not be considered in the first naive association testing step of the pipeline.

Importantly, metadeconfoundR will always rank-transform the features during analysis.

Figure 2: detailed status labelling decision tree.

MetaDeconfound()

data(reduced_feature)
data(metaMatMetformin)

# check correct ordering
all(rownames(metaMatMetformin) == rownames(reduced_feature))
#> [1] TRUE
all(order(rownames(metaMatMetformin)) == order(rownames(reduced_feature)))
#> [1] TRUE

example_output <- MetaDeconfound(
  featureMat = reduced_feature,
  metaMat = metaMatMetformin,
  returnLong = TRUE,
  logLevel = "ERROR"
)

RandDataset_output <- MetaDeconfound(
  featureMat = reduced_feature,
  metaMat = metaMatMetformin,
  randomVar = c("Dataset"),
  returnLong = TRUE,
  logLevel = "ERROR"
)
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')

BuildHeatmap()

Minimal input consists only of an output object from the main MetaDeconfound() function. This will return in a ggplot2 heatmap with effect size as tile color and black asterisks or grey circles indicating significant and not confounded or confounded associations based on corrected p-values. The plot is clustered on both axes and features as well as meta-variables without any associations passing effect size and significance cutoffs (default: q_cutoff = 0.1, d_cutoff = 0.01) are removed (Fig. 3).

left <- BuildHeatmap(example_output)
right <- BuildHeatmap(RandDataset_output)
grid.arrange(left, right, ncol = 2)

Figure 3: default output of the BuildHeatmap() function

For both this default heatmap, as well as the alternative cuneiform plot (cuneiform = TRUE), a range of customizations are available. In Fig. 4 meta-variables not passing the effect size and significance cutoffs are manually kept in the plot (keepMeta), and the shown range of effect sizes is set from -1 to +1 (d_range = "full"). For a full list of options, again, refer to the help page.

BuildHeatmap(
  example_output,
  cuneiform = TRUE,
  keepMeta = colnames(metaMatMetformin),
  d_range = "full"
)

Figure 4: alternative cuneiform output of the BuildHeatmap() function

The BuildHeatmap() function returns a ggplot2 object. This makes it possible to perform some easy alterations manually (Fig. 5).

BuildHeatmap(example_output) +
  theme(
    legend.position = "none",
    axis.text.y = element_text(face = "italic"),
    plot.title = element_blank(),
    plot.subtitle = element_blank()
  )

Figure 5: post-plotting ggplot2 alterations

importLongPrior()

The newly added function ImportLongPrior() of the metadeconfoundR package allows for easy integration of a list of feature-metadata pairings into a new run of the main Metadeconfound() function.

These pairings must be supplied in the form of a long-format data.frame with the first column called “feature” containing names of omics features and the second column called “metaVariable” containing names of the respective associated metadata variable. When only these two columns are used as input for ImportLongPrior(), all given pairings are assumed to be significant.

Additional columns called “Qs” and “status”, as produced by metadeconfoundR:

print(example_output[101:105, ])

can also be supplied. Whenever the “status” column is part of the input, only feature–metadata pairings having a status other than “NS” are assumed to be significant.

To assure only features and metadata present in the current dataset are imported, the two input data.frames for the Metadeconfound() function, featureMat and metaMat, also need to be supplied to ImportLongPrior(). It is therefore important to have consistent naming of variables between datasets.

minQValues <- ImportLongPrior(longPrior = example_output,
                                featureMat = reduced_feature,
                                metaMat = metaMatMetformin)

In any case, ImportLongPrior() will return a wide-format data.frame containing significance values for the supplied pairings. If a “Qs” column was part of the longPrior input, these multiple-testing corrected p-values will be used for the output. If no “Qs” was supplied, artificial “significant” values of -1 will be used to mark significant associations.

print(minQValues[1:5, 1:5])

This data.frame can now be used as the minQValues input parameter of MetaDeconfound().

example_output2 <- MetaDeconfound(featureMat = reduced_feature,
                                  metaMat = metaMatMetformin,
                                  minQValues = minQValues)

This way, all significant feature–metadata pairings in the minQValues data.frame will be treated as potential confounders in the MetaDeconfound().

Note that this example is only to demonstrate the general process of integrating prior knowledge into a MetaDeconfound() analysis. Using the output of a MetaDeconfound() run as minQValues input for a second run with the exact same features and metadata will not lead to any new insights since the set of QValues calculated by MetaDeconfound() and the set supplied using the minQValues parameter are identical in this case.