Is the function: get.exp.stat() using p.value or FDR as default
Hi Danko-Lab, thank you for this awesome tool!
I am using the function: get.exp.stat() and i was wondering if it automatically corrects for multiple testing on the p-values for each gene or if I have to change something in the code in order to get the FDR from the function: scran::combineMarkers()
I have tried to read the doc https://rdrr.io/bioc/scran/man/combineMarkers.html and have made some changes to the code marked with " #******* "
The code:::::::::::::::
get.exp.stat.1 <- function(sc.dat, cell.type.labels, cell.state.labels, cell.count.cutoff = 50, n.cores = 1) {
ct.to.cst <- unique(cbind(cell.type = cell.type.labels, cell.state = cell.state.labels)) cst.count.table <- table(cell.state.labels) low.count.cst <- names(cst.count.table)[cst.count.table < cell.count.cutoff]
lib.size <- rowSums(sc.dat) lib.size <- lib.size / median(lib.size) dat.tmp <- sc.dat / lib.size dat.tmp <- log2(dat.tmp + 0.1) - log2(0.1) #******* I change the pseudo.count parm an set it to 0.1 for 10x data
fit.up <- scran::pairwiseTTests(x = t(dat.tmp),
groups = cell.state.labels,
direction = "up",
BPPARAM = BiocParallel::MulticoreParam(n.cores))
pairs.celltype.first <- ct.to.cst[match(fit.up$pairs$first, ct.to.cst[,"cell.state"]),"cell.type"] pairs.celltype.second <- ct.to.cst[match(fit.up$pairs$second, ct.to.cst[,"cell.state"]),"cell.type"]
filter.idx <- pairs.celltype.first != pairs.celltype.second & !fit.up$pairs$second %in% low.count.cst
fit.up[[1]] <- fit.up[[1]][filter.idx] fit.up[[2]] <- fit.up[[2]][filter.idx,]
output.up <- scran::combineMarkers(fit.up$statistics, fit.up$pairs, pval.type = "all", min.prop = NULL, log.p.in = FALSE, log.p.out = FALSE, full.stats = FALSE, pval.field = "p.value", #******* could potentially changed from "p.value" to "FDR" effect.field = "logFC", sorted = FALSE, BPPARAM = BiocParallel::MulticoreParam(n.cores)) #******* I added MulticoreParam
all.ct <- unique(ct.to.cst[,"cell.type"])
ct.stat.list <- lapply(all.ct, function(ct.i){ cst.i <- ct.to.cst[ct.to.cst[,"cell.type"] == ct.i,"cell.state"] output.up.i <- output.up[cst.i]
pval.up.i <- do.call(cbind,lapply(output.up.i, '[', "FDR")) #******* I changed from "p.value" to "FDR" pval.up.min.i <- apply(pval.up.i, 1, min)
lfc.i <- apply(do.call(cbind,lapply(output.up.i, function(output.up.i.j) { apply(output.up.i.j[,grepl("logFC",colnames(output.up.i.j)),drop = FALSE], 1, min) })), 1, max)
data.frame(pval.up.min = pval.up.min.i, min.lfc = lfc.i) })
names(ct.stat.list) <- all.ct return(ct.stat.list) }
End off the code :::::::::::::::