neurorestore / Libra

MIT License
153 stars 25 forks source link

I change the param in edgeR-LRT, and get a better result. #18

Open huanananan opened 2 years ago

huanananan commented 2 years ago

When I use your code in "pseudobulk_de.R" to analysis the DE in my single cell RNA-seq data, I found that the result always not corresponding with pheatmap, featureplot and count file. Then I modified the code as below, it became a correct result.

I'm not sure it is correct, but I think we should discuss about it.

y <- DGEList(counts=CompareMatrix,group=CompareGroup)

keep <- filterByExpr(y)

y <- y[keep,,keep.lib.sizes=FALSE]

y <- calcNormFactors(y,method = 'TMM')

#design = model.matrix(~CompareGroup, data = CompareMetaData)

# change the design as below:
# control treatment
# 1 0
# 1 0
# 0 1
# 0 1
design <- data.frame(shame=append(rep(1,control_cell)),
                                  rep(0,length(treatment_cell))),
                     clp=append(rep(0,length(control_cell)),
                                rep(1,length(treatment_cell))))
design$shame = as.array(design$shame)
design$clp = as.array(design$clp)
design=as.matrix(design)

y <- estimateDisp(y,design)

fit <- glmFit(y,design = design)
# to compare 2 vs 1
# the coef=2 or coef=2:1 not return a correct result. And whatever I change, I couldn't change the order of compare (the gene up-regulate in treatment cell always get a minus logFC). But contrast=c(-1,1) return a result corresponding with pheatmap, featureplot and count file.
lrt <- glmLRT(fit, contrast=c(-1,1))

topTags(lrt)
jordansquair commented 2 years ago

If I understand correctly you are saying the direction of your logFC is backwards? If so you can just set your factor levels prior to run_de - or flip the sign on your logFC.