Error reading Velocity matrix

[romanhaa]

Hi Alex,

Following your question on Twitter I started a comparison for the spliced/unspliced counts from velocyto, STARsolo and kallisto. Unfortunately, when trying to read in the "Velocyto" matrix generated by STARsolo, I'm getting this error:

Matrix::readMM(".../STARsolo/Velocyto/raw/matrix.mtx")
# Error: readMM(): row     values 'i' are not in 1:nr

Any idea what the issue could be? The mtx file is ~1 GB in size.

Using scanpy's read_mtx function works fine.

import scanpy as sc
sc.read_mtx(".../STARsolo/Velocyto/raw/matrix.mtx")
# AnnData object with n_obs × n_vars = 11720 × 60609

For those running into the same problem, you can load the .mtx file generated by STARsolo and then save it again using scanpy/scipy functions as shown below. Afterwards, it will be readable by Matrix::readMM() in R.

import scanpy as sc
import scipy
t = sc.read_mtx('.../STARsolo/Velocyto/raw/matrix.mtx')
scipy.io.mmwrite('.../STARsolo/Velocyto/raw/matrix_new_format.mtx', t.X)

Matrix::readMM(".../STARsolo/Velocyto/raw/matrix_new_format.mtx") %>% str()
# Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
#   ..@ i       : int [1:48570060] 1 1 1 1 1 1 1 1 1 1 ...
#   ..@ j       : int [1:48570060] 6962 54995 69331 89585 121706 145759 228188 337667 383554 383896 ...
#   ..@ Dim     : int [1:2] 60609 6794880
#   ..@ Dimnames:List of 2
#   .. ..$ : NULL
#   .. ..$ : NULL
#   ..@ x       : num [1:48570060] 0 0 0 0 0 0 0 0 0 0 ...
#   ..@ factors : list()

Best,
Roman

# R version 3.6.1 (2019-07-05)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Debian GNU/Linux 9 (stretch)
# 
# Matrix products: default
# BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
# 
# locale:
#  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
#  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] Matrix_1.2-17   forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3    
#  [5] purrr_0.3.3     readr_1.3.1     tidyr_1.0.0     tibble_2.1.3   
#  [9] ggplot2_3.2.1   tidyverse_1.2.1
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.2       cellranger_1.1.0 pillar_1.4.2     compiler_3.6.1  
#  [5] tools_3.6.1      zeallot_0.1.0    jsonlite_1.6     lubridate_1.7.4 
#  [9] lifecycle_0.1.0  gtable_0.3.0     nlme_3.1-141     lattice_0.20-38 
# [13] pkgconfig_2.0.3  rlang_0.4.1      cli_1.1.0        rstudioapi_0.10 
# [17] haven_2.1.1      withr_2.1.2      xml2_1.2.2       httr_1.4.1      
# [21] generics_0.0.2   vctrs_0.2.0      hms_0.5.2        grid_3.6.1      
# [25] tidyselect_0.2.5 glue_1.3.1       R6_2.4.0         readxl_1.3.1    
# [29] modelr_0.1.5     magrittr_1.5     backports_1.1.5  scales_1.0.0    
# [33] rvest_0.3.4      assertthat_0.2.1 colorspace_1.4-1 stringi_1.4.3   
# [37] lazyeval_0.2.2   munsell_0.5.0    broom_0.5.2      crayon_1.3.4

[romanhaa]

Now that I was able to load the data, I'm asking myself what this matrix is. Shouldn't there have been spliced, unspliced, and ambiguous? Sorry, maybe I'm completely on the wrong track.

[alexdobin]

Hi Roman,

thanks for looking into it... it seems you are going to be the alpha-tester for this feature! :)
Sorry for the lack of documentation, and for the breaking of the matrix market format.
Basically, I dumped the spliced, unspliced and ambiguous counts as three columns in this file
i.e. col3=spliced, col4=unspliced, col5=ambiguous

Cheers
Alex

[romanhaa]

Ohhhh, I see. I guess in that case I can reproduce the sparse matrix like this:

t <- readr::read_delim('matrix.mtx', delim = ' ', skip = 3, col_names = FALSE)

spliced <- Matrix::sparseMatrix(
  i = t$X1,
  j = t$X2,
  x = t$X3,
  dims = c(60609, 6794880)
)
Matrix::writeMM(spliced, 'spliced.mtx')

unspliced <- Matrix::sparseMatrix(
  i = t$X1,
  j = t$X2,
  x = t$X4,
  dims = c(60609, 6794880)
)
Matrix::writeMM(unspliced, 'unspliced.mtx')

The matrix dimensions are taken from the header of the matrix.mtx file. It's looking good so far. No idea what the scanpy method gave me then :D

[romanhaa]

I see high correlation between the results of velocyto (input was the Cell Ranger .bam file) and STARsolo, both for spliced and unspliced counts. Less so for kallisto. There are some interesting differences in the number of transcripts and detected genes.

Number of transcripts:

Number of detected genes:

Correlation of spliced counts (by genes; 0's are set to 0.1 for log-scale):

Correlation of unspliced counts (by genes; 0's are set to 0.1 for log-scale):

[alexdobin]

Hi Roman,

thanks! Looks similar to what I was getting for STARsolo.
I could not fully reproduce the velocyto logic, unfortunately. I was communicating with Gioele La Manno about it, but I guess we both got busy and dropped the ball. I am planning to start working on it again. If you make these results public, it will put a bit more pressure on us. :)

When I tried to run kallisto with velocyto calculation, it used 68GB of RAM.
I am going to ask Lior about it on tweeter. :)

Cheers
Alex

[romanhaa]

Yes, I'd like to know Lior's opinion as well. kallisto took a lot longer to run this time. I have no reliable numbers (will have to repeat) but kallisto took a significant hit in performance for this analysis. Generating the index alone took 2.5 hours (and it is 42GB in size). STARsolo didn't need a new index and finished the analysis in 1 hour with 20 threads (I get 9.5 CPU hours in total).

If you make these results public, it will put a bit more pressure on us.

I'll do my best but I'm not sure which form would be appropriate.

Cheers,
Roman

[alexdobin]

Hi Roman,

If without velocyto STARsolo took 35 min, and with 1 hour - it seems reasonable. It's a bit slower than my expectation (~50% increase in time) - of course, the actual results depend strongly on the hardware. The velocyto counting is not parallelized yet, I hope to improve it in the future.

I am afraid to give advice about the form or forum. Tweeter discussions tend to become overheated... some people just cannot take others criticizing their work (not me : ) - I am happy to hear any kind of critique, as it always helps to improve the results).
I think you may have enough material for a benchmarking paper of your own. It would probably be best if you contacted all Authors for optimal parameters etc.
We can also have an e-mail discussion with Velocyto Authors. Gioele La Manno was the 1st Author, and Peter Kharchenko and Sten Linnarson the senior Authors.

Cheers
Alex

[romanhaa]

Yes the run time is absolutely reasonable. Actually, for the spliced/unspliced counts STARsolo appears to be the fastest of all three methods. But further improvements would of course be appreciated :)

I appreciate your transparency and open-minded attitude towards these comparisons. Indeed, Twitter discussions are not ideal. I don't want to side with anybody and just present my perspective as a user. I'll have to think about the possibility of putting together a paper. Thanks for your suggestions and the offer to reach out to Gioele La Manno and the other authors of Velocyto. Optimizing the parameters for each tool would probably be an important thing to do before publishing performance comparisons.

Btw, what do you think about the difference in ratio (spliced/unspliced counts) between STARsolo vs Velocyto + kallisto? I found that to be quite striking.

Cheers,
Roman

[alexdobin]

Hi Roman,

sorry for the belayed reply, I was swamped by the Genome Informatics conference.

I appreciate your transparency and open-minded attitude towards these comparisons. Indeed, Twitter discussions are not ideal. I don't want to side with anybody and just present my perspective as a user. I'll have to think about the possibility of putting together a paper. Thanks for your suggestions and the offer to reach out to Gioele La Manno and the other authors of Velocyto. Optimizing the parameters for each tool would probably be an important thing to do before publishing performance comparisons.

Absolutely!
And I would welcome any feedback from the user perspective.

Btw, what do you think about the difference in ratio (spliced/unspliced counts) between STARsolo vs Velocyto + kallisto? I found that to be quite striking.

Indeed, the difference seems on the large sides. I have not yet looked at kallisto/bustools results carefully. I have made my comparisons with velocyto, and I remember that the total counts: spliced+unspliced+ambiguous were very close, but there was some exchange of counts between spliced and ambiguous. Just to clarify, by the number of transcripts you mean "number of molecules" - i.e. UMIs? And you are comparing the counts for the ~10k cells that were filtered by velocyto?

Cheers
Alex

[jenzopr]

Just for reference/completion, here is a python way to do that:

import numpy as np
from scipy import sparse
mtx = np.loadtxt('matrix.mtx', skiprows=3, delimiter=' ')

spliced = sparse.coo_matrix((mtx[:,2], (mtx[:,0]-1, mtx[:,1]-1)), shape = (60558,737280))
unspliced = sparse.coo_matrix((mtx[:,3], (mtx[:,0]-1, mtx[:,1]-1)), shape = (60558,737280))
ambiguous = sparse.coo_matrix((mtx[:,4], (mtx[:,0]-1, mtx[:,1]-1)), shape = (60558,737280))

The shape (60558,737280) is taken from the matrix header
Thanks to @msbentsen :-)

[droplet-lab]

Thanks @jenzopr, I am still quite new to python, but was wondering how you merged this to the adata gene count matrix for velocity analysis using scVelo?

[jenzopr]

@Jachimdejonghe scvelo has a nice utility function for this: scvelo.utils.merge, find the help here