Skip to content
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.

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

Is the function: get.exp.stat() using p.value or FDR as default #61

Open
Aeget1000 opened this issue Oct 25, 2023 · 0 comments
Open

Is the function: get.exp.stat() using p.value or FDR as default #61

Aeget1000 opened this issue Oct 25, 2023 · 0 comments

Comments

@Aeget1000
Copy link

Aeget1000 commented Oct 25, 2023

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 :::::::::::::::

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant