Skip to content

Robust handling of root-finding errors in fitpt.m0: skip genes with uniroot failures and return NA instead of stopping analysis #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Ye1203 wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Robust handling of root-finding errors in fitpt.m0: skip genes with uniroot failures and return NA instead of stopping analysis #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
79a1a25
debug btye
Ye1203 Jun 18, 2025
7532c7c
updated 250624
Ye1203 Jun 24, 2025
9e55698
updated 250624
Ye1203 Jun 24, 2025
9c08288
fix not full rank perdesign
Ye1203 Jun 25, 2025
ce1b11c
updated 250702
Ye1203 Jul 2, 2025
dd4c294
change default sd.adjust
Ye1203 Jul 18, 2025
f7c7afb
updated
Ye1203 Sep 22, 2025
9e52a34
updated
Ye1203 Sep 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 18 additions & 12 deletions R/fitfunc.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,13 @@
#' @param ncores the number of cores to be used. If ncores > 1, it will be implemented in parallel mode.
#' @param verbose.output logical. If TRUE, print intermediate information.

fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type = 'Time', testvar=testvar, EMmaxiter=100, EMitercutoff=0.05, ncores=1, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design, verbose.output = FALSE) {
fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type = 'Time', testvar=testvar, maxknotallowed = 10, EMmaxiter=100, EMitercutoff=0.05, ncores=1, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design, verbose.output = FALSE) {
expr <- expr[gene, ,drop=FALSE]
if (verbose.output) print(paste0('iter ', iter, '\n'))
if (toupper(test.type)=='TIME') {
if (iter == 1){
fitres.full <- fitpt(expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design[,1,drop=FALSE], testvar=testvar,EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model=-1)
fitres.null <- fitpt.m0(expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design[,1,drop=FALSE], EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff)
fitres.full <- fitpt(expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design[,1,drop=FALSE], testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model=-1)
fitres.null <- fitpt.m0(expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design[,1,drop=FALSE], EMmaxiter=EMmaxiter, maxknotallowed = maxknotallowed, EMitercutoff=EMitercutoff)
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
} else {
perpsn <- lapply(rownames(design), function(s){
Expand All @@ -45,8 +45,8 @@ fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type
percellanno <- cellanno[sampcell,,drop=F]
perpsn <- perpsn[sampcell]
colnames(perexpr) <- percellanno[,1] <- names(perpsn) <- paste0('cell_',1:length(perpsn))
tryCatch(fitres.full <- fitpt(expr=perexpr, cellanno=percellanno, pseudotime=perpsn, design=design,testvar=testvar, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = -1), warning = function(w){}, error = function(e) {})
tryCatch(fitres.null <- fitpt.m0(expr=perexpr, cellanno=percellanno, pseudotime=perpsn, design=design, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff), warning = function(w){}, error = function(e) {})
tryCatch(fitres.full <- fitpt(expr=perexpr, cellanno=percellanno, pseudotime=perpsn, design=design,testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = -1), warning = function(w){}, error = function(e) {})
tryCatch(fitres.null <- fitpt.m0(expr=perexpr, cellanno=percellanno, pseudotime=perpsn, design=design, EMmaxiter=EMmaxiter, maxknotallowed = maxknotallowed, EMitercutoff=EMitercutoff), warning = function(w){}, error = function(e) {})
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand All @@ -68,8 +68,8 @@ fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type
mod.null = 2
}
if (iter == 1){
fitres.full <- fitpt(expr, cellanno, pseudotime, design,testvar=testvar, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.full)
fitres.null <- fitpt(expr, cellanno, pseudotime, design, testvar=testvar, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.null, knotnum = fitres.full[[2]])
fitres.full <- fitpt(expr, cellanno, pseudotime, design, testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.full)
fitres.null <- fitpt(expr, cellanno, pseudotime, design, testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.null, knotnum = fitres.full[[2]])
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand All @@ -80,11 +80,17 @@ fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type
} else {
dn <- paste0(design[,testvar],collapse = '_')
perdn <- dn
while(perdn==dn) {
# make sure perdesign is full rank
while (TRUE) {
perid <- sample(1:nrow(design))
perdesign <- design
perdesign[,testvar] <- design[perid,testvar]
perdn <- paste0(perdesign[,testvar],collapse = '_')
perdesign[, testvar] <- design[perid, testvar]
rnk <- as.integer(Matrix::rankMatrix(perdesign))
perdn_new <- paste0(perdesign[, testvar], collapse = '_')
if (perdn_new != dn && rnk == ncol(perdesign)) {
perdn <- perdn_new
break
}
}
row.names(perdesign) <- row.names(design)
sampcell <- sample(1:ncol(expr),replace=T) ## boostrap cells
Expand All @@ -95,8 +101,8 @@ fitfunc <- function(iter, diffType = 'overall', gene = rownames(expr), test.type
colnames(perexpr) <- percellanno[,1] <- names(psn) <- paste0('cell_',1:length(psn))


fitres.full <- fitpt(perexpr, percellanno, psn, perdesign,testvar=testvar, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.full)
fitres.null <- fitpt(perexpr, percellanno, psn, perdesign,testvar=testvar, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.null, knotnum = fitres.full[[2]])
fitres.full <- fitpt(perexpr, percellanno, psn, perdesign,testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.full)
fitres.null <- fitpt(perexpr, percellanno, psn, perdesign,testvar=testvar, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.null, knotnum = fitres.full[[2]])
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand Down
31 changes: 19 additions & 12 deletions R/fitfunc_h5.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,12 @@
#' @param EMitercutoff a numeric number indicating the log-likelihood cutoff applied to stop the EM algorithm
#' @param verbose.output logical. If TRUE, print intermediate information.
#' @param ncores the number of cores to be used. If ncores > 1, it will be implemented in parallel mode.
fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testvar, test.type = 'Time', expr = expr, cellanno = cellanno, pseudotime = pseudotime, design = design, EMmaxiter = 100, EMitercutoff = 0.05, verbose.output = F, ncores = 1) {
fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testvar, maxknotallowed = 10, test.type = 'Time', expr = expr, cellanno = cellanno, pseudotime = pseudotime, design = design, EMmaxiter = 100, EMitercutoff = 0.05, verbose.output = F, ncores = 1) {
if (verbose.output) print(paste0('iter ', iter, '\n'))
if (toupper(test.type)=='TIME') {
if (iter == 1) {
fitres.full <- fitpt_h5(expr=expr, pseudotime=pseudotime, design=design[,1,drop=FALSE],testvar=testvar,targetgene=gene, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model=-1)
fitres.null <- fitpt_m0_h5(expr=expr, pseudotime=pseudotime, design=design[,1,drop=FALSE],targetgene=gene, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff)
fitres.full <- fitpt_h5(expr=expr, pseudotime=pseudotime, design=design[,1,drop=FALSE],testvar=testvar, maxknotallowed = maxknotallowed, targetgene=gene, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model=-1)
fitres.null <- fitpt_m0_h5(expr=expr, pseudotime=pseudotime, design=design[,1,drop=FALSE], targetgene=gene, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff)
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
} else {
perpsn <- lapply(rownames(design), function(s){
Expand All @@ -42,8 +42,8 @@ fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testva
perpsn <- perpsn[sampcell]
boot <- data.frame(percellanno[,1],paste0('cell_',1:length(perpsn)),stringsAsFactors = F) #### rename cells
percellanno[,1] <- names(perpsn) <- paste0('cell_',1:length(perpsn)) ## save the original cell name and permuted cell names relation, not used here actually, because hdf5 file already save the cells for each sample seperately
fitres.full <- fitpt_h5(expr=expr, pseudotime=perpsn, design=design, boot=boot,targetgene=gene,EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = -1,testvar=testvar)
fitres.null <- fitpt_m0_h5(expr=expr, pseudotime=perpsn, design=design, boot=boot,targetgene=gene,EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff)
fitres.full <- fitpt_h5(expr=expr, pseudotime=perpsn, design=design, boot=boot,targetgene=gene, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = -1,testvar=testvar)
fitres.null <- fitpt_m0_h5(expr=expr, pseudotime=perpsn, design=design, boot=boot,targetgene=gene, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff)
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand All @@ -65,8 +65,8 @@ fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testva
mod.null = 2
}
if (iter == 1){
fitres.full <- fitpt_h5(expr, pseudotime, design,targetgene=gene, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.full,testvar=testvar)
fitres.null <- fitpt_h5(expr, pseudotime, design,targetgene=gene, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.null, knotnum = fitres.full[[2]],testvar=testvar)
fitres.full <- fitpt_h5(expr, pseudotime, design,targetgene=gene, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.full, testvar=testvar)
fitres.null <- fitpt_h5(expr, pseudotime, design,targetgene=gene, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=1, model = mod.null, knotnum = fitres.full[[2]],testvar=testvar)
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand All @@ -80,10 +80,17 @@ fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testva
} else {
dn <- paste0(as.vector(design),collapse = '_')
perdn <- dn
while(perdn==dn) {
# make sure perdesign is full rank
while (TRUE) {
perid <- sample(1:nrow(design))
perdesign <- design[perid,,drop=F]
perdn <- paste0(as.vector(perdesign),collapse = '_')
perdesign <- design
perdesign[, testvar] <- design[perid, testvar]
rnk <- as.integer(Matrix::rankMatrix(perdesign))
perdn_new <- paste0(perdesign[, testvar], collapse = '_')
if (perdn_new != dn && rnk == ncol(perdesign)) {
perdn <- perdn_new
break
}
}
row.names(perdesign) <- row.names(design)
sampcell <- sample(1:length(pseudotime),replace=T) ## boostrap cells
Expand All @@ -94,8 +101,8 @@ fitfunc_h5 <- function(iter, diffType = 'overall', gene = NULL, testvar = testva
percellanno[,1] <- names(psn) <- paste0('cell_',1:length(psn))


fitres.full <- fitpt_h5(expr, psn, perdesign, boot=boot,targetgene=gene,maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.full,testvar=testvar)
fitres.null <- fitpt_h5(expr, psn, perdesign, boot=boot,targetgene=gene,maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.null, knotnum = fitres.full[[2]],testvar=testvar)
fitres.full <- fitpt_h5(expr, psn, perdesign, boot=boot,targetgene=gene,maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.full,testvar=testvar)
fitres.null <- fitpt_h5(expr, psn, perdesign, boot=boot,targetgene=gene,maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores, model = mod.null, knotnum = fitres.full[[2]],testvar=testvar)
if (exists('fitres.full') & exists('fitres.null')) {
if (verbose.output) print(paste0('iter ', iter, ' success!'))
return(list(fitres.full = fitres.full, fitres.null = fitres.null))
Expand Down
5 changes: 1 addition & 4 deletions R/fitpt.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -281,7 +281,4 @@ fitpt <- function(expr, cellanno, pseudotime, design, testvar=testvar,maxknotall
}
}
list(parameter=para[rownames(expr)],knotnum=knotnum[rownames(expr)])
}



}
22 changes: 18 additions & 4 deletions R/fitpt.m0.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,10 +95,24 @@ fitpt.m0 <- function(expr, cellanno, pseudotime, design, EMmaxiter=100, EMitercu
Jsolve[,s,drop=F] + N[,s,drop=F]*JK[[s]] * JK[[s]] ## [if debug, here]
}), ncol=length(as), dimnames = list(gidr, as)))

eta[gidr] <- sapply(gidr,function(g) {
meanN <- mean(N[g,])
meanA <- mean(A[g,])
uniroot(function(eta) {digamma(eta * meanN)-log(eta)+meanA},c(1e-10,1e10))$root
eta[gidr] <- sapply(gidr, function(g) {
meanN <- mean(N[g, ])
meanA <- mean(A[g, ])
f <- function(eta) digamma(eta * meanN) - log(eta) + meanA
f_lower <- f(1e-10)
f_upper <- f(1e10)
if (f_lower * f_upper > 0) {
result <- tryCatch(
uniroot(f, c(1e-20, 1e20), tol = 1e-12)$root,
error = function(e) {
warning(paste0("No root found for gene: ", g))
NA
}
)
} else {
result <- uniroot(f, c(1e-10, 1e10))$root
}
result
})
alpha[gidr] <- eta[gidr] * rowMeans(N)
para <- list(beta = B, alpha = alpha, eta = eta, omega = omega)
Expand Down
19 changes: 16 additions & 3 deletions R/fitpt_h5.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,11 +252,24 @@ fitpt_h5 <- function(expr,

rN <- sumN/length(as)
rA <- sumA/length(as)
eta[gidr] <- sapply(gidr,function(g) {
eta[gidr] <- sapply(gidr, function(g) {
meanN <- rN[g]
meanA <- rA[g]
# uniroot(function(eta) {digamma(eta * meanN)-log(eta)+meanA},c(1e-10,1e10))$root
optim(eta[g],fn = function(eta) {(digamma(eta * meanN)-log(eta)+meanA)^2},gr = function(eta) {2*(digamma(eta * meanN)-log(eta)+meanA)*(trigamma(eta*meanN)*meanN-1/eta)},lower = 1e-10,method = 'L-BFGS-B')$par
f <- function(eta) digamma(eta * meanN) - log(eta) + meanA
f_lower <- f(1e-10)
f_upper <- f(1e10)
if (f_lower * f_upper > 0) {
result <- tryCatch(
uniroot(f, c(1e-20, 1e20), tol = 1e-12)$root,
error = function(e) {
warning(paste0("No root found for gene: ", g))
NA
}
)
} else {
result <- uniroot(f, c(1e-10, 1e10))$root
}
result
})
alpha[gidr] <- eta[gidr] * rN

Expand Down
22 changes: 18 additions & 4 deletions R/fitpt_m0_h5.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,10 +125,24 @@ fitpt_m0_h5 <- function(expr, pseudotime, design, targetgene=NULL, boot=NULL, EM
Jsolve[,s,drop=F] + N[,s,drop=F]*JK[[s]] * JK[[s]] ## debug here !!!
}), ncol=length(as), dimnames = list(gidr, as)))

eta[gidr] <- sapply(gidr,function(g) {
meanN <- mean(N[g,])
meanA <- mean(A[g,])
uniroot(function(eta) {digamma(eta * meanN)-log(eta)+meanA},c(1e-10,1e10))$root
eta[gidr] <- sapply(gidr, function(g) {
meanN <- mean(N[g, ])
meanA <- mean(A[g, ])
f <- function(eta) digamma(eta * meanN) - log(eta) + meanA
f_lower <- f(1e-10)
f_upper <- f(1e10)
if (f_lower * f_upper > 0) {
result <- tryCatch(
uniroot(f, c(1e-20, 1e20), tol = 1e-12)$root,
error = function(e) {
warning(paste0("No root found for gene: ", g))
NA
}
)
} else {
result <- uniroot(f, c(1e-10, 1e10))$root
}
result
})
alpha[gidr] <- eta[gidr] * rowMeans(N)

Expand Down
2 changes: 1 addition & 1 deletion R/getCovariateGroupDiff.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ getCovariateGroupDiff <- function(testobj,
num.timepoint = 1e3,
testvar = 2) {
knotnum = testobj$knotnum[gene]
pseudotime = seq(1, max(testobj$pseudotime), length.out = min(num.timepoint, max(testobj$pseudotime)))
pseudotime = seq(1, max(testobj$pseudotime), length.out = num.timepoint)
if ('testvar' %in% names(testobj)) testvar = testobj$testvar
beta <- lapply(gene, function(g) {
tmp = matrix(testobj$parameter[[g]]$beta, ncol = knotnum[g]+4)
Expand Down
8 changes: 4 additions & 4 deletions R/getPopulationFit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ getPopulationFit <- function(testobj,
} else {
testvar = testobj$testvar
}

if (type == 'VARIABLE'){
design = testobj$design[, c(1, testobj$testvar)] ## design for multi
} else {
Expand All @@ -33,14 +33,14 @@ getPopulationFit <- function(testobj,
knotnum = testobj$knotnum
pseudotime = testobj$pseudotime
pseudotime = pseudotime[order(pseudotime)]
pt <- round(seq(1, max(pseudotime), length.out = min(num.timepoint, max(pseudotime)))) ## downsample
pt <- seq(min(pseudotime), max(pseudotime), length.out = num.timepoint)

if (sum(design[, 1]) != nrow(design)){
print("The first column of design matrix should be all 1s (intercept)! Using the first column as the variable column ...")
design = cbind(intercept = 1, design)
colnames(design)[1] <- 'intercept'
}

if (is.null(gene)) gene <- rownames(testobj$statistics)

if (type == 'TIME') {
Expand All @@ -58,7 +58,7 @@ getPopulationFit <- function(testobj,
} else {
beta = as.vector(tmp[1, ]) ### subset the beta values of the intercept and the test covariate for multi
}

x <- sapply(row.names(design), function(i) {
kronecker(diag(knotnum[g] + 4), design[i, , drop = FALSE]) ###
}, simplify = FALSE)
Expand Down
Loading