Hi, I like to confirm if the statistical model of Deseq2 is appropriate for my experiment design. Here’s my code.
> dds<-DESeqDataSetFromMatrix(countData=countsTable,colData,formula(~ Transfection*Enrichment )) > colData(dds)$Transfection<-factor(colData(dds)$Transfection,levels=c("gfp","kd")) > colData(dds)$Enrichment<-factor(colData(dds)$Enrichment,levels=c("total","nascent")) > dds$group<-factor(paste0(dds$Transfection, dds$Enrichment)) > design(dds)<-~group > dds<-DESeq(dds,test = "Wald") > resultsNames(dds) [1] "Intercept" "groupgfpnascent" "groupgfptotal" "groupkdnascent" "groupkdtotal"
I generated the following results.
> res_KDvsGFPinNascent<-results(dds, contrast=list(c("groupkdnascent"),c("groupgfpnascent"))) > res_KDvsGFPinTotal<-results(dds, contrast=list(c("groupkdtotal"),c("groupgfptotal")))
But I also want to calculate the fold change and the pvalue of (KDvsGFPinNascent)/(KDvsGFPinTotal). I use the following codes.
res_KDnGFPt_GFPnKDt<-results(dds, contrast=list(c("groupkdnascent","groupgfptotal"),c("groupgfpnascent","groupkdtotal")))
I have confirmed that the fold change of “res_KDnGFPt_GFPnKDt” is equal to (“res_KDvsGFPinNascent”)/ (“res_KDvsGFPinTotal”).
But I got confused, because I found that moslty p<0.05 genes of “res_KDnGFPt_GFPnKDt” were not overlapping with p<0.05 genes in "res_KDvsGFPinNascent" and "res_KDvsGFPinTotal".
I am not sure which one is correct. Please enlighten me how do I calculated the pvalue based on (KDvsGFPinNascent)/(KDvsGFPinTotal) in DEseq2 correctly.
Thank you very much.
And the opposite can be true as well: you can have weak, positive LFC in Nascent and weak, negative LFC in Total, where neither passes your FDR threshold, although the test of whether the LFC is different could pass your FDR threshold.