Question

comparing 2 models with confounders with MAST (scRNA-seq DEA)

0

Entering edit mode

ramonmassoni ▴ 10

@ramonmassoni-19642

Last seen 6 months ago

Spain

Hi,

I want to conduct a differential expression analysis with MAST between two groups. In the analysis, I want to correct for a continuous confounder in a manner similar to this article: https://science.sciencemag.org/content/sci/suppl/2019/05/15/364.6441.685.DC1/aav8130VelmeshevSM.pdf . They say that "To identify genes differentially expressed due to the disease effect, likelihood ratio test (LRT) was performed by comparing the model with and without the diagnosis factor". Analogously, I would constract the following two models:

zlmsignal <- zlm(~ iscancer + biasscore, sca) zlmbackground <- zlm(~ bias_score, sca)

Thus, comparing both models I could avoid finding genes that vary due to "bias_score", and only detect those that are important in cancer. Could I use a likelihood test for that? If so, how?

Thanks a lot in advance.

scRNA-seq MAST confounders • 1.6k views

ADD COMMENT • link updated 5.0 years ago by Andrew_McDavid ▴ 270 • written 5.0 years ago by ramonmassoni ▴ 10

score 1 · Answer 1 · 2020-01-08

1

Entering edit mode

Andrew_McDavid ▴ 270

@andrew_mcdavid-11488

Last seen 11 weeks ago

United States

Edit: 10-Jan-2020 to fix syntax

You are 99% of the way there. You can compare zlmsignal to zlmbackground with

lrTest(zlmsignal, CoefficientHypothesis('is_richterTRUE'))

which just provides P values or

summary(zlmsignal, doLRT = 'is_richterTRUE')

which provides p values and coefficient estimates. The is_richterTRUE comes from how is_richter was named in the model matrix after being expanded as a contrast. If unsure, you can see how these were named with colnames(coef(zlmsignal, 'D')), or the error message also provides details. See section 4 of the MAIT vignette for some more examples.

ADD COMMENT • link 5.0 years ago Andrew_McDavid ▴ 270

0

Entering edit mode

Thanks a lot for your help!

So, to be more specific, this is how I'm building the SCA object and creating the models:

counts(richter_sce) <- as.matrix(counts(richter_sce))
logcounts(richter_sce) <- as.matrix(logcounts(richter_sce))
richter_sca <- SceToSingleCellAssay(richter_sce)
zlm_signal <- zlm(~ is_richter + cold_shock_score1, richter_sca)
zlm_background <- zlm(~ cold_shock_score1, richter_sca)

Then, I try both of your options, but I get the following errors:

lr_test <- lrTest(zlm_signal, zlm_background)
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘lrTest’ for signature ‘"ZlmFit", "ZlmFit"’

and

lr_test_summary <- summary(zlm_signal, doLRT = "is_richter")
Combining coefficients and standard errors
Calculating log-fold changes
Calculating likelihood ratio tests
Error in generateHypothesis(hypothesis, colnames(object@coefD)) : 
  Term(s) 'is_richter ,' not found.
Terms available: (Intercept) , is_richterTRUE , cold_shock_score1 , 
In addition: Warning messages:
1: In melt(coefAndCI, as.is = TRUE) :
  The melt generic in data.table has been passed a array and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(coefAndCI). In the next version, this warning will become an error.
2: In melt(lfc) :
  The melt generic in data.table has been passed a list and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the next version, this warning will become an error.

Which is weird since is_richter is definetely in the object:

table(richter_sca@colData$is_richter)
FALSE  TRUE 
 6254  1900

Do you have an idea of what might be wrong?

Thanks a lot in advance

R

ADD REPLY • link 5.0 years ago ramonmassoni ▴ 10

0

Entering edit mode

Oops, I made a typo in the syntax. Please see my updated answer above.

ADD REPLY • link 5.0 years ago Andrew_McDavid ▴ 270