diff --git a/R/fitfunc.R b/R/fitfunc.R index 1d3446a..433b9cb 100644 --- a/R/fitfunc.R +++ b/R/fitfunc.R @@ -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){ @@ -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)) @@ -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)) @@ -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 @@ -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)) diff --git a/R/fitfunc_h5.R b/R/fitfunc_h5.R index 94510b8..b699ca9 100644 --- a/R/fitfunc_h5.R +++ b/R/fitfunc_h5.R @@ -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){ @@ -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)) @@ -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)) @@ -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 @@ -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)) diff --git a/R/fitpt.R b/R/fitpt.R index e9b1118..77bfae8 100644 --- a/R/fitpt.R +++ b/R/fitpt.R @@ -281,7 +281,4 @@ fitpt <- function(expr, cellanno, pseudotime, design, testvar=testvar,maxknotall } } list(parameter=para[rownames(expr)],knotnum=knotnum[rownames(expr)]) -} - - - +} \ No newline at end of file diff --git a/R/fitpt.m0.R b/R/fitpt.m0.R index 8ee8f73..fcf04a0 100644 --- a/R/fitpt.m0.R +++ b/R/fitpt.m0.R @@ -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) diff --git a/R/fitpt_h5.R b/R/fitpt_h5.R index 86bd186..d8f132c 100644 --- a/R/fitpt_h5.R +++ b/R/fitpt_h5.R @@ -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 diff --git a/R/fitpt_m0_h5.R b/R/fitpt_m0_h5.R index 22c4168..687a2a1 100644 --- a/R/fitpt_m0_h5.R +++ b/R/fitpt_m0_h5.R @@ -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) diff --git a/R/getCovariateGroupDiff.R b/R/getCovariateGroupDiff.R index 2c974a2..811567e 100644 --- a/R/getCovariateGroupDiff.R +++ b/R/getCovariateGroupDiff.R @@ -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) diff --git a/R/getPopulationFit.R b/R/getPopulationFit.R index da10f4d..ec6da46 100644 --- a/R/getPopulationFit.R +++ b/R/getPopulationFit.R @@ -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 { @@ -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') { @@ -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) diff --git a/R/lamian_test.R b/R/lamian_test.R index e22bc50..4bffad3 100644 --- a/R/lamian_test.R +++ b/R/lamian_test.R @@ -27,7 +27,7 @@ #' data(mandata) #' a = lamian_test(expr = mandata$expr[seq(1,3),], cellanno = mandata$cellanno, pseudotime = mandata$pseudotime, design = mandata$design, test.method = 'chisq', permuiter = 10, EMmaxiter = 10, EMitercutoff = 10, verbose.output = FALSE) -lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, permuiter=100, EMmaxiter=100, EMitercutoff=0.05, verbose.output = FALSE, ncores=detectCores(), test.type='Time', fit.resolution = 1000, return.all.data = TRUE, overall.only = FALSE, test.method = 'permutation', ncores.fit = 1, fix.all.zero = TRUE, cutoff = 1e-5, sd.adjust = 1) { +lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, permuiter=100, maxknotallowed = 10, EMmaxiter=100, EMitercutoff=0.05, verbose.output = FALSE, ncores=detectCores(), test.type='Time', fit.resolution = 1000, return.all.data = TRUE, overall.only = FALSE, test.method = 'permutation', ncores.fit = 1, fix.all.zero = TRUE, cutoff = 1e-5, sd.adjust = 1e-5) { if (test.method == 'permutation') ncores.fit = 1 set.seed(12345) cellanno = data.frame(Cell = as.character(cellanno[,1]), Sample = as.character(cellanno[,2]), stringsAsFactors = FALSE) @@ -60,6 +60,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm pseudotime, design = design[, 1, drop = FALSE], testvar = testvar, + maxknotallowed = maxknotallowed, maxknotallowed = 10, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, @@ -92,7 +93,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm pseudotime, design, testvar = testvar, - maxknotallowed = 10, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, ncores = ncores.fit, @@ -107,7 +108,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm pseudotime, design, testvar = testvar, - maxknotallowed = 10, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, ncores = ncores.fit, @@ -123,7 +124,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm pseudotime, design, testvar = testvar, - maxknotallowed = 10, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, ncores = ncores.fit, @@ -165,6 +166,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'overall', test.type = test.type, testvar = testvar, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -182,6 +184,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'overall', test.type = test.type, testvar = testvar, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -219,14 +222,15 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm llr.overall <- ll.full - ll.null pval.overall <- sapply(seq_len(nrow(llr.overall)), function(i) { z <- llr.overall[i, seq(2, ncol(llr.overall))] + z <- z[!is.na(z)] den <- density(z)$bw - mean(pnorm(llr.overall[i,1], z, sd=den, lower.tail = F)) + mean(pnorm(llr.overall[i, 1], z, sd = den, lower.tail = F)) }) - log.pval <- sapply(seq_len(nrow(llr.overall)), function(i) { z <- llr.overall[i, seq(2, ncol(llr.overall))] + z <- z[!is.na(z)] den <- density(z)$bw - max(pnorm(llr.overall[i,1], z, sd=den, lower.tail = F, log.p = T)) + max(pnorm(llr.overall[i, 1], z, sd = den, lower.tail = F, log.p = T)) }) fdr.overall <- p.adjust(pval.overall,method='fdr') names(pval.overall) <- names(fdr.overall) <- row.names(llr.overall) @@ -245,6 +249,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall < 0.05], test.type = test.type, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -263,6 +268,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall < 0.05], test.type = test.type, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -281,6 +287,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm llr <- llr[complete.cases(llr), ] if (sum(fdr.overall<0.05) == 1){ z <- llr[seq(2, length(llr))] + z <- z[!is.na(z)] den <- density(z)$bw fdr <- pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F)) log.pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F,log.p=T)) @@ -289,11 +296,13 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm } else { pval <- sapply(seq_len(nrow(llr)), function(i) { z <- llr[i, seq(2, ncol(llr))] + z <- z[!is.na(z)] den <- density(z)$bw mean(pnorm(llr[i,1], z, sd=den,lower.tail = F)) }) log.pval <- sapply(seq_len(nrow(llr)), function(i) { z <- llr[i, seq(2, ncol(llr))] + z <- z[!is.na(z)] den <- density(z)$bw max(pnorm(llr[i,1], z, sd=den,lower.tail = F, log.p = T)) }) @@ -311,6 +320,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall < 0.05], test.type = test.type, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -329,6 +339,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall < 0.05], test.type = test.type, + maxknotallowed = maxknotallowed, EMmaxiter = EMmaxiter, EMitercutoff = EMitercutoff, verbose.output = verbose.output, @@ -347,6 +358,7 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm llr <- llr[complete.cases(llr), ] if (sum(fdr.overall<0.05) == 1){ z <- llr[seq(2, length(llr))] + z <- z[!is.na(z)] den <- density(z)$bw fdr <- pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F)) log.pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F,log.p=T)) @@ -355,11 +367,13 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm } else { pval <- sapply(seq_len(nrow(llr)), function(i) { z <- llr[i, seq(2, ncol(llr))] + z <- z[!is.na(z)] den <- density(z)$bw mean(pnorm(llr[i,1], z, sd=den,lower.tail = F)) }) log.pval <- sapply(seq_len(nrow(llr)), function(i) { z <- llr[i, seq(2, ncol(llr))] + z <- z[!is.na(z)] den <- density(z)$bw max(pnorm(llr[i,1], z, sd=den,lower.tail = F, log.p = T)) }) @@ -390,4 +404,3 @@ lamian_test <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, perm return(c(reslist, list(test.type = test.type, test.method = test.method, testvar = testvar))) } } - diff --git a/R/lamian_test_h5.R b/R/lamian_test_h5.R index 08dc4a9..4a66b88 100644 --- a/R/lamian_test_h5.R +++ b/R/lamian_test_h5.R @@ -27,7 +27,7 @@ #' data(expdata) #' res <- lamian_test_h5(expr='data/multi.h5', cellanno=expdata$cellanno, pseudotime=expdata$pseudotime, design=expdata$design, testvar=2, test.type = 'Variable', overall.only = F, test.method = 'chisq') -lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, permuiter=100, EMmaxiter=100, EMitercutoff=0.05, verbose.output=F, ncores=detectCores(), test.type='Time', fit.resolution = 1000, return.all.data = TRUE, overall.only = F, test.method = 'permutation', ncores.fit = 1, fix.all.zero = TRUE, cutoff = 1e-5) { +lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, permuiter=100, maxknotallowed = 10, EMmaxiter=100, EMitercutoff=0.05, verbose.output=F, ncores=detectCores(), test.type='Time', fit.resolution = 1000, return.all.data = TRUE, overall.only = F, test.method = 'permutation', ncores.fit = 1, fix.all.zero = TRUE, cutoff = 1e-5) { cellanno <- cellanno[match(names(pseudotime), cellanno[,1]), ] if (test.method == 'permutation') ncores.fit = 1 set.seed(12345) @@ -38,9 +38,9 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p if (test.method == 'chisq'){ if (toupper(test.type) == 'TIME'){ - res1 <- fitpt_h5(expr, pseudotime, design=design[,1,drop=FALSE], maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 1,testvar=testvar)##save 13% + res1 <- fitpt_h5(expr, pseudotime, design=design[,1,drop=FALSE], maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 1,testvar=testvar)##save 13% ll1 <- sapply(res1$parameter,function(i) i$ll) - res0 <- fitpt_h5(expr, pseudotime, design[,1,drop=FALSE],EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff) ## + res0 <- fitpt_h5(expr, pseudotime, design[,1,drop=FALSE], maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff) ## ll0 <- sapply(res0[[1]], function(i) i$ll) paradiff10 <- sapply(res1[[1]], function(i) length(unlist(i[1:4]))) - sapply(res0[[1]], function(i) length(unlist(i[1:4]))) pval.chisq.constantTest <- pchisq(2*(ll1-ll0),df=paradiff10,lower.tail = F) @@ -52,11 +52,11 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p stringsAsFactors = FALSE) reslist = list(statistics = res, parameter = res1$parameter, knotnum = res1$knotnum) ## function return } else if (toupper(test.type) == 'VARIABLE'){ - res1 <- fitpt_h5(expr, pseudotime, design, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 1,testvar=testvar)## save 13% + res1 <- fitpt_h5(expr, pseudotime, design, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 1,testvar=testvar)## save 13% ll1 <- sapply(res1$parameter,function(i) i$ll) - res2 <- fitpt_h5(expr, pseudotime, design, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 2, knotnum = res1[[2]],testvar=testvar)## save 13% + res2 <- fitpt_h5(expr, pseudotime, design, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 2, knotnum = res1[[2]],testvar=testvar)## save 13% ll2 <- sapply(res2$parameter,function(i) i$ll) - res3 <- fitpt_h5(expr, pseudotime, design, maxknotallowed=10, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 3, knotnum = res1[[2]],testvar=testvar) + res3 <- fitpt_h5(expr, pseudotime, design, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, ncores=ncores.fit, model = 3, knotnum = res1[[2]],testvar=testvar) ll3 <- sapply(res3$parameter,function(i) i$ll) paradiff31 <- sapply(res3$parameter,function(i) length(i$beta))-sapply(res1$parameter,function(i) length(i$beta)) @@ -84,8 +84,8 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p } else if (test.method == 'permutation'){ if (verbose.output) print('fitting model: overall: CovariateTest (Model 3 vs.1) or ConstantTest (Model 1) ...') if (ncores == 1){ - fit <- lapply(1:(permuiter+1),function(i) fitfunc_h5(iter = i, diffType = 'overall', test.type = test.type, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar))} else { - fit <- mclapply(1:(permuiter+1),function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'overall', test.type = test.type, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) + fit <- lapply(1:(permuiter+1),function(i) fitfunc_h5(iter = i, diffType = 'overall', test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar))} else { + fit <- mclapply(1:(permuiter+1),function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'overall', test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) } if (verbose.output) { print('The length of fit is ...') ## @@ -108,22 +108,23 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p } } - # saveRDS(fit, '/home-4/whou10@jhu.edu/scratch/Wenpin/trajectory_variability/debug/fitoverall.rds') + # saveRDS(fit, '/projectnb/wax-es/Bingtian/Lamian/test_for_design/4sample_TCPO/fitoverall.RDS') knotnum <- fit[[1]]$fitres.full$knotnum parameter <- fit[[1]]$fitres.full$parameter ll.full <- sapply(1:(length(fit)),function(i) sapply(fit[[i]]$fitres.full$parameter,function(j) unname(j$ll),USE.NAMES = F)[gn]) ll.null <- sapply(1:(length(fit)),function(i) sapply(fit[[i]]$fitres.null$parameter,function(j) unname(j$ll),USE.NAMES = F)[gn]) llr.overall <- ll.full - ll.null - pval.overall <- sapply(1:nrow(llr.overall), function(i) { - z <- llr.overall[i,2:ncol(llr.overall)] + pval.overall <- sapply(seq_len(nrow(llr.overall)), function(i) { + z <- llr.overall[i, seq(2, ncol(llr.overall))] + z <- z[!is.na(z)] den <- density(z)$bw - mean(pnorm(llr.overall[i,1], z, sd=den, lower.tail = F)) + mean(pnorm(llr.overall[i, 1], z, sd = den, lower.tail = F)) }) - - log.pval <- sapply(1:nrow(llr.overall), function(i) { - z <- llr.overall[i,2:ncol(llr.overall)] + log.pval <- sapply(seq_len(nrow(llr.overall)), function(i) { + z <- llr.overall[i, seq(2, ncol(llr.overall))] + z <- z[!is.na(z)] den <- density(z)$bw - max(pnorm(llr.overall[i,1], z, sd=den, lower.tail = F, log.p = T)) + max(pnorm(llr.overall[i, 1], z, sd = den, lower.tail = F, log.p = T)) }) fdr.overall <- p.adjust(pval.overall,method='fdr') names(pval.overall) <- names(fdr.overall) <- row.names(llr.overall) @@ -135,9 +136,9 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p if (sum(fdr.overall <0.05 ) > 0 & toupper(test.type) == 'VARIABLE' & !overall.only){ if (verbose.output) print('meanDiff pvalues: Model 2 vs. model 1...') if (ncores == 1){ - fit <- lapply(seq_len(permuiter+1),function(i) fitfunc_h5(iter = i, diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)) + fit <- lapply(seq_len(permuiter+1),function(i) fitfunc_h5(iter = i, diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)) } else { - fit <- mclapply(seq_len(permuiter+1),function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type,EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) + fit <- mclapply(seq_len(permuiter+1),function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'meanDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) } fit <- fit[!sapply(fit,is.null)] ll.full <- sapply(1:(length(fit)),function(i) sapply(fit[[i]]$fitres.full$parameter,function(j) unname(j$ll),USE.NAMES = F)[gn]) @@ -146,6 +147,7 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p llr <- llr[complete.cases(llr), ] if (sum(fdr.overall<0.05) == 1){ z <- llr[2:length(llr)] + z <- z[!is.na(z)] den <- density(z)$bw fdr <- pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F)) log.pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F,log.p=T)) @@ -154,11 +156,13 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p } else { pval <- sapply(1:nrow(llr), function(i) { z <- llr[i,2:ncol(llr)] + z <- z[!is.na(z)] den <- density(z)$bw mean(pnorm(llr[i,1], z, sd=den,lower.tail = F)) }) log.pval <- sapply(1:nrow(llr), function(i) { z <- llr[i,2:ncol(llr)] + z <- z[!is.na(z)] den <- density(z)$bw max(pnorm(llr[i,1], z, sd=den,lower.tail = F, log.p = T)) }) @@ -170,9 +174,9 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p if (verbose.output) print('trendDiff pvalues: Model 3 vs. model 2...') if (ncores == 1){ - fit <- lapply(1:(permuiter+1), function(i) fitfunc_h5(iter = i, diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)) + fit <- lapply(1:(permuiter+1), function(i) fitfunc_h5(iter = i, diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)) } else { - fit <- mclapply(1:(permuiter+1), function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) ## return a list of (permuiter + 1) where the first is a list of fitres.full and fitres.null + fit <- mclapply(1:(permuiter+1), function(i){set.seed(i); fitfunc_h5(iter = i, diffType = 'trendDiff', gene = names(fdr.overall)[fdr.overall<0.05], test.type = test.type, maxknotallowed = maxknotallowed, EMmaxiter=EMmaxiter, EMitercutoff=EMitercutoff, verbose.output=verbose.output, expr=expr, cellanno=cellanno, pseudotime=pseudotime, design=design,testvar=testvar)}, mc.cores = ncores) ## return a list of (permuiter + 1) where the first is a list of fitres.full and fitres.null } fit <- fit[!sapply(fit,is.null)] ll.full <- sapply(1:(length(fit)),function(i) sapply(fit[[i]]$fitres.full$parameter,function(j) unname(j$ll),USE.NAMES = F)[gn]) @@ -181,6 +185,7 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p llr <- llr[complete.cases(llr), ] if (sum(fdr.overall<0.05) == 1){ z <- llr[2:length(llr)] + z <- z[!is.na(z)] den <- density(z)$bw fdr <- pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F)) log.pval <- mean(pnorm(llr[1], z, sd=den,lower.tail = F,log.p=T)) @@ -189,11 +194,13 @@ lamian_test_h5 <- function(expr, cellanno, pseudotime, design=NULL, testvar=2, p } else { pval <- sapply(1:nrow(llr), function(i) { z <- llr[i,2:ncol(llr)] + z <- z[!is.na(z)] den <- density(z)$bw mean(pnorm(llr[i,1], z, sd=den,lower.tail = F)) }) log.pval <- sapply(1:nrow(llr), function(i) { z <- llr[i,2:ncol(llr)] + z <- z[!is.na(z)] den <- density(z)$bw max(pnorm(llr[i,1], z, sd=den,lower.tail = F, log.p = T)) }) diff --git a/R/plotGenePopulation.R b/R/plotGenePopulation.R index 9a3bc8c..96c2594 100644 --- a/R/plotGenePopulation.R +++ b/R/plotGenePopulation.R @@ -89,9 +89,10 @@ plotGenePopulation <- function(testobj, if (is.na(gene)) gene <- rownames(fit[[1]]) pd <- sapply(seq_len(length(fit)), function(i) { - tmp <- reshape2::melt(fit[[i]][gene, , drop = FALSE]) + pt_i <- seq(1, max(testobj$pseudotime), length.out = ncol(fit[[i]])) + colnames(fit[[i]]) <- as.character(pt_i) + if (!is.na(subSampleNumber)) { - ## set.seed(12345) id <- sample(seq_len(nrow(fit[[i]])), subSampleNumber) tmp <- reshape2::melt(fit[[i]][gene, id , drop = FALSE]) } else { @@ -99,10 +100,10 @@ plotGenePopulation <- function(testobj, } colnames(tmp) <- c('gene', 'pseudotime', 'expression') + tmp$pseudotime <- as.numeric(as.character(tmp$pseudotime)) + if (!is.na(subSampleNumber)) { - ## set.seed(12345) - tmp <- - tmp[sample(seq_len(nrow(tmp)), subSampleNumber), , drop = FALSE] + tmp <- tmp[sample(seq_len(nrow(tmp)), subSampleNumber), , drop = FALSE] } tmp <- data.frame(tmp, type = names(fit)[i], @@ -132,6 +133,7 @@ plotGenePopulation <- function(testobj, xlab('Pseudotime') + ylab('Expression') + labs(color = '') + if(length(gene) != 1){ if (is.na(ncol)) { p <- p + facet_wrap( ~ gene, nrow = nrow, scales = a) } else { @@ -142,7 +144,7 @@ plotGenePopulation <- function(testobj, } else { p <- p + scale_color_manual(values = grDevices::colorRampPalette(brewer.pal(8, palette))(length(unique(pd$type)))) - } + }} } if (!is.na(ylim)[1]) @@ -156,5 +158,4 @@ plotGenePopulation <- function(testobj, p <- p + scale_x_continuous(breaks = c(min(pd$pseudotime), max(pd$pseudotime))) } - print(p) } diff --git a/R/plotSingleGene.R b/R/plotSingleGene.R new file mode 100644 index 0000000..988929d --- /dev/null +++ b/R/plotSingleGene.R @@ -0,0 +1,63 @@ +#' Plot gene(s) by showing both of the original cellcular gene expression and the sample-level, population-level fitting values. +#' +#' This function is used to plot gene(s) by showing both of the original cellcular gene expression and the sample-level, population-level fitting values. +#' +#' @import ggplot2 RColorBrewer splines gridExtra viridis +#' @importFrom grDevices colorRampPalette +#' @return a plot +#' @author Wenpin Hou +#' @export +#' @param testobj object returned from lamian_test(). +#' @param gene a character vector of gene names. It can be of length 1 or > 1. +#' @param type One of c('Time', 'Variable'). +#' @param variable character, the variable (covariate) to color the samples, should be null or one of the column names of design matrix. Default is NULL, meaning each sample is colored differently. Otherwise, samples are colored by the variable (covariate) values. +#' @param variable.text a character vector. The text for the legend of the plot, corresponding to each variable values. +#' @param facet.sample logical. If TRUE (default), facet_wrap the samples. +#' @param plot.point point size +#' @param line.alpha alpha value of the curves +#' @param continuous if TRUE, samples are colored using viridis continuous colors. If FALSE, RColorBrewer "Dark2" discrete palette. +#' @param cellProp logical. If FALSE (default), plot gene expression. If TRUE, it is cell proportion. +#' @param x.lab a string to indicates x-axis label +#' @param y.lab a string to indicates y-axis label +#' @param point.size the size value of points. +#' @param point.alpha the alpha value of points. +#' @param ylim y-axis limits to be passed to ggplot. +#' @param xlim x-axis limits to be passed to ggplot. +#' @param sep a string in the gene names that needs to replaced with blank. +#' @param free.scale logical. If TRUE, the y-axis is on free scale. +#' @param palette a RColorBrewer palette name. +#' @param ncol number of colums for organizing all genes' plot. +#' @param line.size the size of the curves. +#' @param axis.text.blank logical. If TRUE, leave axis text as blank. +#' @examples +#' data(mantestobj) +#' plotGene(testobj = mantestobj, gene = rownames(mantestobj$populationFit[[1]])[1], variable = 'gender') +plotSingleGene <- + function(testobj, + gene, + variable = NULL, + variable.text = NULL, + free.scale = TRUE, + facet.sample = FALSE, + plot.point = FALSE, + line.alpha = 1, + line.size = 1, + point.alpha = 1, + point.size = 0.5, + continuous = TRUE, + sep = NA, + palette = 'Dark2', + ncol = NULL, + axis.text.blank = FALSE, + x.lab = 'Pseudotime', + y.lab = 'Expression') { + pseudotime <- testobj[['pseudotime']] + cellanno <- testobj[['cellanno']] + colnames(cellanno) <- c('Cell', 'Sample') + expression <- testobj[['expr']][gene, , drop = FALSE] + + predict.values <- + predict_fitting(testobj, gene = gene, test.type = testobj$test.type) + + + }