DESeq2 coutoff or LFC

Question

[anna]

This may be a simple question but i don't know how to answer (I'm new in bioinformatics). Briefly, a couple of months ago I did my analysis of differential expression using DESeq2 and my last command line to obtain de differentially expressed genes between my two conditions was: results <- subset(results, padj < 0.05). After all the next data interpretation i wrote my final work and in the section of materials and methods i described it this way: Differential expression analyses were performed using DESeq2 R-package (v.1.22.1) with a P < 0.05 (Love, 2014). But then one of my reviewers asked me why hadn't i considered a cutoff of fold change to get DEG which is typical in this kind of analysis?

I have to mention that in the results of the DESeq2 i have something like this: LFC > 0 (up) : 138, 100% LFC < 0 (down) : 0, 0%

Does that mean that my LFC was =0? Is this value a default parameter? because i never set it.

I honestly have no idea what to answer because after reading the paper of DESeq2 i thought that the padj (or P value) was like setting the fold change cutoff. Can someone explain this to me? or my results are wrong? if the p value is different to the fold change cutoff how can i justify that i used the p value? Thanks in advance.

Answer 1

[prabin.dm]

This is discussed in the paper Love, Genome Biology. 2014

They mention the default approach of DESeq uses a null hypothesis of 0 LFC. In this case, each gene is compared across treatment for this null hypothesis, and considered significantly different if smaller than the FDR. When you use a fold change cut off (of lets say 1.5) or LFC 0.584, your null hypothesis is that genes between the experiment differ more than 1.5 fold change, and you will be testing that hypothesis with FDR cut off.

Please see the paper and DESeq vignette LFC threshold for more details. The reviewers are correct, as you would see differentially expressed genes with LFC 0 would may not be biologically significant.

[anna]

Thanks for your answer. I just checked all my analysis and i found out that if i set the LFC to 0.5 or 1, i get the same results as if i set the LFC to 0. Do you think that there's a mistake or this could be right? Because this way i can say that i used the LFC=0.5 or 1 which will have more sense than saying it was 0. right? Also i was thinking that setting an LFC=2 would give much fewer genes than i already have (which are not so much).

[prabin.dm]

To obtain meaningful comparison at specific LFC, you have to use lfcThreshold and altHypothesis in the result. Example resGA <- results(dds, lfcThreshold=.5, altHypothesis="greaterAbs") . It might be useful if you added your codes to for us to see.

[anna]

Table1 <- read.csv("htseqcount.csv", header=T, sep=",", row.names="Genes")

colData <- data.frame(condition=factor(c(rep("14",2), rep("9",2))))

colData$condition <- relevel(colData$condition, ref="14")

dds14 <- DESeqDataSetFromMatrix(countData=countsTable1, colData = colData, design=~condition)

dds14 <- DESeq(dds14)

results = results(dds14)

resultsAlpha <- subset(results, padj < 0.05)

#####To get up and down regulated genes:
up14 <- subset(results, log2FoldChange > 1)

down14 <- subset(resultados14vs9prueba8alpha, log2FoldChange < 1)

[prabin.dm]

the above example I have cited is the correct way of getting LFC cut off. Please see Michael Love's response here explanation and see the section "Hypothesis tests with thresholds on effect size" in the DESeq2 paper

Although, the approach you have taken is more common than expected, you might want to use the correct approach.