Error in chol.default

[grenkoca]

Hello,

I've been trying to run Lamian for the last week on a single-cell data set and have been running into some trouble. I have a normalized log gene-expression matrix with about 32,000 cells with 17,000 genes. Each cell has a precomputed pseudotime in the range [0,1]. I am trying to run both Lamian's TDE and XDE tests on this data.

When I initially ran either test, I was trying to multithread the job. It got the point where it would run the test and print fitting model: overall: CovariateTest (Model 3 vs.1), or ConstantTest (Model 1) ..., and then throw the following error:

Error: $ operator is invalid for atomic vectors
In addition: Warning messages:
1: In mclapply(seq_len((permuiter + 1)), function(i) { :
  scheduled cores 1, 2, 3, 5, 6, 8, 9, 10, 12, 15, 17, 20, 35, 36, 37, 39, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64 did not deliver results, all values of the jobs will be affected
2: In mclapply(seq_len((permuiter + 1)), function(i) { :
  scheduled cores 40, 33, 25, 34, 4, 21, 13, 26, 22, 28, 31, 14, 27, 24, 29, 18, 32, 7, 23, 19, 30, 11, 38, 16 encountered errors in user code, all values of the jobs will be affected
Execution halted

To get around using mcapply I reduced ncores to 1, and then ran it again. Now it is giving the following error after getting to the same point:

Error in chol.default(matrix(omega[g, ], nrow = nb)) : 
  the leading minor of order 2 is not positive definite
Calls: lamian.test ... lapply -> FUN -> <Anonymous> -> <Anonymous> -> chol.default
Execution halted

What exactly is going wrong here? When I subsample the data and push it through the same pipeline (see attached) it works perfectly. run_lamian_XDE.R.txt

I would be happy to share some data with you to help troubleshoot.

Some additional information that might be helpful: R version = 3.6.2 OS = Scientific Linux 6 RAM = 256Gb

Thanks, Caleb

Tagging @letaylor

[Winnie09]

Hello Caleb,

Thanks for the interest in Lamian! Could you apply str to these objects which pass to the function Lamian.test() and show the results here? design_mat, adata_pseudotime_u, adata_cellanno, adata_expr_t

For example, str(design_mat).

Thank you!

Best regards, Wenpin

[grenkoca]

Hi Wenpin,

Thanks for your quick response. Absolutely, please see here:

adata_expr_t:

 num [1:17523, 1:21383] 0.991 0.353 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:17523] "SAMD11" "NOC2L" "KLHL17" "PLEKHN1" ...
  ..$ : chr [1:21383] "AGAAATGAGGCTATCT-1-E12-01" "ATCGGATCAGTCTCTC-1-E12-01" "GATGCTAAGAGAGGGC-1-E12-01" "GGAGGTATCGGTCTGG-1-E12-01" ...

adata_cellanno:

Classes ‘data.table’ and 'data.frame':  21383 obs. of  2 variables:
 $ Cell         : chr  "AGAAATGAGGCTATCT-1-E12-01" "ATCGGATCAGTCTCTC-1-E12-01" "GATGCTAAGAGAGGGC-1-E12-01" "GGAGGTATCGGTCTGG-1-E12-01" ...
 $ Experiment_ID: chr  "E12-01" "E12-01" "E12-01" "E12-01" ...
 - attr(*, ".internal.selfref")=<externalptr> 

adata_pseudotime_u:

Named num [1:21383] 0 0.00213 0.0214 0.00605 0.00649 ...
 - attr(*, "names")= chr [1:21383] "AGAAATGAGGCTATCT-1-E12-01" "ATCGGATCAGTCTCTC-1-E12-01" "GATGCTAAGAGAGGGC-1-E12-01" "GGAGGTATCGGTCTGG-1-E12-01" ...

design_mat:

 num [1:42, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:42] "E12-01" "E12-02" "E12-03" "E12-04" ...
  ..$ : chr [1:2] "intercept" "group"

Thanks, Caleb

[Winnie09]

Thank you, Caleb! It seems that the input data type is what Lamian needs. Could you try these:

1. for the object adata_expr_t: filter out genes with 0 expression in too many cells. For example, filtering by adata_expr_t <- adata_expr_t[rowMeans(adata_expr_t > 0) > 0.01,]. Here, 0.01 (genes should have expression in at least 1% of the cells) can be replaced by 0.001 or other numbers you perfer.
2. make sure there is no cell with 0 expressions across all genes.
3. for the object adata_pseudotime_u: use the current pseudotime times 1000, and make sure the minimum is 1.
4. make sure there is enough memory space. In troubleshooting, setting ncores=1 is a good strategy to rule of memory issues.

If the above doesn't work, please try to reorganize the pseudotime so that it is a vector of consecutive integers, i.e., 1,2,3,... since this is the typical output from TSCAN. Please let me know if this could work (if yes, I should double-check the codes). Thank you!

Wenpin

[grenkoca]

Hi Wenpin,

I tried all your suggestions and it still failed. I set ncores=1 for debugging and it failed at the same stage with the same error.

Here's the new structure of adata_pseudotime_u:

 Named int [1:21383] 1 737 734 732 728 726 150 151 721 715 ...
 - attr(*, "names")= chr [1:21383] "AGAAATGAGGCTATCT-1-E12-01" "ATCGGATCAGTCTCTC-1-E12-01" "GATGCTAAGAGAGGGC-1-E12-01" "GGAGGTATCGGTCTGG-1-E12-01" ...

This data was generated in python by replacing the float pseudotime value with np.argsort(pseudotimes) which gave them sequential values.

Here is the error it throws at the same point as before:

Error in chol.default(crossprod(phiss)) : 
  the leading minor of order 2 is not positive definite
Calls: lamian.test ... lapply -> FUN -> <Anonymous> -> <Anonymous> -> chol.default
Execution halted

Please see here for the script I'm using.

Thank you! Caleb

[Winnie09]

Hi Caleb,

The computation cluster I am using is still down today. I will work on it when the cluster is back to normal. Thanks for the patience!

Best, Wenpin

[grenkoca]

Hi Wenpin,

No worries. Let me know how I can help resolve this. If you want I'd be happy to share our data privately to recreate the error.

Thanks, Caleb

[jr-leary7]

Hi, what was the resolution of this issue ? I'm having the same problem today, & have made sure my inputs match the required format (filtered counts matrix, integer pseudotime, etc.). Thanks !

[Sayyam-Shah]

Hello, I am also getting the same error @Winnie09. I would love your support in debugging this.