Entering edit mode
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
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....
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 :
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.
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).
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...
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
I solved the problem. It should have been
resLFCalso outside of bracket like this: