comparing 2 models with confounders with MAST (scRNA-seq DEA)
1
0
Entering edit mode
ramonmassoni ▴ 10
@ramonmassoni-19642
Last seen 12 months ago

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 • 273 views
ADD COMMENT
1
Entering edit mode
@andrew_mcdavid-11488
Last seen 4 months ago

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

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

ADD REPLY

Login before adding your answer.

Traffic: 335 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