how long will "fitting gamma coefficients" step take?

[wulabupenn]

Hi Team velocyto,

I know the computational time depends on a lot of factors: CPU, memory and size of data set. Our HPC have 48 processors and 256 GB memory. I run the example listed in README : http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html. The emat is 736 common genes x 2600 cells, which takes only few seconds to finish the "gene.relative.velocity.estimates" step. Here I have a data set comprised by 944 common genes x 5000 cells. It has already run about 2 hours but still stay at the step "fitting gamma coefficients ...". I am wondering whether something wrong with my installation or it just need longer time to finish. Any suggestion will be helpful.

Thanks, Peng

[pkharchenko]

It shouldn’t take that long - the scaling is essentially linear with respect to the number of genes and something like n*log(n) for cells. I would first check to make sure it’s actually running (I’ve seen parallel processing hang once in a while). Are there any more detailed messages? Are you running with verbose=T?

-peter.

On Aug 22, 2018, at 20:26, wulabupenn notifications@github.com wrote:

Hi Team velocyto,

I know the computational time depends on a lot of factors: CPU, memory and size of data set. Our HPC have 48 processors and 256 GB memory. I run the example listed in README : http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html. The emat is 736 common genes x 2600 cells, which takes only few seconds to finish the "gene.relative.velocity.estimates" step. Here I have a data set comprised by 944 common genes x 5000 cells. It has already run about 2 hours but still stay at the step "fitting gamma coefficients ...". I am wondering whether something wrong with my installation or it just need longer time to finish. Any suggestion will be helpful.

Thanks, Peng

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[wulabupenn]

@pkharchenko Thank you for the prompt reply. Here is detail information about the data set:

`> str(emat1) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:8088910] 0 1 2 3 4 6 7 8 9 11 ... ..@ p : int [1:5722] 0 2242 4576 6716 9034 11191 13480 15724 17943 20229 ... ..@ Dim : int [1:2] 3139 5721 ..@ Dimnames:List of 2 .. ..$ : chr [1:3139] "AASDHPPT" "AB019441.29" "ABBA01017803.1" "ABCE1" ... .. ..$ : chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_ATTTTTCTTTTN" "H7-CM-batch7-d0_CACTTACATACT" ... ..@ x : num [1:8088910] 2 19 1 5 1 5 1 19 1 3 ... ..@ factors : list()

str(nmat1) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:5032887] 0 1 5 6 9 14 16 18 20 22 ... ..@ p : int [1:5722] 0 2020 3865 5920 7976 10072 11793 13979 16081 18112 ... ..@ Dim : int [1:2] 5464 5721 ..@ Dimnames:List of 2 .. ..$ : chr [1:5464] "A2ML1-AS1" "A2MP1" "AACS" "AAK1" ... .. ..$ : chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_ATTTTTCTTTTN" "H7-CM-batch7-d0_CACTTACATACT" ... ..@ x : num [1:5032887] 1 3 1 1 1 1 1 2 1 1 ... ..@ factors : list() str(cell.dist) Class 'dist' atomic [1:16362060] 0.112 0.225 0.171 0.192 0.123 ... ..- attr(, "Labels")= chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_CACTTACATACT" "H7-CM-batch7-d0_CTGGAGCTGGCC" ... ..- attr(, "Size")= int 5721 ..- attr(, "call")= language as.dist.default(m = 1 - armaCor(t(pca.mat))) ..- attr(, "Diag")= logi FALSE ..- attr(*, "Upper")= logi FALSE

rvel.cd <- gene.relative.velocity.estimates(emat1,nmat1,deltaT=1,kCells=5,cell.dist=cell.dist,fit.quantile=0.02,verbose=T,n.cores=10) matching cells between cell.dist and emat/nmat ... done calculating cell knn ... done calculating convolved matrices ... done fitting gamma coefficients ... `

Thanks a lot, Peng

[pkharchenko]

Not sure why it’s hanging (is it using CPU?) I would also try more stringent gene filtering (down to ~900 genes), which should run as fast as the tutorial.

-peter.

On Aug 22, 2018, at 20:43, wulabupenn notifications@github.com wrote:

@pkharchenko Thank you for the prompt reply. Here is detail information about the data set:

`> str(emat1) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:8088910] 0 1 2 3 4 6 7 8 9 11 ... ..@ p : int [1:5722] 0 2242 4576 6716 9034 11191 13480 15724 17943 20229 ... ..@ Dim : int [1:2] 3139 5721 ..@ Dimnames:List of 2 .. ..$ : chr [1:3139] "AASDHPPT" "AB019441.29" "ABBA01017803.1" "ABCE1" ... .. ..$ : chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_ATTTTTCTTTTN" "H7-CM-batch7-d0_CACTTACATACT" ... ..@ x : num [1:8088910] 2 19 1 5 1 5 1 19 1 3 ... ..@ factors : list()

str(nmat1) Formal class 'dgCMatrix' [package "Matrix"] with 6 slots ..@ i : int [1:5032887] 0 1 5 6 9 14 16 18 20 22 ... ..@ p : int [1:5722] 0 2020 3865 5920 7976 10072 11793 13979 16081 18112 ... ..@ Dim : int [1:2] 5464 5721 ..@ Dimnames:List of 2 .. ..$ : chr [1:5464] "A2ML1-AS1" "A2MP1" "AACS" "AAK1" ... .. ..$ : chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_ATTTTTCTTTTN" "H7-CM-batch7-d0_CACTTACATACT" ... ..@ x : num [1:5032887] 1 3 1 1 1 1 1 2 1 1 ... ..@ factors : list() str(cell.dist) Class 'dist' atomic [1:16362060] 0.112 0.225 0.171 0.192 0.123 ... ..- attr(, "Labels")= chr [1:5721] "H7-CM-batch7-d0_CGTCACGTTCCN" "H7-CM-batch7-d0_CGCGAATTGTGA" "H7-CM-batch7-d0_CACTTACATACT" "H7-CM-batch7-d0_CTGGAGCTGGCC" ... ..- attr(, "Size")= int 5721 ..- attr(, "call")= language as.dist.default(m = 1 - armaCor(t(pca.mat))) ..- attr(, "Diag")= logi FALSE ..- attr(*, "Upper")= logi FALSE

rvel.cd <- gene.relative.velocity.estimates(emat1,nmat1,deltaT=1,kCells=5,cell.dist=cell.dist,fit.quantile=0.02,verbose=T,n.cores=10) matching cells between cell.dist and emat/nmat ... done calculating cell knn ... done calculating convolved matrices ... done fitting gamma coefficients ... `

Thanks a lot, Peng

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[wulabupenn]

Not, it is not use CPU at all. Do you think it is possible you can help me trouble shoot if I can provide the link for downloading the input dataset. The size of data set is ~ 200 Mb. Thank you!

[pkharchenko]

Most likely it’s stalling due to some multithreading lock. Normally saving the workplace, restarting the session and trying to run the fit command again fixes the issue. Have you tried that?

On Aug 22, 2018, at 2:10 PM, wulabupenn notifications@github.com wrote:

Not, it is not use CPU at all. Do you think it is possible you can help me trouble shoot if I can provide the link for downloading the input dataset. The size of data set is ~ 200 Mb. Thank you!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/velocyto-team/velocyto.R/issues/38#issuecomment-415126891, or mute the thread https://github.com/notifications/unsubscribe-auth/ALT78pnteOOnRfFGiV_lZu34ndJbdrrsks5uTZ6cgaJpZM4WIFnO.

[wulabupenn]

I save the data and restart session. Still stop at that step:

> load("0822_velocity_run.Rdata")
> library("velocyto.R")
Loading required package: Matrix
> rvel.cd <- gene.relative.velocity.estimates(emat1,nmat1,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=0.02)
matching cells between cell.dist and emat/nmat ... done
calculating cell knn ... done
calculating convolved matrices ... done
fitting gamma coefficients ... 

[hahaschool]

UPD: After testing various methods, I finally solved the problem by temporally set thread number of MKL to 1 inside the function body inside mclapply. It seems this is a legacy problem of MRO: https://groups.google.com/forum/#!msg/rropen/MG98lXsFepo/3h6XIKl_CwAJ .

---

Hi,

I have also encountered the same issue when I run gene.relative.velocity.estimates() on the dataset processed from https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR6854062 . I am using MRO 3.5.1.

I precalculated loom file by using velocyto CLI, and followed "Dentate Gyrus / loom" tutorial. Then gene.relative.velocity.estimates() stalls if n.cores > 1 (it won't stall if runs sequentially).

The dataset information is like:

> str(emat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:4333037] 5 6 25 27 38 46 53 56 57 61 ...
  ..@ p       : int [1:4262] 0 715 1442 2306 3159 4155 5031 6602 7376 8316 ...
  ..@ Dim     : int [1:2] 4269 4261
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:4269] "Vcpip1" "Snhg6" "Cops5" "Arfgef1" ...
  .. ..$ : chr [1:4261] "10x81_1_vanilla:AACGTTGTCAGATAAGx" "10x81_1_vanilla:AAAGTAGCATCTACGAx" "10x81_1_vanilla:AAAGATGGTCGTTGTAx" "10x81_1_vanilla:AAAGATGAGATAGCATx" ...
  ..@ x       : num [1:4333037] 4 1 1 2 2 1 1 2 1 1 ...
  ..@ factors : list()

By adding some debugging statement inside the mclapply function block of gene.relative.velocity.estimates() method, I have found that only the calculation of n.cores genes will be done before stalling. I have restarted R numerous times and the stalling persists.

Your help would be really appreciated.

[MichaelPeibo]

UPD: After testing various methods, I finally solved the problem by temporally set thread number of MKL to 1 inside the function body inside mclapply. It seems this is a legacy problem of MRO: https://groups.google.com/forum/#!msg/rropen/MG98lXsFepo/3h6XIKl_CwAJ .

Hi,

I have also encountered the same issue when I run gene.relative.velocity.estimates() on the dataset processed from https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR6854062 . I am using MRO 3.5.1.

I precalculated loom file by using velocyto CLI, and followed "Dentate Gyrus / loom" tutorial. Then gene.relative.velocity.estimates() stalls if n.cores > 1 (it won't stall if runs sequentially).

The dataset information is like:

> str(emat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:4333037] 5 6 25 27 38 46 53 56 57 61 ...
  ..@ p       : int [1:4262] 0 715 1442 2306 3159 4155 5031 6602 7376 8316 ...
  ..@ Dim     : int [1:2] 4269 4261
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:4269] "Vcpip1" "Snhg6" "Cops5" "Arfgef1" ...
  .. ..$ : chr [1:4261] "10x81_1_vanilla:AACGTTGTCAGATAAGx" "10x81_1_vanilla:AAAGTAGCATCTACGAx" "10x81_1_vanilla:AAAGATGGTCGTTGTAx" "10x81_1_vanilla:AAAGATGAGATAGCATx" ...
  ..@ x       : num [1:4333037] 4 1 1 2 2 1 1 2 1 1 ...
  ..@ factors : list()

By adding some debugging statement inside the mclapply function block of gene.relative.velocity.estimates() method, I have found that only the calculation of n.cores genes will be done before stalling. I have restarted R numerous times and the stalling persists.

Your help would be really appreciated.

Hi @hahaschool , same issue with me. how exactly did you fix it?

My understanding is my Ubuntu's problem with parallel or something(now I have same issue on ScaleDataand JackStraw in Seurat, just stuck when I use these functions, since I tried to install a software and tried to install some dependecies...)

[hahaschool]

UPD: After testing various methods, I finally solved the problem by temporally set thread number of MKL to 1 inside the function body inside mclapply. It seems this is a legacy problem of MRO: https://groups.google.com/forum/#!msg/rropen/MG98lXsFepo/3h6XIKl_CwAJ . Hi, I have also encountered the same issue when I run gene.relative.velocity.estimates() on the dataset processed from https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR6854062 . I am using MRO 3.5.1. I precalculated loom file by using velocyto CLI, and followed "Dentate Gyrus / loom" tutorial. Then gene.relative.velocity.estimates() stalls if n.cores > 1 (it won't stall if runs sequentially). The dataset information is like:

> str(emat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:4333037] 5 6 25 27 38 46 53 56 57 61 ...
  ..@ p       : int [1:4262] 0 715 1442 2306 3159 4155 5031 6602 7376 8316 ...
  ..@ Dim     : int [1:2] 4269 4261
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:4269] "Vcpip1" "Snhg6" "Cops5" "Arfgef1" ...
  .. ..$ : chr [1:4261] "10x81_1_vanilla:AACGTTGTCAGATAAGx" "10x81_1_vanilla:AAAGTAGCATCTACGAx" "10x81_1_vanilla:AAAGATGGTCGTTGTAx" "10x81_1_vanilla:AAAGATGAGATAGCATx" ...
  ..@ x       : num [1:4333037] 4 1 1 2 2 1 1 2 1 1 ...
  ..@ factors : list()

By adding some debugging statement inside the mclapply function block of gene.relative.velocity.estimates() method, I have found that only the calculation of n.cores genes will be done before stalling. I have restarted R numerous times and the stalling persists. Your help would be really appreciated.

Hi @hahaschool , same issue with me. how exactly did you fix it?

My understanding is my Ubuntu's problem with parallel or something(now I have same issue on ScaleDataand JackStraw in Seurat, just stuck when I use these functions, since I tried to install a software and tried to install some dependecies...)

Hi @MichaelPeibo,

A dirty but simple solution might be calling setMklthreads(1) each time before your invoking any of those functions, or more thoroughly, prepare a vanilla R environment for those problematic packages (should not be hard if you use anaconda).

(I'm a little bit surprised that this problem is still not fixed in this package at this time, it seems most people still prefer vanilla R so the problem is not so prevalent.)

[MichaelPeibo]

Hi @hahaschool Thanks for your response. I install RevoUtilsMath, and command this in R:

library(RevoUtilsMath) setMKLthreads(1)

But I still got stuck at fitting gamma coefficients step. I do not see ant program is taking CPU when I command top.

I got size of emat and namt here.

dim(emat) [1] 2900 9349 dim(nmat) [1] 8864 9349

Did I do the right command as you described above?

[hahaschool]

Hi @hahaschool Thanks for your response. I install RevoUtilsMath, and command this in R:

library(RevoUtilsMath) setMKLthreads(1)

But I still got stuck at fitting gamma coefficients step. I do not see ant program is taking CPU when I command top.

I got size of emat and namt here.

dim(emat) [1] 2900 9349 dim(nmat) [1] 8864 9349

Did I do the right command as you described above?

Hi @MichaelPeibo,

It seems that for each slave R process there will still be non-1 MKL threads.

You could try:

1. export MKL_NUM_THREADS=1 in shell to make it a environment variable to disable MKL.
2. Do a very dirty hack to the code of velocyto.R: suppose you are using the method gene.relative.velocity.estimates, then you could find its code (getAnywhere() may help), and find blocks using mclapply and add setMKLthreads(1) inside those blocks. Then you could override the hacked code back to velocyto package (assignInNamespace() may help).

An example:

...
if (is.null(old.fit)) {
      cat("fitting smat-based offsets ... \n")
      sfit <- data.frame(do.call(rbind, pbmclapply(velocyto.R:::sn(rownames(conv.emat.norm)), 
                                                           function(gn) {
                                                             setMKLthreads(1); # Hack here
                                                             df <- data.frame(n = (conv.nmat.norm[gn, steady.state.cells]), 
                                                                              e = (conv.emat.norm[gn, steady.state.cells]), 
                                                                              s = conv.smat.norm[gn, steady.state.cells])
                                                             sd <- lm(n ~ s, data = df)
                                                             r <- with(df[df[["s"]] > 0, ], cor(n, s, method = "spearman"), 
                                                                       3)
                                                             return(c(o = pmax(0, as.numeric(sd[["coef"]][1])), 
                                                                      s = as.numeric(sd[["coef"]][2]), r = r))
                                                           }, mc.cores = n.cores, mc.preschedule = F)))
      cat("done\n")
    }
    ...

I actually used 2. to solve the issue (also you could hack more, e.g. add a progress bar by using package pbmcapply) but as it is too dirty you may not wish to really use it if 1. could work.

[MichaelPeibo]

Hi @hahaschool , Thanks for your detailed demonstration.

After I installed RevoUtilsMath(some intel R?), I got crashed Seurat and velocyto.R. After I tried to reinstall them, I lost RevoUtilsMath, it seems these packages have some crash? So for your setMKLthreads(1) solution unfortunately did not work for me.

Actually, I got through by setting n.cores=1, I do not understand deep reason of it, it seems that my linux system are not working with parallel? like you referred that google group issue, it will make threads sleep not work.

Thanks!

[hahaschool]

Hi @MichaelPeibo,

Maybe your MRO is not using MKL as BLAS so MKL settings are not working. But anything banning multithreading in BLAS should work in this case. It seems to be some historical problems that mclapply is not compatible with multi-core BLAS and I am also not knowing about low-level reasons in detail.

Would be happy to see if you have any updates investigating such problems. Thanks.

[ice4prince]

I found this issue has been open for 3 years. I encountered the same issue. Anyone has already fixed it?