Question: How do you test a contrast in MAST?
gravatar for Shian Su
2.1 years ago by
Shian Su20
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Shian Su20 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 2.1 years ago by Andrew_McDavid190 • written 2.1 years ago by Shian Su20
Answer: How do you test a contrast in MAST?
gravatar for Andrew_McDavid
2.1 years ago by
Andrew_McDavid190 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 2.1 years ago • written 2.1 years ago by Andrew_McDavid190

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 2.1 years ago by Shian Su20

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 2.1 years ago • written 2.1 years ago by Andrew_McDavid190
Please log in to add an answer.


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