lusystemsbio / NetAct

Network construction Algorithm
MIT License
17 stars 5 forks source link

bug in RNAseqDegs_DESeq function #5

Open yeungj234 opened 1 year ago

yeungj234 commented 1 year ago

Hello, I have encountered a bug when using DESeq instead of limma for differential gene analysis step.

Here is the error message with batch information included: image

Here is the error message without batch information included: image

If I use the same exact code and phenoData information but with RNAseqDegs_limma, the function works fine and I don't encounter this issue.

I tried to fix the issue by looking at your source code for RNAseqDegs_DESeq and changing it to this:

RNAseqDegs_DESeq = function(counts, phenodata, complist, qval = 0.05) { celltype = phenodata$celltype compList <- complist DEresult = vector(mode = "list", length = length(compList) + 1) names(DEresult) = c(compList, "Overall")

dds <- DESeqDataSetFromMatrix(countData = counts, colData =phenoData@data , design = ~ celltype) 
dds <- estimateSizeFactors(dds)
e <- counts(dds, normalized = TRUE)

for (comp in compList) {
  cs = strsplit(comp, split = "-")[[1]]
  ids = which(celltype %in% cs)
  tmpcounts = counts[, ids]
  tmpph = as(phenodata[ids, ], "data.frame")

  if (is.null(phenodata$batch) | length(unique(phenodata$batch)) == 1){
    dds = DESeqDataSetFromMatrix(countData = data.frame(tmpcounts), colData = tmpph, design = ~ celltype)
  }else{
    dds = DESeqDataSetFromMatrix(countData = data.frame(tmpcounts), colData = tmpph, design = ~ batch + celltype)
  }
  dds = DESeq(dds)
  DEtable = results(dds, alpha = qval)
  DEtable =  data.frame(DEtable[complete.cases(DEtable), ])
  DEtable = DEtable[order(DEtable$pvalue), ]
  rank_vector = abs(DEtable$stat); names(rank_vector) = rownames(DEtable)
  degs = rownames(DEtable[DEtable$padj < qval, ])
  tmpDErslt <- list(table = DEtable, rank_vector = rank_vector, degs = degs)
  DEresult[[comp]] = tmpDErslt
}
DEresult$Overall = list(table = DEtable, rank_vector = rank_vector, degs = degs, e = e)
return(DEresult) }

This code then allows for DESeq to run properly if only celltype information is provided in phenoData but not batch information. I'm not sure how to fix the problem at this point but since my other analyses have been using DESeq and not limma, I would like to use NetAct with DESeq preferably.

yeungj234 commented 1 year ago

note: doing DESEq with the code I provided doesn't allow you to proceed further with downstream steps since I get this error when trying to run the TF_Activity function(). Please advise. image

iamcorrinne commented 1 year ago

I also receive the same error message when trying to use RNAseqDegs_DESeq, and I am interested in a fix for this.

iamcorrinne commented 1 year ago

@vivekkohar @lusystemsbio any insight here?

YukaiYou commented 1 year ago

note: doing DESEq with the code I provided doesn't allow you to proceed further with downstream steps since I get this error when trying to run the TF_Activity function(). Please advise. image

We have updated this function, and it should work normally after reinstallation. By the way, I think you should use cell type without batch

YukaiYou commented 1 year ago

@vivekkohar @lusystemsbio any insight here?

I also receive the same error message when trying to use RNAseqDegs_DESeq, and I am interested in a fix for this.

After reinstalling, the function should work properly since we have implemented updates.