Question: How do you test a contrast in MAST?
gravatar for Shian Su
8 days 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 8 days ago by Andrew_McDavid80 • written 8 days ago by Shian Su10
gravatar for Andrew_McDavid
8 days ago by
Andrew_McDavid80 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 8 days ago • written 8 days ago by Andrew_McDavid80

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 7 days 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 7 days ago • written 7 days ago by Andrew_McDavid80
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: 271 users visited in the last hour