CSS scaling factors in fitZig model matrix

[mcalgaro93]

Hello, I'm using metagenomeSeq and many other packages for comparisons. My issue is about the way model matrix for method fitZig is obtained. The code says:

# Initializing the model matrix
    if (useCSSoffset == TRUE){
        if (any(is.na(normFactors(obj)))) {
          stop("Calculate the normalization factors first!")
        }
        mmCount <- cbind(mod, log2(normFactors(obj)/1000 + 1))
        colnames(mmCount)[ncol(mmCount)] <- "scalingFactor"
    } else { 
        mmCount <- mod
    }

So by default model matrix is made by intercept, coefficients if they are present in mod and log2(normFactors(obj)/1000 + 1). But in the paper is also written that:

So I understand the 1000 at denominator. For my data I prefer to use the median of the scaling factor as denominator (as suggested above). So I put useCSSoffset = FALSE and I manually calculate offsets... The thing I don't get is the fact that in supplementary notes is reported a different log2 of the scaling factor.

So which is the correct way to include the scaling factor in the model?

• log2(normFactors(obj)/median(normFactors(obj)) + 1);
• log2((normFactors(obj)+1)/median(normFactors(obj));

Wouldn't be better to put one of the solutions above as default instead of scaling per 1000? Thanks in advance, best regards.

Matteo

[hcorrada]

Thank you for writing. The equation on the supplementary material is incorrect, log2((s/N) + 1) is correct, as in the code you link.

We've vacilated between the 1000 vs. median scale factor in our analyses and resort to setting the offsets manually for the latter. With N=1000 interpretation is "count per thousands", which median scale factor doesn't have.