determine.timing: Error in peak.rle[peak, "start"]:peak.rle[peak, "end"] : NA/NaN argument

[grimwoo]

Hi，I am a freshman in URD. And thank you for the development of such a nice package!

Here, I want to plot a figure like Figure 2B in "Molecular logic of cellular diversification in the mouse cerebral cortex":

And I found a tutorial to plot a figure like this in Cell.

However, the code resulted in error after the "geneSmoothFit" function.

spline <- geneSmoothFit(axial, method = "spline", pseudotime = "pseudotime", cells = cellsInCluster(axial, "stage", c("W5", "W6", "E10", "E11")), genes = var.genes, moving.window = 1, cells.per.window = 5, spar = 0.9 )

determine.timing <- function(s, genes = rownames(s$mean.expression)) { s$timing <- as.data.frame(do.call("rbind", lapply(genes, function(g) { sv <- as.numeric(s$scaled.smooth[g, ]) pt <- as.numeric(colnames(s$scaled.smooth)) min.val <- max(min(sv), 0) peak.val <- ((1 - min.val)/2) + min.val exp.val <- ((1 - min.val)/5) + min.val peak.rle <- rle(sv >= peak.val) peak.rle <- data.frame(lengths = peak.rle$lengths, values = peak.rle$values) peak.rle$end <- cumsum(peak.rle$lengths) peak.rle$start <- head(c(0, peak.rle$end) + 1, -1) exp.rle <- rle(sv >= exp.val) exp.rle <- data.frame(lengths = exp.rle$lengths, values = exp.rle$values) exp.rle$end <- cumsum(exp.rle$lengths) exp.rle$start <- head(c(0, exp.rle$end) + 1, -1) peak <- which(peak.rle$values) peak <- peak[order(peak.rle[peak, "lengths"], decreasing = T)][1:min(2, length(peak))] peak <- peak[order(peak.rle[peak, "start"], decreasing = T)][1] peak <- which.max(sv[peak.rle[peak, "start"]:peak.rle[peak, "end"]]) + peak.rle[peak, "start"] - 1 exp.start <- exp.rle[which(exp.rle$end >= peak & exp.rle$start <= peak), "start"] exp.end <- exp.rle[which(exp.rle$end >= peak & exp.rle$start <= peak), "end"] smooth.start <- sv[exp.start] smooth.end <- sv[exp.end] exp.start <- pt[exp.start] exp.end <- pt[exp.end] peak <- pt[peak] v <- c(exp.start, peak, exp.end, smooth.start, smooth.end) names(v) <- c("pt.start", "pt.peak", "pt.end", "exp.start", "exp.end") return(v) }))) rownames(s$timing) <- genes s$gene.order <- rownames(s$timing)[order(s$timing$pt.peak, s$timing$pt.start, s$timing$pt.end, s$timing$exp.end, decreasing = c(F, F, F, T), method = "radix")] return(s) }

spline <- determine.timing(s = spline)

Error here

Error in peak.rle[peak, "start"]:peak.rle[peak, "end"] : NA/NaN argument

The following code is : filter.heatmap.genes <- function(genes) { mt.genes <- grep("^mt-", ignore.case = T, genes, value = T) many.genes <- grep("\(1 of many\)", ignore.case = T, genes, value = T) ribo.genes <- grep("^rpl|^rps", ignore.case = T, genes, value = T) cox.genes <- grep("^cox", ignore.case = T, genes, value = T) return(setdiff(genes, c(mt.genes, many.genes, ribo.genes, cox.genes))) } order <- filter.heatmap.genes(spline$gene.order)

spline$scaled.smooth[spline$scaled.smooth < 0] <- 0

gplots::heatmap.2(x = as.matrix(spline$scaled.smooth[order, ]), Rowv = F, Colv = F, dendrogram = "none", trace = "none", density.info = "none", key = F, cexCol = 0.8, cexRow = 0.15, margins = c(8, 8), lwid = c(0.3, 4), lhei = c(0.3,

### Could any help me to fix this error or teach me how to get a similar heatmap with gene ordered according to the pseudotime?