UCLouvain-CBIO / scp

Single cell proteomics data processing
https://uclouvain-cbio.github.io/scp/index.html
19 stars 2 forks source link

Differential analysis: non-conformable arguments #65

Closed samgregoire closed 4 days ago

samgregoire commented 3 weeks ago

I ran into this issue when using the scpDifferentialAnalysis() function.

> table(sce$CellType)

CD38_neg CD38_pos      LSC     prog 
      38       49       46       53 
> 
> (daRes <- 
+     scpDifferentialAnalysis(
+       sce, contrast = list(c("CellType", "CD38_neg", "CD38_pos"),
+                            c("CellType", "CD38_neg", "LSC"),
+                            c("CellType", "CD38_neg", "prog"),
+                            c("CellType", "CD38_pos", "LSC"),
+                            c("CellType", "CD38_pos", "prog"),
+                            c("CellType", "LSC", "prog"))))
Error in t(l) %*% mm : non-conformable arguments

This is what I get when I run the traceback.

> traceback()
10: .multiLevelsToContrastMatrix(contrast, levels)
9: .levelsToContrastMatrix(contrast, levs)
8: FUN(X[[i]], ...)
7: lapply(X = X, FUN = FUN, ...)
6: sapply(fits, function(fit) {
       coef <- scpModelFitCoefficients(fit)
       vcov <- scpModelFitVcov(fit)
       levs <- scpModelFitLevels(fit)[[contrast[[1]]]]
       if (length(levs) <= 1) 
           return(c(logFc = NA, se = NA))
       contrastMat <- .levelsToContrastMatrix(contrast, levs)
       sel <- grepl(contrast[[1]], names(coef))
       logFc <- contrastMat %*% coef[sel]
       se <- sqrt(contrastMat %*% vcov[sel, sel] %*% t(contrastMat))
       if (length(contrastMat) == 1) {
           logFc <- logFc * 2
           se <- se * 2
       }
       c(logFc = logFc, se = se)
   })
5: .contrastToEstimates(object, contrast, name)
4: FUN(X[[i]], ...)
3: lapply(contrasts, function(contrast) {
       est <- .contrastToEstimates(object, contrast, name)
       .computeTTest(est$logFc, est$se, df)
   })
2: .scpDifferentialAnalysisOnContrast(object, contrasts, name)
1: scpDifferentialAnalysis(sce, contrast = list(c("CellType", "CD38_neg", 
       "CD38_pos"), c("CellType", "CD38_neg", "LSC"), c("CellType", 
       "CD38_neg", "prog"), c("CellType", "CD38_pos", "LSC"), c("CellType", 
       "CD38_pos", "prog"), c("CellType", "LSC", "prog")))

I couldn't find a solution myself. Do you have any idea what's going on here? If it helps, the dataset I'm using is from scpdata::petrosius2023_AstralAML().

cvanderaa commented 3 weeks ago

This is indeed unexpected... Did you try to use the latest code from the pending PR (#63)? I found a few bugs while unit testing, so it may be solved there. If you still get an error with the latest code, could you send me your modelled sce object? That would facilitate debugging.

lgatto commented 3 weeks ago

PR #63 has now been merged. Use

BiocManager::install("UCLouvain-CBIO/scp")

to install it.

It has also been pushed to Bioc (devel) and should become available in a day or two.

samgregoire commented 3 weeks ago

I still get a similar (altought slightly different) error with the updated version of the code: Error in contrastMat %*% vcov[sel, sel] : non-conformable arguments

> traceback()
8: FUN(X[[i]], ...)
7: lapply(X = X, FUN = FUN, ...)
6: sapply(fits, function(fit) {
       coef <- scpModelFitCoefficients(fit)
       vcov <- scpModelFitVcov(fit)
       levs <- scpModelFitLevels(fit)[[contrast[[1]]]]
       if (length(levs) <= 1) 
           return(c(logFc = NA, se = NA))
       contrastMat <- .levelsToContrastMatrix(contrast, levs)
       sel <- grepl(contrast[[1]], names(coef))
       logFc <- contrastMat %*% coef[sel]
       se <- sqrt(contrastMat %*% vcov[sel, sel] %*% t(contrastMat))
       c(logFc = logFc, se = se)
   })
5: .contrastToEstimates(object, contrast, name)
4: FUN(X[[i]], ...)
3: lapply(contrasts, function(contrast) {
       est <- .contrastToEstimates(object, contrast, name)
       daTable <- .computeTTest(est$logFc, est$se, df, weights)
       metadata(daTable)$contrast <- contrast
       daTable[!is.na(daTable$Estimate), ]
   })
2: .scpDifferentialAnalysisOnContrast(object, contrasts, name)
1: scpDifferentialAnalysis(sce, contrast = list(c("CellType", "CD38_neg", 
       "CD38_pos"), c("CellType", "CD38_neg", "LSC"), c("CellType", 
       "CD38_neg", "prog"), c("CellType", "CD38_pos", "LSC"), c("CellType", 
       "CD38_pos", "prog"), c("CellType", "LSC", "prog")))

You can download my modelled sce object here.

cvanderaa commented 6 days ago

Hello Sam!

Sorry, I finally took the time to look into the issue, and this was indeed due to a nasty little bug :sweat_smile:. The solution is provided through this PR.

You can install it with

BiocManager::install("UCLouvain-CBIO/scp", ref = "issue65")

Could you confirm it works for you as well?

samgregoire commented 6 days ago

It does! Thank you