Corrected protein data matrix

[Arthfael]

This is not an issue, more like a question/suggestion:
Suppose I have run msqrob2 successfully. I have my significant proteins and am a happy little researcher. Good. But now, I also want for whichever reason to perform further inference down the road with another algorithm (SAINTexpress for instance, which like MSqRob is quite well suited to missing values, but doesn't have formula-based experiment structure modeling), or to plot a heatmap, write per samples expression data in a table for end users. For all of these I would need a proteins x samples expression matrix. Ideally, I would like this to reflect the modeling which MSqRob has based its inference on, and which would make for more robust estimates and remove some undesirable effects.
My question would be: if I wanted to use MSqRob as a modeling/data correction step prior to these further steps, could I currently get a predicted per samples quantitative matrix using the models? I know it would not have all of the additional modeling information build into MSqRob, but it would still be useful for these purposes - unless I am mistaken!
I know how to get it from the QFeatures object, but could verify that after I run msqrob() on it the values are not updated. I also explored the models created in the object but they do not appear to have any per samples prediction.

So my questions are:

• is this currently already feasible?
• ... if not, is it because it would be an absolute no go, an absurdity?...
• ... or is it something which could be done (with the obvious caveats), and if so would you consider adding this in the future?

[cvanderaa]

Hello @Arthfael,

Below you can find some code that I use for getting corrected data. I'm sorry, it is implemented very quickly, so there may be bugs and I haven't documented the code. The main function is msqrobRemoveBatchEffect() which takes the same arguments as msqrob(), except batchVars which is a vector of characters providing the names of the variables in your formula for which you wish to correct for (e.g. formula = ~ bioFactorA + bioFactorB + techFactorC + techFactorD would take batchVars = c("techFactorC" + "techFactorD")). The function takes an object of class SummarizedExperiment (you can easily write a wrapper yourself for QFeatures object) and returns an SummarizedExperiment object where the assay contains the batch corrected quants.

⚠️ Disclaimer ⚠️ this is super experimental and has limitations because the implementation:

• Cannot properly handle complex designs that involve interactions, or random effects other than a random slope.
• Cannot handle models with feature-level variables, that is models fit using msqrobAggregate().

I cannot see a way to implement an automated batch correction approach that is sufficiently general to account for these cases. I therefore don't think this will/should be integrated into msqob2.

I hope however this already fulfills your needs.

The code

.parseLMMFormula <- function(object, f) {
    cd <- data.frame(colData(object))
    cd$..lhs.. <- runif(nrow(cd))
    rhs <- as.character(f)[[2]]
    lme4::lFormula(formula(paste("..lhs.. ~", rhs)), data = cd)
}

.getX <- function(lf, ridge) {
    X <- lf$X
    if (ridge) {
        notIntercept <- colnames(X) != "(Intercept)"
        newParamNames <- paste0("ridge", colnames(X)[notIntercept])
        colnames(X)[notIntercept] <- newParamNames
    }
    return(X)
}

.getZ <- function(lf) {
    re <- lf$reTrms
    Z <- t(re$Zt)
    if (!all(unlist(re$cnms) == "(Intercept)")) 
        stop("Batch correction can only handle random intercepts.")
    prefix <- paste0(re$cnms[re$Lind], names(re$cnms)[re$Lind])
    colnames(Z) <- paste0(prefix, colnames(Z))
    return(Z)
}

.getCoefs <- function(models, paramNames) {
    coefs <- matrix(
        NA, ncol = length(paramNames), nrow = length(models),
        dimnames = list(
            names(models),
            paramNames
        )   
    )
    for (i in rownames(coefs)) {
        coefsi <- models[[i]]@params$coefficients
        coefs[i, ] <- unname(coefsi[colnames(coefs)])
    }
    return(coefs)
}

.predictBatchEffect <- function(X, beta, batchVars) {
    if (any(grepl(":", colnames(X))))
        warning("Batch correction with interactions is experimental (hence unreliable).")
    batchPattern <- paste0("(^", batchVars, ")")
    batchPattern <- paste0(batchVars, collapse = "|")
    batchParams <- grep(batchPattern, colnames(beta))
    t(tcrossprod(X[, batchParams], beta[, batchParams]))
}

msqrobRemoveBatchEffect <- function(
    object, formula, ridge, batchVars, modelName = "msqrobModels"
) {
    models <- rowData(object)[[modelName]]
    lf <- .parseLMMFormula(object, formula)
    XZ <- cbind(.getX(lf, ridge), .getZ(lf))
    betau <- .getCoefs(models, colnames(XZ))
    yHatBc <- .predictBatchEffect(XZ, betau, batchVars)
    yBc <- assay(object) - yHatBc[rownames(object), colnames(object)]
    rowData(object)[[modelName]] <- NULL
    SummarizedExperiment(
        assays = list(batchCorrected = as.matrix(yBc)), 
        colData = colData(object),
        rowData = rowData(object)
    )
}