How do you test a contrast in MAST?
1
0
Entering edit mode
Shian Su ▴ 40
@shian-su-9869
Last seen 2 days ago
Walter and Eliza Hall Institute of Medi…

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?

 

MAST single-cell differential gene expression • 2.9k views
ADD COMMENT
3
Entering edit mode
@andrew_mcdavid-11488
Last seen 6 weeks ago
United States
  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 COMMENT
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 835 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6