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?