This is a test version of Biostars. For the public version, visit https://www.biostars.org.
DESeq2 not producing the correct log2fc
#packages here 
library(matrixStats) 
library(DESeq2)
library(dplyr)
#### 
> Expdesign <- read.csv(file = "design.csv", sep = ",", header = TRUE) 
> head(Expdesign[1,]) 
> colnames(Expdesign)[1] <- "SampleName"
> head(Expdesign[1,])

> options(scipen=999)
> InputData <- read.csv(file = "counts.csv", sep = ",", header = TRUE, row.names = 1) 
> head(InputData[1,]) 
> InputData <- as.matrix(as.integer(InputData), stringsAsFactors = TRUE)


> DESeqTable <- DESeqDataSetFromMatrix(countData = InputData, colData = Expdesign, design = ~treatment)

> DESeqTable$treatment <- relevel(DESeqTable$treatment, ref = "CTRL")

> dds <- DESeq(DESeqTable, betaPrior = FALSE, minRep = Inf)
> res <- results(dds)

> write.csv(res, file = "TREATEDvsCTRL4.csv")

Above is the code that I have written for DESeq2. When I manually do the math for log2fc, it is different from the values that was produced through this code. I have set the correct colData and rowNames by using a TFmatch code to make sure. I have no idea why the values are not producing correctly.

>TFmatch <- all(DESeqTable@colData@rownames == DESeqTable$SampleName

>TRUE

Here I am checking the control variable

>levels(DESeqTable$expvar)

>CTRL TREAT

For example, this is the first row:

  • KLF6
  • Mean of CTRL: 43.97683921
  • Mean of TREAT: 136.555425
  • log2FoldChange: 0.324728758

This should be 1.245. Please help.

bioconductors deseq2 r
• 1,580 views • link
written by marianne • marianne

The fold changes are not the naive difference between the group means. There are additional factors that play a role. Can you use plotCounts for that gene then I can try to explain with an example? In short, if genes have low counts and/or lots of variability then the fold changes are shrunken towards zero as the observed mean has large standard errors and the given information (sample size) is not sufficient to confidently give the "calculated" (by mean) fold change. I bet that gene has a low baseMean. Check with the results table and also use plotMA().

• 0 views • link
written by ATpoint ATpoint

That gene is not the only gene that has a incorrect log2FoldChange, the whole matrix is producing the wrong values. The base mean for the KLF6 gene is 87.40739865. That's rather high than the rest of the genes in the matrix. Warning though as this is my first production of anything in R, so I know I am tackling a hard task.

• 0 views • link
written by marianne • marianne

0 answers

No answers yet.

Log in to answer this question.

More posts like this