Hi

I'm using DESeq2 to analyze raw read counts. I have 5 groups and 3 batches (different projects). In order to find out which genes are expressed differentially in all groups and also adjust for batch effects I conducted the LRT test step by step from the guideline. This is my code:

rna.o=DESeqDataSetFromMatrix(data,colD,design = ~projects+subtypes) rna.d=DESeq(rna.o,parallel = T,reduced = ~projects,test = "LRT") rna.res=results(rna.d,parallel = T)

as a result this code gave me thousands of highly significant genes. I took normalized counts and removed batch effects by comBat in SVA package to plot a heatmap of these significant genes. no significant differences was observable even in genes with adjusted p-values less than 10^-100. it is evident in this heatmap: https://drive.google.com/file/d/1EhIWQXYAhH1LNy5YJlLX6EG0RSA4Qpiq/view?usp=sharing

below link is a heatmap before the batch effect removal and the batch effect is completely evident: https://drive.google.com/open?id=1RSq3AaR9ucL_a03y4KAmgy17J0iFyePC

so i decided to use plotCounts to find out how much is the difference in groups for my top genes. this is the result: https://drive.google.com/open?id=18ue-oiTmdlFWrNbRXRQDTpUD4hKOkuGh this is the code:

plotCounts(rna.d, gene=rownames(rna.res[1,]), intgroup="Subtypes") . . . plotCounts(rna.d, gene=rownames(rna.res[5,]), intgroup="Subtypes")

It is completely evident that there are minor difference over groups in top 5 genes (p-value<10^-200). In the next step I plotted counts like this:

plotCounts(rna.d, gene=rownames(rna.res[1,]), intgroup=c("projects","Subtypes")) . . . plotCounts(rna.d, gene=rownames(rna.res[5,]), intgroup=c("projects","Subtypes"))

this is the result: https://drive.google.com/open?id=18ue-oiTmdlFWrNbRXRQDTpUD4hKOkuGh

there are more differences in combination of groups and projects but I can't see enough difference for a p-value like that. what is the problem and what should I do? should i just normalize by DESeq2, remove batch effect and move to limma?

For heatmaps, we have recommended code in the workflow. It’s hard to see anything on the count scale in a heatmap, we have recommended transformations to use.

In your first counts plot, I see differences across groups. In particular group 2 often appears different than the rest. The pvalue is for the null that all groups are equal. So I’m not sure what’s the issue exactly.

If you want to find genes where groups 1, 3-5 are different from each other you could remove group 2 and perform a similar LRT.