Guidance on cov_pca and EM usage + error running cov_flash

[lucia-ramirez]

Hello! Thank you very much for creating and maintaining this package.

I'm currently using mashR to combine differentially expressed genes across multiple cell types and conditions. Since the cell types come from the same patients, I’m including a correlation matrix and testing different approaches (both null tests and EM), as well as comparing data-driven matrices (cov_pca and cov_flash).

I have a few questions and an issue to report:

1) Choosing the number of PCs in cov_pca

Do you have a general rule to decide how many pcs to use? I have up to 44 conditions (columns) in my bhat / shat matrices, though I might not include all of them in a single analysis.

In your simulations, you had 5 conditions and used 5 PCs, does that mean n_conditions = n_PCs? Or is 5 PCs generally sufficient even for larger numbers of conditions?

2) EM method for correlated measurements

When running mash_estimate_corr_em, should one include the canonical matrices (as in the vignette), or is it preferable to use only data-driven matrices, or both?

3) Error when running cov_flash

I’ve been trying to follow the vignette for flash but I’m encountering the following error/output:

cov_flash(D.simple$strong, factors="nonneg", tag="non_neg")
Adding factor 1 to flash object...
Adding factor 2 to flash object...
Adding factor 3 to flash object...
Adding factor 4 to flash object...
Adding factor 5 to flash object...
Adding factor 6 to flash object...
Adding factor 7 to flash object...
Adding factor 8 to flash object...
Adding factor 9 to flash object...
Adding factor 10 to flash object...
Adding factor 11 to flash object...
Adding factor 12 to flash object...
Adding factor 13 to flash object...
Adding factor 14 to flash object...
Adding factor 15 to flash object...
Adding factor 16 to flash object...
Adding factor 17 to flash object...
Adding factor 18 to flash object...
Adding factor 19 to flash object...
Adding factor 20 to flash object...
Adding factor 21 to flash object...
Adding factor 22 to flash object...
Adding factor 23 to flash object...
Adding factor 24 to flash object...
Adding factor 25 to flash object...
Adding factor 26 to flash object...
Adding factor 27 to flash object...
Adding factor 28 to flash object...
Factor doesn't significantly increase objective and won't be added.
Wrapping up...
Done.
Backfitting 27 factors (tolerance: 1.61e-03)...
  Difference between iterations is within 1.0e+03...
  Difference between iterations is within 1.0e+02...
  Difference between iterations is within 1.0e+01...
Error in if (any(s == 0)) { : missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In simpute.als(x, J, thresh, lambda, maxit, trace.it, warm.start,  :
  Convergence not achieved by 100 iterations
2: In handle_standard_errors(x, s) :
  Nonpositive SEs have been replaced by small positive SEs.

Both my bhat and shat matrices have no missing values, and shat contains no negative entries, so I’m unsure what might be triggering this issue. The cov_pca function runs successfully with the same input, suggesting the issue may be specific to cov_flash.

I'm using R 4.5.0, mashr_0.2.79 and flashier_1.0.7.

My data takes some time to run with flash but I could try to find a smaller subset or share the full data if needed.

Thanks!

[pcarbo]

Hi @lucia-ramirez I would suggest consulting our more recent work on Poisson mash for some more ideas on performing multi-condition analyses with (single-cell) RNA-seq data. It turns out that mashr works reasonably well, but perhaps more relevant to you is that the papers gives practical guidance on how to perform multi-condition analyses of scRNA-seq data (either with mash or with Poisson mash). Also, there are some tricky questions that arise when analyzing differences in expression across multiple conditions (as opposed to the simpler case of control vs. treatment). So hopefully this paper will give you some practical ideas about how to proceed with your analysis.

Feel free to reply here if you have additional questions.

[stephens999]

@willwerscheid do you have any ideas on this flash error?

[willwerscheid]

Hi @lucia-ramirez . I see you are using the CRAN version of flashier. Please try updating flashier to the most recent GitHub version - we implemented a fix for at least some instances of the error you're seeing.

[pcarbo]

Thanks @willwerscheid !

[lucia-ramirez]

Hi @pcarbo and @willwerscheid,

Thank you very much for the replies and for sharing the new paper. I read it and I think it's great idea, though it’s a bit of a shame that the current implementation isn’t yet adapted for datasets with replicates. In our case, we have between 7 and 30+ samples per condition and work with disease cohorts, which can be quite heterogeneous. I’m not sure collapsing the replicates for the Poisson model would be a good idea for us.

That said, I think as more large-scale single-cell atlases and meta-analyses become available, there will definitely be growing interest in extending mash-poisson to handle replicates and would be a very valuable addition for the single-cell community!

Regarding the cov_flash issue, I was initially testing it with a small number of conditions (N = 4), just to get the workflow running. Now that I’ve rerun it with all 44 conditions, it seems to be working smoothly, so the earlier issue might have been related to the small test size.

Finally, I was wondering if you might be able to clarify the two questions from my original post?

1. Choosing the number of PCs in cov_pca
2. Using data-driven matrices as input for mash_estimate_corr_em

Many thanks again for your time and for your replies!

[pcarbo]

@lucia-ramirez That's great news! Partly my reasoning for pointing you to the Poisson mash paper was to mention that in that paper we used mash for a similar type of data to you (single-cell RNA-seq data with different conditions), and some of the details of that analysis may be useful to you. In particular, note that we used the "common baseline" variant of mash (see the supplement of that paper for explanations and technical details).

For Q2, yes you would want to include all the covariance matrices in the call to mash_estimate_corr_em.

For Q1, essentially, you want a large number of covariance matrices — as many as you can handle. The main downside to having "too many" covariance matrices is that it slows down the calculations. A priori, 100 seems like a reasonable number, but you may want to try a few different numbers (e.g., 50, 100, 200) to see how it impacts your mash results.

Your questions about estimating the covariances also brings to mind that we have a newer R package udr (and a paper on the new methods) that should improve the estimates of the covariance matrices, especially considering these are high-dimensional (44 x 44) matrices. The key idea was to "regularize" the estimates. (The new methods are also more robust to mis-specifying the number of components, which will help with Q1.)

[lucia-ramirez]

Thank you very much for the detailed reply, that’s really helpful!

From a quick look at the udr paper, it does seem like a good idea to use it to improve the estimation of the covariance matrices, especially in our high-dimensional setting. Are there any plans to release a vignette or workflow showing how to integrate mashR and udr?

I had a look at the application and code for the GTEx cis-eQTLs example to see if we could adapt that for our project, but I am not entirely sure how to extract the posterior summaries when combining the two packages. Any guidance or pointers would be greatly appreciated.

[pcarbo]

@lucia-ramirez We do not have a vignette yet — this is on our to-do list. (We also would like to better integrate udr and mashr.) I suppose one complication of using udr might be choosing the penalty strength (lambda), which previously we did by cross-validation (CV). So in short we need to improve the analysis interface. :)

If you are motivated, the key steps might look something like this:

library(mashr)
library(udr)
lambda <- 1 # Arbitrary; ideally choose via CV, but 1 is better than no penalty (0).
mash_data <- mash_set_data(Z,V = V)
U_canonical <- cov_canonical(mash_data)
fit0 <- ud_init(Z,V = V,U_scaled = U_canonical,n_rank1 = 0)
fit_ed <- ud_fit(fit0,verbose = TRUE,
                 control = list(unconstrained.update = "ed",
                                lambda = lambda,penalty.type = "iw",
                                maxiter = 1000))
U_ed      <- lapply(fit_ed$U,"[[","mat")
fit_mash  <- mash(mash_data,U_ed)

You could compare the output from udr::ud_fit() to what you get using mashr::cov_ed().