Model covariates

[rdrummond]

I would like to consider sample covariates (sequencing target and platform used, PCR protocol, etc.) when estimating taxon biases to normalize metagenomics data.
I am working with mock samples of known compositions, but presence/absence patterns vary across samples. My goal is to calibrate a model to later normalize cancer samples of unknown composition, in which I am searching microbiota patterns.

Could you please advise me on how to add covariates to the Metacal model?
Thanks in advance.

[mikemc]

Hi @rdrummond, this is a good question. You have a few different options right now, and which you will prefer will likely come down to how comfortable you are with different statistical frameworks.

For some intuition about the various solutions you have available, it is helpful to first talk about the case where all taxa are in all samples. Then you can simply use a log-ratio transformation (ALR, CLR, whatever), and perform multivariate linear regression with the built in lm() function, using R's formula interface to supply a formula such as observed ~ 1 + pcr_protocol + offset(actual), where observed is the log-ratio transformed counts, and actual is the matrix of log-ratio-transformed nominal ("actual") abundances. In this case, the coefficients for the intercept term will be the bias in the first protocol, and the coefficients of the pcr_protocol terms will be the differential bias of the other protocols relative to the first. (The confidence intervals etc from lm() probably won't be very reliable, but bootstrapping or various alternatives to lm() could help deal with this.)

Unfortunately, this only works when all taxa are in all samples. For the case where presence/absence varies, some options are:

1. A simple option for discrete covariates is to split your data into subsets corresponding to different covariate values, separately estimate bias on each subset, and compare the estimates. For example, suppose you have two PCR protocols. You can estimate bias on the data from Protocol 1 and the data from Protocol 2 separately, and compare the values. And if you compositionally subtract the estimates (by element-wise division), that tells you the effect of using one protocol versus another. This is essentially what we did in our 2019 paper to compare the bias in samples started as cell mixtures and samples that stated as DNA mixtures, and compute the implied impact of DNA extraction.

2. A more general and (for some) natural approach would be to perform linear regression model on the log-ratio-transformed abundances like described with lm() above. I have developed an extension of the metacal estimation method that can do this, using bootstrapping for uncertainty estimation, but have not released it yet mainly because I need to do some refactoring and testing for the case where not all taxa are in all samples. I will try to get this out but it might take a couple weeks.

3. Ni Zhao and Glen Satten have made an R package https://github.com/zhaoni153/MicroBias for doing linear regression on log proportions in such a way as to be consistent with the metacal model and allows having different taxa in different samples. The model shown in the Readme right now is actually not compatible with metacal but I think this is due to a typo (The model in the Readme may be misspecified zhaoni153/MicroBias#2). I don't think their manuscript is publically available yet but you might email them a request, and the Readme and function doc's might be enough to get started. (I haven't tried it yet).

4. Another option I recommend considering if you are comfortable with Bayesian analysis is @jsilve24's fido package, and specifically the Multinomial Logistic Normal models that can be fit with the pibble() function. With this approach, you must take all taxa to be in all samples, but they can be present at very small frequencies (so small that you are very likely to observe zero reads). Provided that you set your model matrix and priors correctly, you can get coefficient estimates that will be compatible with what you would get from the metacal approach. Supplying the nominal ("true") abundances for the mock communities is a bit tricky and I've experimented with two approaches that both seem to work. I will try to post a bare-bones example next week to show how to do this.

I'm hoping to post some example workflows doing this type of analysis in the coming 1-2 months that will demonstrate the basic mechanics and perhaps more importantly, how to interpret the coefficients, which is in fact quite tricky due to the compositional nature of the bias parameters.

[mikemc]

Note, my reply above addresses the question of how to estimate coefficients that describe how a sample covariate (such as a protocol choice, or an experimental batch) changes the set of relative efficiencies associated with a measurement. However, doing so is unnecessary if you wish only to calibrate samples that are measured by a protocol that you also used to measure mocks. In that case, you can just estimate bias for that protocol on the subset of mock samples in which that protocol was applied (suggestion 1. above), and use that bias estimate for your calibration.

[mikemc]

I've finally published a demonstration of using fido for bias inference from mocks here. I give some idea of how to include covariates at the end, but still need to flesh this part out.