Question: How do you test a contrast in MAST?
gravatar for Shian Su
15 months ago by
Shian Su10
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Shian Su10 wrote:

I am trying to analyse some single cell RNA seq data with MAST, currently I'm using an edgeR pipeline and I cannot work out how to properly test a contrast in MAST.

The experimental setup is a 3x2 with 3 strains and 2 conditions, I'd like to test differential expression between the conditions within each strain. I have set up a "strain-condition" style factor for my samples and I'm using a `~0 + group` design. My understanding was that I could use

    zlm_output <- zlm(~0 + group, x)

    summary_contr <- summary(zlm_output, doLRT = c(1, -1, 0, 0, 0, 0))

To test my contrast of interest. However this gives me the following error

    Error in logFC(zlmfit, contrast0, contrast1) : 
      Assuming comparision to intercept, but I can't figure out what coefficient that corresponds to.  Provide `contrast0`.

Could someone guide me in the correct usage of this package for testing contrasts?


ADD COMMENTlink modified 15 months ago by Andrew_McDavid150 • written 15 months ago by Shian Su10
gravatar for Andrew_McDavid
15 months ago by
Andrew_McDavid150 wrote:
  1. You can fit such a model a "cellmean model" (zlm_cellmean = zlm(~ 0 + group, x)) though it is possibly a little less statistically efficient [1] than using contrasts that explicitly includes an intercept:

    ## assuming you manually crossed strain and condition to get group
    zlm_dummy = zlm(~strain*condition, x)
    ## Not so friendly
    contrasts(colData(x)$strain) = "contr.sum"
    contrasts(colData(x)$condition) = "contr.sum"
    ## classical anova -- maybe this is what you want?
    zlm_sum_to_zero = zlm(~strain*condition, x)
  2. The summary method (?'summary,ZlmFit-method') supports contrasts offered as a character vector, which seemingly doesn't cover the zlm_cellmean case.

  3. To test contrasts symbolically or numerically you may offer them to lrTest
    lr_anova <- lrTest(zlm_cellmean, Hypothesis('groupLevel2 - groupLevel1'))
    # columns give a set of contrasts. more than one column means a joint test.
    lr_anova2 <- lrTest(zlm_cellmean, as.matrix(c(1, -1, 0, 0, 0, 0)))
    all.equal(lr_anova, lr_anova2)

​[1] This is due to a heuristic in the coefficient shrinkage that shrinks "intercept" columns less than "treatment" columns.

ADD COMMENTlink modified 15 months ago • written 15 months ago by Andrew_McDavid150

Thanks, I guess having an intercept model is fine. Is this interface for contrasts documented anywhere? I couldn't quite work it out from the vignette and help pages.

ADD REPLYlink written 15 months ago by Shian Su10

Can you be more specific when you say "interface for contrasts" The contrasts(colData(x)$strain) = "contr.sum" business is base R and determines the coding for your factors, hence interpretation of your coefficients.  Venables & Ripley is the canonical reference for that.  

Specifying linear combinations of coefficients to test is documented in ?"ZlmFit-class".

ADD REPLYlink modified 15 months ago • written 15 months ago by Andrew_McDavid150
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 210 users visited in the last hour