Open Jemkon opened 6 years ago
I am getting the same error. I dug into the code and the internalcc
variable calculated is all 'not a number', i.e. NaN
. Check out the output from my debugger
Browse[2]> str(cc) num [1:6004, 1:6004] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... Browse[2]> all(is.nan(cc)) [1] TRUE
Note the cc
variable is calculated by colDeltaCor(em, (sqrt(abs(nd)) * sign(nd)), nthreads = n.cores)
. The transition probability tp
variable is a function of cc
so it too becomes a NaN
. These errors eventually cause the seq(from, to, length.out)
command to fail because the range variables rx
are NaN
.
I am stuck too. Help please.
Hi guys,
I had the same error and then realise cell names in the embedding were not matching those in the velocity object. In this case, cc is NaN because em and emb are empty. That was it.
Best, Ivan
Thanks for the hint, Ivan!
If others keep encountering this error, please share an example dataset to help us replicate it.
Thank you, -peter.
On May 23, 2018, at 15:23, rosshandler notifications@github.com wrote:
Hi guys,
I had the same error and then realise cell names in the embedding were not matching those in the velocity object. In this case, cc is NaN because em and emb are empty. That was it.
Best, Ivan
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.
Ivan, thank you for your reporting; however, this is not the same problem I have. Further investigation shows the deltaE
object (unscaled transcriptional change ) returned from gene.relative.velocity.estimates
is virtually all zeros (values are about 2e-15 at largest). This result is used as nd
in the colDeltaCor(em, (sqrt(abs(nd)) * sign(nd)), nthreads = n.cores)
call. The NaN
values are a result of nd
having a standard deviation of zero for most of its columns when evaluated for its correlation to em
.
Hello all,
I am getting the same error :
`
emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.2) nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05) length(intersect(rownames(emat),rownames(nmat))) [1] 1357 fit.quantile <- 0.02 rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,fit.quantile=fit.quantile) calculating cell knn ... done calculating convolved matrices ... done fitting gamma coefficients ... done. succesfful fit for 1357 genes filtered out 61 out of 1357 genes due to low nmat-emat correlation filtered out 79 out of 1296 genes due to low nmat-emat slope calculating RNA velocity shift ... done calculating extrapolated cell state ... done show.velocity.on.embedding.cor(emb,rvel.cd,n=200,scale='sqrt',cex=0.8,arrow.scale=3,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1) delta projections ... sqrt knn ... transition probs ... done calculating arrows ... done grid estimates ... Error in seq.default(rx[1], rx[2], length.out = grid.n) : 'from' must be a finite number`
Any help is greatly appreciated. Let me know if you need me to share my data and what would be the best way to do so.
Thanks, Fiona
Same problem here.
show.velocity.on.embedding.cor(fourh.emb,fourh.rvel.cd,n=300,scale='sqrt',cell.colors=ac(fourh.cell.colors,alpha=0.5),cex=0.8,arrow.scale=5,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1) delta projections ... sqrt knn ... transition probs ... done calculating arrows ... done grid estimates ... Error in seq.default(rx[1], rx[2], length.out = grid.n) : 'from' must be a finite number
Let me know how best to share my data so you can troubleshoot.
Hi,
we can do this: emat <- emat[,colSums(is.na(emat))<nrow(emat)]
to remove all NAs...
The script runs smoothly then.
I have successfully run velocyto on one dataset. When I switch to another dataset, I get the same error...
Let me know if any other information you need me to share. Thank you so much!
Hi, I'm also getting the same error. I've run velocyto successfully on another dataset but I keep getting an error with this dataset. I split it in half and one subset works where as the other subset has the same error. This dataset has 30 000 cells.
Error: delta projections ... sqrt knn ... transition probs ... done calculating arrows ... done grid estimates ... Error in seq.default(rx[1], rx[2], length.out = grid.n) : 'from' must be a finite number Calls: show.velocity.on.embedding.cor -> seq -> seq.default
I tried the NA remove but it's not resolving the problem:
emat2 <- emat[,colSums(is.na(emat))<nrow(emat)] emat3 <- emat[,colSums(is.infinite(emat))<nrow(emat)] emat4 <- emat[,colSums(is.finite(emat))<nrow(emat)]
dim(emat) [1] 723 29943 dim(emat2) [1] 723 29943 dim(emat3) [1] 723 29943 dim(emat4) [1] 723 0
All numbers are finite, there are no Na, NaN, or infinites to remove. I also double checked the nmat object doesn't have na errors.
Thanks you so much in advanced, please let me know if any other information can help trouble shoot.
Best, -Frances
I got the same error.
But once I simply changed the average
parameter of filter.genes.by.cluster.expression
, It worked.
though I do not know why...
Takahiro
Digging in: The error occurs because cc <- colDeltaCorLog10(...)
returns a matrix with a values that are NA, which causes Matrix::colSums(tp)
to be all-NaN as well, ultimately causing embArrows(coords, tp, ...)
to be all-NA.
The relevant lines:
The NA values occur for me because in vel$current
, a cell had all-NaN values: vel$current[, 'SCG71:AAGTTGTCACCTTGCC'] == c(NaN, NaN, ...)
OK, there’s more problems. If I eliminate this cell, I still get NAs in the result of colDeltaCorLog10
and friends, which propagate to tp
and then to the result of arsd <- embArrows(...)
.
I can’t say I really understand colDeltaCor*
, and if there’s a bug…
colDeltaCor functions calculate for each cell (column E_i of a matrix E) correlation between E-E_i and D_i (difference vector expected from velocity of the cell I). colDeltaSqrt and Log10 variants calculate correlation on the sqrt or log10 transformed values of E-E_i. This is done to figure out towards what cells in the neighborhood of a cell i should the velocity D_i be projected.
I am not sure what’s causing the issue here ... perhaps you have some duplicated cells. Is it just the diagonal values of the return matrix that are NaN? Can you share an example?
On Jul 17, 2019, at 06:07, Philipp A. notifications@github.com wrote:
OK, there’s more problems. If I eliminate this cell, I still get NAs in the result of colDeltaCorLog10, which propagate to tp and then to
I can’t say I really understand colDeltaCorLog10, and if there’s a bug…
— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.
Not just the diagonals. The cells where those NAs occur are mostly all-zero in E, and the few that aren’t are equal:
na_idx <- as.data.frame(which(is.na(cc), arr.ind = TRUE))
cbind(
apply(as.matrix(vel$current)[, na_idx$row], 2, sum),
apply(as.matrix(vel$current)[, na_idx$col], 2, sum)
)
Yes, the if there are two identical cells (i.e. two all-0 E columns), that would lead to NAs … most likely those entries in the output matrix can be simply replaced with 0s, but let me look at it.
On Jul 17, 2019, at 06:53, Philipp A. notifications@github.com wrote:
Not just the diagonals. The cells where those NAs occur are mostly all-zero in E:
na_idx <- as.data.frame(which(is.na(cc), arr.ind = TRUE )) cbind( apply(as.matrix( vel$current)[, na_idx$row], 2, sum ), apply(as.matrix( vel$current)[, na_idx$col], 2, sum ) )
Output — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.
In the code, em
is gene.relative.velocity.estimates(...)$current
, the “current normalized expression state”.
I’m using your example data from http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom In the raw data, there’s no cells with identical expression profile.
Therefore the cells are only identical after gene.relative.velocity.estimates
did its thing. How does that happen? Did I miss a preprocessing step?
In any case, the function should warn if it’s returning junk that can’t be further processed. I spent far too much time debugging this when gene.relative.velocity.estimates
could simply raise a warning
telling me about $current
containing all-NA cells and duplicate cells now, and that I should use $valid_cells
to filter emat
, nmat
, and cell.dist
and try again.
If I hackily remove duplicated and NA cells in a loop, it will work. It looks weird though:
Original | Mine |
---|---|
![]() |
![]() |
I am having the same issue.
However, if I pre-filter out cells more stringently based on nFeatures and nUMI the problem has gone away. Thanks!
I have the same issue here. I've just started using velocyto.R with the seurat wrapper, not sure how to trouble shoot as it is all automated
I gave up and switched to scvelo. I don’t want to advertise as it’s from my lab and I’m biased, so I’m only saying that it worked immediately, and with more reasonable results (Probably because I forgot lib size normalization when trying velocyto). If someone wants to try my hack to eliminate problematic cells, here it is:
velocity_estimates_hack <- function(spliced, unspliced, avg_spliced, avg_unspliced, dists, ...)
emat <- filter.genes.by.cluster.expression(spliced, cluster_labels, min.max.cluster.average = avg_spliced)
nmat <- filter.genes.by.cluster.expression(unspliced, cluster_labels, min.max.cluster.average = avg_unspliced)
dists <- as(dists, 'sparseMatrix')
while (TRUE) {
vel <- gene.relative.velocity.estimates(emat, nmat, cell.dist = as.dist(dists), ...)
idx_avail <- apply(vel$current, 2, function(cell) all(!is.na(cell)))
idx_uniqu <- !duplicated(as.matrix(vel$current), MARGIN = 2)
keep_cells <- idx_avail & idx_uniqu
if (all(keep_cells)) break
if (!all(idx_avail)) message('Eliminating NA cells: ', paste(colnames(emat)[!idx_avail], collapse = ', '))
if (!all(idx_uniqu)) message('Eliminating duplicate cells: ', paste(colnames(emat)[!idx_uniqu], collapse = ', '))
emat <- emat[, keep_cells, drop = FALSE]
nmat <- nmat[, keep_cells, drop = FALSE]
dists <- dists[keep_cells, keep_cells]
}
list(emat = emat, nmat = nmat, dists = dists, vel = vel)
}
Thank you, Philipp for the suggested patch.
Unfortunately, I can’t seem to replicate the error. As suggested, I’ve ran the SCG71 tutorial example (http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html) on several R versions that I had locally (R3.6.1, R3.5.3 and R3.5.1) and in each case it worked fine.
If someone could point me at an example that may be easier to reproduce, I would appreciate it!
Thanks!
On Jul 19, 2019, at 03:25, Philipp A. notifications@github.com wrote:
I gave up and switched to scvelo. I don’t want to advertise as it’s from my lab and I’m biased, so I’m only saying that it worked immediately, and with more reasonable results. If someone wants to try my hack to eliminate problematic cells, here it is:
velocity_estimates_hack <- function(spliced, unspliced, avg_spliced, avg_unspliced, dists, ... )
emat <- filter.genes.by.cluster.expression(spliced, cluster_labels, min.max.cluster.average = avg_spliced )
nmat <- filter.genes.by.cluster.expression(unspliced, cluster_labels, min.max.cluster.average = avg_unspliced )
dists <- as(dists, 'sparseMatrix' )
while (TRUE ) {
vel <- gene.relative.velocity.estimates(emat, nmat, cell.dist = as.dist(dists), ... )
idx_avail <- apply(vel$current, 2, function(cell) all(!is.na(cell )))
idx_uniqu <- !duplicated(as.matrix(vel$current), MARGIN = 2 )
keep_cells <- idx_avail & idx_uniqu
if (all(keep_cells)) break
if (!all(idx_avail)) message('Eliminating NA cells: ', paste(colnames(emat)[!idx_avail], collapse = ', ' ))
if (!all(idx_uniqu)) message('Eliminating duplicate cells: ', paste(colnames(emat)[!idx_uniqu], collapse = ', ' ))
emat <- emat[, keep_cells, drop = FALSE ]
nmat <- nmat[, keep_cells, drop = FALSE ]
dists <- dists[keep_cells, keep_cells ] }
list(emat = emat, nmat = nmat, dists = dists, vel = vel ) }
— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.
Having the same issue and would appreciate any help.
@pkharchenko, happy to share the data confidentially to help you track down what the issue could be. Im trying to get velocyto to work on an indrop dataset.
Can you please paste the error traceback() and sessionInfo(). Thank you, -peter.
On Jul 22, 2019, at 17:33, mbassalbioinformatics notifications@github.com wrote:
Having the same issue and would appreciate any help. @pkharchenko, happy to share the data confidentially to help you track down what the issue could be. Im trying to get velocyto to work on an indrop dataset.
— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.
traceback()
4: stop("'from' must be a finite number")
3: seq.default(rx[1], rx[2], length.out = grid.n)
2: seq(rx[1], rx[2], length.out = grid.n)
1: show.velocity.on.embedding.cor(emb_matrix, rvel.cd, n = 200,
scale = "sqrt", cell.colors = ac(cell.colors, alpha = 0.5),
cex = 0.8, arrow.scale = 3, show.grid.flow = TRUE, min.grid.cell.mass = 0.5,
grid.n = 40, arrow.lwd = 1, do.par = F, cell.border.alpha = 0.1,
n.cores = 35)
sessioninfo()
Rversion 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
[6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] SparseData_1.0.0 Biobase_2.44.0 BiocGenerics_0.30.0 forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
[8] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3 ggplot2_3.2.0 tidyverse_1.2.1 pagoda2_0.1.0 igraph_1.2.4.1
[15] velocyto.R_0.6 Matrix_1.2-17
loaded via a namespace (and not attached):
[1] httr_1.4.0 bit64_0.9-7 jsonlite_1.6 splines_3.6.1 hdf5r_1.2.0 modelr_0.1.4 assertthat_0.2.1 triebeard_0.3.0
[9] urltools_1.7.3 cellranger_1.1.0 remotes_2.1.0 pillar_1.4.2 backports_1.1.4 lattice_0.20-38 glue_1.3.1 rvest_0.3.4
[17] colorspace_1.4-1 pkgconfig_2.0.2 broom_0.5.2 haven_2.1.1 scales_1.0.0 processx_3.4.1 brew_1.0-6 mgcv_1.8-28
[25] generics_0.0.2 withr_2.1.2 lazyeval_0.2.2 cli_1.1.0 magrittr_1.5 crayon_1.3.4 readxl_1.3.1 ps_1.3.0
[33] fansi_0.4.0 nlme_3.1-140 MASS_7.3-51.4 xml2_1.2.0 pkgbuild_1.0.3 Rook_1.1-1 tools_3.6.1 prettyunits_1.0.2
[41] hms_0.5.0 munsell_0.5.0 cluster_2.1.0 irlba_2.3.3 callr_3.3.1 pcaMethods_1.76.0 compiler_3.6.1 rlang_0.4.0
[49] grid_3.6.1 rstudioapi_0.10 rjson_0.2.20 base64enc_0.1-3 gtable_0.3.0 curl_3.3 R6_2.4.0 lubridate_1.7.4
[57] bit_1.1-14 zeallot_0.1.0 utf8_1.1.4 rprojroot_1.3-2 stringi_1.4.3 Rcpp_1.0.1 vctrs_0.2.0 tidyselect_0.2.5
@pkharchenko, any updates regarding this?
I would also appreciate any updates about this error. Using the Seurat wrapper for velocyto.R. @pkharchenko
getting same error, any news?
@pkharchenko
I am also having same error. Any suggestions ?
emat <- emat[,colSums(is.na(emat))<nrow(emat)] is not solving the issue
Keep having the same error, even though there is no NaNs or NAs, can it be related to 5'chemistry?
Hi, I was having the same error when running 'show.velocity.on.embedding.cor'
delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done
grid estimates ... Error in seq.default(rx[1], rx[2], length.out = grid.n) :
'from' must be a finite number
From the different comments above I thought about eliminating any duplicated columns from emat and nmat. I did the following:
emat<-emat[,-which(duplicated(as.matrix(emat), MARGIN = 2))]
nmat<-nmat[,-which(duplicated(as.matrix(nmat), MARGIN = 2))]
For the newbies like me,
Before running I would recommend to check what
which(duplicated(as.matrix(emat), MARGIN = 2))
gives back for emat and nmat. In my case it was a dozen columns=cells and it was exactly the same columns for both emat and nmat.
I hope this can help some of you but use with caution since I am not an expert and I have not fully checked what these columns are.
Just the same issue. I just modified the code of function show.velocity.on.embedding.cor to remove NA of ars and arsd. And it worked. ars = na.omit(ars) arsd = na.omit(arsd)
No idea if there are other problems.
I've found the same error, even though there are no NANs or NAs in my data. Any suggestions, please?
I had the same issue and solved it by filtering cells by nFeature_spliced more stringently before running SCTransform. Hope it helps!
@vickyrev I found the same to help for me. Thanks!
@pkharchenko In my case, all the matrices become NA after the following statements. I changed two things. (1) Matrix::colSums(tp) This function needs to be called with a "na.rm = TRUE" parameter in order to avoid NA-only result. (2) Any NA values in the sparse matrix should be converted to zero.
ORIGINAL
cat("transition probs ... ")
tp <- exp(cc/corr.sigma) * emb.knn
tp <- t(t(tp)/Matrix::colSums(tp)) # Matrix::colSums(tp) makes every element NA
tp <- as(tp, "dgCMatrix")
cat("done\n")
...
# without replacing NAs with zeros, arsd becomes an NA matrix
arsd <- data.frame(t(**embArrows(emb, tp, arrow.scale,
n.cores)**))
CHANGED
cat("transition probs ... ")
tp <- exp(cc/corr.sigma) * emb.knn
tp <- t(t(tp)/Matrix::colSums(tp, na.rm =TRUE))
tp <- as.matrix(tp)
tp[is.na(tp)] <- 0
tp <- as(tp, "dgCMatrix")
cat("done\n")
I hope that this helps anyone who encounters this problem.
I am having the same issue.
However, if I pre-filter out cells more stringently based on nFeatures and nUMI the problem has gone away. Thanks!
This works for me. Thanks!
Very weird...
Nothing from the above answers works for me but I accidentally solve the problem simply by adding the argument "n.cores = 10" in show.velocity.on.embedding.cor()
......
@pkharchenko In my case, all the matrices become NA after the following statements. I changed two things. (1) Matrix::colSums(tp) This function needs to be called with a "na.rm = TRUE" parameter in order to avoid NA-only result. (2) Any NA values in the sparse matrix should be converted to zero.
ORIGINAL
cat("transition probs ... ") tp <- exp(cc/corr.sigma) * emb.knn tp <- t(t(tp)/Matrix::colSums(tp)) # Matrix::colSums(tp) makes every element NA tp <- as(tp, "dgCMatrix") cat("done\n") ... # without replacing NAs with zeros, arsd becomes an NA matrix arsd <- data.frame(t(**embArrows(emb, tp, arrow.scale, n.cores)**))
CHANGED
cat("transition probs ... ") tp <- exp(cc/corr.sigma) * emb.knn tp <- t(t(tp)/Matrix::colSums(tp, na.rm =TRUE)) tp <- as.matrix(tp) tp[is.na(tp)] <- 0 tp <- as(tp, "dgCMatrix") cat("done\n")
I hope that this helps anyone who encounters this problem. This works for me. Thanks!
Dear all, we also faced this issue during our analysis. After detailed debugging of the problem, I found the issue which comes from the part where correlation coefficients are calculated. Please check your data to see if any of the samples from spliced and unspliced matrixes have 0 standard deviations (most probably nonfiltered cells, which have all 0s). This 0 SD samples lead to NaN values from correlation calculation, which leads to the same error "'from' must be a finite number".
You need to either remove those 0 SD cells or set all NaN values to 0 in arsd table in the body of show.velocity.on.embedding.cor function.
Dear all, we also faced this issue during our analysis. After detailed debugging of the problem, I found the issue which comes from the part where correlation coefficients are calculated. Please check your data to see if any of the samples from spliced and unspliced matrixes have 0 standard deviations (most probably nonfiltered cells, which have all 0s). This 0 SD samples lead to NaN values from correlation calculation, which leads to the same error "'from' must be a finite number".
You need to either remove those 0 SD cells or set all NaN values to 0 in arsd table in the body of show.velocity.on.embedding.cor function.
In response to this suggestion, I have modified the source code and made it available for everyone to use. source("./show.velocity.on.embedding.cor.my.R") environment(show.velocity.on.embedding.cor.my) <- environment(show.velocity.on.embedding.cor) p1 <- show.velocity.on.embedding.cor.my(ldatumap,rvel.cd,n=30,scale='sqrt',cell.colors=ac(colors,alpha=0.5), cex=0.8,arrow.scale=2,show.grid.flow=T,min.grid.cell.mass=1.0,grid.n=50,arrow.lwd=1,do.par=F, cell.border.alpha =0.1,USE_OPENMP=1,n.cores=50,main="Cell Velocity")
show.velocity.on.embedding.cor.my <- function (emb, vel, n = 100, cell.colors = NULL, corr.sigma = 0.05, show.grid.flow = FALSE, grid.n = 20, grid.sd = NULL, min.grid.cell.mass = 1, min.arrow.size = NULL, arrow.scale = 1, max.grid.arrow.length = NULL, fixed.arrow.length = FALSE, plot.grid.points = FALSE, scale = "log", nPcs = NA, arrow.lwd = 1, xlab = "", ylab = "", n.cores = defaultNCores(), do.par = T, show.cell = NULL, cell.border.alpha = 0.3, cc = NULL, return.details = FALSE, expression.scaling = FALSE, ...) { randomize <- FALSE if (do.par) par(mfrow = c(1, 1), mar = c(3.5, 3.5, 2.5, 1.5), mgp = c(2, 0.65, 0), cex = 0.85) celcol <- "white" if (is.null(show.cell)) { celcol <- cell.colors[rownames(emb)] } plot(emb, bg = celcol, pch = 21, col = ac(1, alpha = cell.border.alpha), xlab = xlab, ylab = ylab, ...) em <- as.matrix(vel$current) ccells <- intersect(rownames(emb), colnames(em)) em <- em[, ccells] emb <- emb[ccells, ] nd <- as.matrix(vel$deltaE[, ccells]) cgenes <- intersect(rownames(em), rownames(nd)) nd <- nd[cgenes, ] em <- em[cgenes, ] if (randomize) { nd <- t(apply(nd, 1, function(x) (rbinom(length(x), 1, 0.5) 2 - 1) abs(sample(x)))) } if (is.null(cc)) { cat("delta projections ... ") if (scale == "log") { cat("log ") cc <- colDeltaCorLog10(em, (log10(abs(nd) + 1) sign(nd)), nthreads = n.cores) } else if (scale == "sqrt") { cat("sqrt ") cc <- colDeltaCorSqrt(em, (sqrt(abs(nd)) sign(nd)), nthreads = n.cores) } else if (scale == "rank") { cat("rank ") cc <- colDeltaCor((apply(em, 2, rank)), (apply(nd, 2, rank)), nthreads = n.cores) } else { cat("linear ") cc <- colDeltaCor(em, nd, nthreads = n.cores) } colnames(cc) <- rownames(cc) <- colnames(em) diag(cc) <- 0 } cat("knn ... ") if (n > nrow(cc)) { n <- nrow(cc) } emb.knn <- balancedKNN(t(emb), k = n, maxl = nrow(emb), dist = "euclidean", n.threads = n.cores) diag(emb.knn) <- 1 cat("transition probs ... ") tp <- exp(cc/corr.sigma) emb.knn tp <- t(t(tp)/Matrix::colSums(tp)) tp <- as(tp, "dgCMatrix") cat("done\n") if (!is.null(show.cell)) { i <- match(show.cell, rownames(emb)) if (is.na(i)) stop(paste("specified cell", i, "is not in the embedding")) points(emb, pch = 19, col = ac(val2col(tp[rownames(emb), show.cell], gradient.range.quantile = 1), alpha = 0.5)) points(emb[show.cell, 1], emb[show.cell, 2], pch = 3, cex = 1, col = 1) di <- t(t(emb) - emb[i, ]) di <- di/sqrt(Matrix::rowSums(di^2)) arrow.scale di[i, ] <- 0 dir <- Matrix::colSums(di tp[, i]) dic <- Matrix::colSums(di (tp[, i] > 0)/sum(tp[, i] > 0)) dia <- dir - dic suppressWarnings(arrows(emb[colnames(em)[i], 1], emb[colnames(em)[i], 2], emb[colnames(em)[i], 1] + dic[1], emb[colnames(em)[i], 2] + dic[2], length = 0.05, lwd = 1, col = "blue")) suppressWarnings(arrows(emb[colnames(em)[i], 1], emb[colnames(em)[i], 2], emb[colnames(em)[i], 1] + dir[1], emb[colnames(em)[i], 2] + dir[2], length = 0.05, lwd = 1, col = "red")) suppressWarnings(arrows(emb[colnames(em)[i], 1] + dic[1], emb[colnames(em)[i], 2] + dic[2], emb[colnames(em)[i], 1] + dir[1], emb[colnames(em)[i], 2] + dir[2], length = 0.05, lwd = 1, lty = 1, col = "grey50")) suppressWarnings(arrows(emb[colnames(em)[i], 1], emb[colnames(em)[i], 2], emb[colnames(em)[i], 1] + dia[1], emb[colnames(em)[i], 2] + dia[2], length = 0.05, lwd = 1, col = "black")) } else { cat("calculating arrows ... ") arsd <- data.frame(t(embArrows(emb, tp, arrow.scale, n.cores))) arsd[is.na(arsd)] <- 0 rownames(arsd) <- rownames(emb) if (expression.scaling) { tpb <- tp > 0 tpb <- t(t(tpb)/colSums(tpb)) es <- as.matrix(em %% tp) - as.matrix(em %% as.matrix(tpb)) pl <- pmin(1, pmax(0, apply(as.matrix(vel$deltaE[, colnames(es)]) es, 2, sum)/sqrt(colSums(es es)))) arsd <- arsd pl } ars <- data.frame(cbind(emb, emb + arsd)) colnames(ars) <- c("x0", "y0", "x1", "y1") colnames(arsd) <- c("xd", "yd") rownames(ars) <- rownames(emb) cat("done\n") if (show.grid.flow) { cat("grid estimates ... ") rx <- range(c(range(ars$x0), range(ars$x1))) ry <- range(c(range(ars$y0), range(ars$y1))) gx <- seq(rx[1], rx[2], length.out = grid.n) gy <- seq(ry[1], ry[2], length.out = grid.n) if (is.null(grid.sd)) { grid.sd <- sqrt((gx[2] - gx[1])^2 + (gy[2] - gy[1])^2)/2 cat("grid.sd=", grid.sd, " ") } if (is.null(min.arrow.size)) { min.arrow.size <- sqrt((gx[2] - gx[1])^2 + (gy[2] - gy[1])^2) 0.01 cat("min.arrow.size=", min.arrow.size, " ") } if (is.null(max.grid.arrow.length)) { max.grid.arrow.length <- sqrt(sum((par("pin")/c(length(gx), length(gy)))^2)) 0.25 cat("max.grid.arrow.length=", max.grid.arrow.length, " ") } garrows <- do.call(rbind, lapply(gx, function(x) { cd <- sqrt(outer(emb[, 2], -gy, "+")^2 + (x - emb[, 1])^2) cw <- dnorm(cd, sd = grid.sd) gw <- Matrix::colSums(cw) cws <- pmax(1, Matrix::colSums(cw)) gxd <- Matrix::colSums(cw arsd$xd)/cws gyd <- Matrix::colSums(cw arsd$yd)/cws al <- sqrt(gxd^2 + gyd^2) vg <- gw >= min.grid.cell.mass & al >= min.arrow.size cbind(rep(x, sum(vg)), gy[vg], x + gxd[vg], gy[vg] + gyd[vg]) })) colnames(garrows) <- c("x0", "y0", "x1", "y1") if (fixed.arrow.length) { suppressWarnings(arrows(garrows[, 1], garrows[, 2], garrows[, 3], garrows[, 4], length = 0.05, lwd = arrow.lwd)) } else { alen <- pmin(max.grid.arrow.length, sqrt(((garrows[, 3] - garrows[, 1]) par("pin")[1]/diff(par("usr")[c(1, 2)]))^2 + ((garrows[, 4] - garrows[, 2]) par("pin")[2]/diff(par("usr")[c(3, 4)]))^2)) suppressWarnings(lapply(1:nrow(garrows), function(i) arrows(garrows[i, 1], garrows[i, 2], garrows[i, 3], garrows[i, 4], length = alen[i], lwd = arrow.lwd))) } if (plot.grid.points) points(rep(gx, each = length(gy)), rep(gy, length(gx)), pch = ".", cex = 0.1, col = ac(1, alpha = 0.4)) cat("done\n") if (return.details) { cat("expression shifts .") scale.int <- switch(scale, log = 2, sqrt = 3, 1) if (!expression.scaling) { tpb <- tp > 0 tpb <- t(t(tpb)/colSums(tpb)) es <- as.matrix(em %% tp) - as.matrix(em %% as.matrix(tpb)) } cat(".") gs <- do.call(cbind, parallel::mclapply(gx, function(x) { cd <- sqrt(outer(emb[, 2], -gy, "+")^2 + (x - emb[, 1])^2) cw <- dnorm(cd, sd = grid.sd) gw <- Matrix::colSums(cw) cws <- pmax(1, Matrix::colSums(cw)) cw <- t(t(cw)/cws) gxd <- Matrix::colSums(cw arsd$xd) gyd <- Matrix::colSums(cw arsd$yd) al <- sqrt(gxd^2 + gyd^2) vg <- gw >= min.grid.cell.mass & al >= min.arrow.size if (any(vg)) { z <- es %% cw[, vg] } else { NULL } }, mc.cores = n.cores, mc.preschedule = T)) if (scale == "log") { nd <- (log10(abs(nd) + 1) sign(nd)) } else if (scale == "sqrt") { nd <- (sqrt(abs(nd)) sign(nd)) } cat(".") gv <- do.call(cbind, parallel::mclapply(gx, function(x) { cd <- sqrt(outer(emb[, 2], -gy, "+")^2 + (x - emb[, 1])^2) cw <- dnorm(cd, sd = grid.sd) gw <- Matrix::colSums(cw) cws <- pmax(1, Matrix::colSums(cw)) cw <- t(t(cw)/cws) gxd <- Matrix::colSums(cw arsd$xd) gyd <- Matrix::colSums(cw arsd$yd) al <- sqrt(gxd^2 + gyd^2) vg <- gw >= min.grid.cell.mass & al >= min.arrow.size if (any(vg)) { z <- nd %% cw[, vg] } else { NULL } }, mc.cores = n.cores, mc.preschedule = T)) cat(". done\n") return(invisible(list(tp = tp, cc = cc, garrows = garrows, arrows = as.matrix(ars), vel = nd, eshifts = es, gvel = gv, geshifts = gs, scale = scale))) } } else { apply(ars, 1, function(x) { if (fixed.arrow.length) { suppressWarnings(arrows(x[1], x[2], x[3], x[4], length = 0.05, lwd = arrow.lwd)) } else { ali <- sqrt(((x[3] - x[1]) par("pin")[1]/diff(par("usr")[c(1, 2)]))^2 + ((x[4] - x[2]) * par("pin")[2]/diff(par("usr")[c(3, 4)]))^2) suppressWarnings(arrows(x[1], x[2], x[3], x[4], length = min(0.05, ali), lwd = arrow.lwd)) } }) } } return(invisible(list(tp = tp, cc = cc))) }
Hello,
I have got the following error while running show.velocity.on.embedding.cor function.
Any idea what's wrong here and how to correct it?
Many thanks