Hi there,

I'm exploring ballgown as a method for differential gene expression, but running into what I believe should be a simple problem in specifying models for testing for differential expression.  The function stattest performs this test and the documentation reports that arguments mod and mod0 can be used to specify custom models with and without desired covariates respectively when testing for differential expression. From the manual:

"mod and mod0 are optional arguments. If mod is specified, you must also specify mod0. If neither is specified, mod0 defaults to the design matrix for a model including only a library-size adjustment, and mod defaults to the design matrix for a model including a library-size adjustment and covariate."

What I would like to do is to use the default mod0 including only the library-size adjustment and test it against a model with two independent covariates (sequentially or together), but I haven't been able to identify the code that implements the default library-size adjustment.  This should be trivial, but I haven't been able to solve it efficiently. Any perspective you can offer would be welcome!

Andrea

So here is the awesome part about the new github mirror. If you go here, and search for ballgown, then click the R directory, scroll down to stattest.R, and click that, you can see all the code with lots of useful comments. About halfway down is this:

 ## library size adjustment
    if(is.null(libadjust)){
        libadjust = apply(expr, 2, function(x){
            lognz = log2(x[x!=0] + 1)
            q3 = quantile(lognz, 0.75)
            sum(lognz[lognz<q3])
        })

Which specifies the default library size (where you restrict to counts for genes in the first three quartiles). And a bit further down, there is this bit:

if(!identical(libadjust, FALSE)){
                eval(parse(text=paste0("mod0 = model.matrix(~ libadjust",
                    variable_list, ")")))

Which is just a programmatic way of doing this at the R prompt:

mod0 <- model.matrix(~ libadjust)

Does that help?

Hi James,

Thanks so much that definitely does help!

For other novices solving the same problem, you will need to define expr appropriately for your analysis, then run the command to generate the function libadjust before running the mod0 model matrix analysis.

Cheers,

Andrea

Hey Andrea,

Sorry for the lack of clarity in the documentation here! The default library size adjustment is (rather loosely) defined in the documentation for the "libadjust" parameter on the stattest help page, but I'll work on making that more obvious over the next few days.

I also wanted to clarify that mod and mod0 need to be nested models, i.e., mod needs to contain all covariates that are included in mod0. So in your example, make sure to include library size in mod (in addition to mod0). The function will give you an error if you try to compare non-nested models, but I'll also try to document the requirement more explicitly. 

Hope this helps!