DESeq2 different log2FC and p-value filtering parameters results comparison
1
0
Entering edit mode
User000 • 0
@ea03770f
Last seen 12 days ago
Italy

Hello,

I have created two datasets with two different filtering parameters as following:

log2folgchange = 0:

contrast = (list(c("treated_1","treated_2"),("control"),listValues = c(1/2,-1)))
resLFC <- lfcShrink(dds, contrast=contrast, res=res, type="ashr")
up <- res[which(resLFC$log2FoldChange > 0 & resLFC$padj < 0.01),]
down <- res[which(resLFC$log2FoldChange < 0 & resLFC$padj < 0.01),]


log2foldchange = 1:

up <- res[which(resLFC$log2FoldChange > 1 & resLFC$padj < 0.01),]
down <- res[which(resLFC$log2FoldChange < -1 & resLFC$padj < 0.01),]


In the dataset at log2FC=0 I see some genes that are not in the dataset of log2FC=1 even thought the log2FC is equal to -1.29 and padj is 0.005, si it should be in my dataset with log2F= 1 I presume. Why does that happen?

gene    baseMean    log2FoldChange  lfcSE   stat    pvalue  padj
ENSG000000XXXX  74,9062387275475    -1,29048864731255   0,377167404326021   -3,4215275034665    0,000622704264736   0,00590226764893

log2 DESeq2 • 257 views
0
Entering edit mode
Basti ▴ 390
@7d45153c
Last seen 1 day ago
France

I think you made something wrong, try setdiff(up2$gene,up$gene)with up <- res[which(resLFC$log2FoldChange > 0 & resLFC$padj < 0.01),] and up2 <- res[which(resLFC$log2FoldChange > 1 & resLFC$padj < 0.01),]

0
Entering edit mode

I checked also with excel before and now with your code, In down log2FC=0 I have 1371 genes and in log2FC=1 616, with your code I get 616, but I do not see the gene in the list anyway....

0
Entering edit mode

Check your objects are correct and re-run your code entirely, there might be an overlap between different objects because you can't have more genes applying the threshold of log2FC< -1 than log2FC<0, this is strictly impossible. If you try with an random dataset, you can see that what you are observing can't happen :

set.seed(1)
dds <- makeExampleDESeqDataSet(n=500,betaSD=1)
dds <- DESeq(dds)
res <- results(dds)

resultsNames(dds)
resLFC <- lfcShrink(dds=dds, coef=2, type="ashr")

up <- res[which(resLFC$log2FoldChange > 1 & resLFC$padj < 0.01),]
down <- res[which(resLFC$log2FoldChange < -1 & resLFC$padj < 0.01),]

up2 <- res[which(resLFC$log2FoldChange > 1 & resLFC$padj < 0.01),]
down2 <- res[which(resLFC$log2FoldChange < -1 & resLFC$padj < 0.01),]

setdiff(up2$gene,up$gene)

0
Entering edit mode

Thanks for your answer, infact I do not have more genes in log2FC=-1 dataset, as I wrote above for log2FC=0 I have 1371 genes and log2FC=-1 I have 616 genes. All log2FC values higher than -1 that are not the in log2FC=-1 dataset (around 125 genes) do not exceed log2FC value more than 1.5/-1.5 for up and down, respectively. I think there is some explanation in the DESeq2 algorithm.

0
Entering edit mode

There is nothing related to DESeq2 in your code, this is a simple filter on a DESeq2 result. Have you ran my example code ?Certainly there is something wrong in your code because what you observe is mathematically impossible (if you add a small example of your dataset in your post I would be glad to help).

0
Entering edit mode

I get 755 genes, exaclty the difference between 1371 and 616. But the problem is that some genes according to my filtering should be in the log2fc=1 dataset and they are missing. So I do not understand how this can solve the problem...

0
Entering edit mode

If you share example dataset we could have a look, otherwise there us something wrong with your code that we can't solve with information available for us

0
Entering edit mode

I solved the problem. It should have been resLFC also outside of bracket like this:

up <- **resLFC**[which(resLFC$log2FoldChange > 0 & resLFC$padj < 0.01),]
down <- **resLFC**[which(resLFC$log2FoldChange < 0 & resLFC$padj < 0.01),]