Question

DESeq2 log2foldchange shows opposite sign than clearly visible in pseudocounts or vst values (correct order in contrast)

0

Entering edit mode

ligia.mateiu • 0

@ligiamateiu-9214

Last seen 4.5 years ago

European Union

I'm analyzing a large dataset of 250 samples and I'm confronted with some results I cannot explain. I've added the code below. The analysis for this data gave only one gene to have significant padj (for which the log2fc sign was clearly correct). Still, we are interested in looking at the lfc for all the genes, so we investigated the logratios file. Maybe it is important to mention that none of the log2foldchanges in the data is >|0.6|

We noticed that the log2foldchange value, for many genes, does not reflect the situation in pseudo-counts (or vst) values. The contrast is respecting the indicated order of the groups, so it is not flipped order. The gene I picked below is one of the most obvious that sth is not right.

(I have tried using lfcShrink function with apeglm as well and it is the same situation). Does anybody has an idea why?

deseqdata <- DESeqDataSetFromMatrix(countData = countsraw,colData = sampleinfo,design = ~Age_categ+RIN_categ+Sex+plate+grouping1)
deseqdata$grouping1 <- factor(deseqdata$grouping1, levels = c("groupA","groupB","control"))

filter<-rowSums(counts(deseqdata) >= 10) >= 5
dds<-deseqdata[filter,]
dds<-DESeq(dds,betaPrior=TRUE)
vst<-vst(dds, blind=F)
res <- results( dds, contrast = c("grouping1","groupA" ,"groupB") ,addMLE=TRUE,alpha = 0.05)

the pseudocounts for the 3 groups are below, the basemean, the shrunken and unshrunken log2FC and pvalue groupA groupB groupC baseMean log2FoldChange lfcMLE lfcSE pvalue padj 1213.9 832.13 671.47 999.3 -0.35 -0.46 0.10 3.09E-04 2.58E-01

(for the vst values it is the same situation) The plot with y=vst values shows the situation https://ibb.co/JK30y8P

deseq2 log2FoldChange • 1.6k views

ADD COMMENT • link updated 4.5 years ago by Michael Love 43k • written 4.5 years ago by ligia.mateiu • 0

score 0 · Answer 1 · 2020-05-29

0

Entering edit mode

Michael Love 43k

@mikelove

Last seen 2 days ago

United States

results is giving you negative, indicating that A is less than B but your violinplot and counts show the opposite.

Just to double check, can you make a plot with plotCounts? It could be that the sample table was not provided to DESeq2 in the order of the count table.

ADD COMMENT • link 4.5 years ago Michael Love 43k

0

Entering edit mode

Hi Mike, thanks a lot for such a swift response. It would have been glad to be that mistake!

 dataplot<-plotCounts(dds,"gene",intgroup = "grouping1",normalized = TRUE,returnData = TRUE)
 p<-ggplot(dataplot,aes(x=grouping1, y=count,fill=grouping1)) +geom_violin(trim=FALSE)+geom_boxplot(width=0.1)

https://ibb.co/NtvhTVp

ADD REPLY • link 4.5 years ago ligia.mateiu • 0

0

Entering edit mode

Can you show the top line when you print:

res

In the R console. What does R say the contrast is?

ADD REPLY • link 4.5 years ago Michael Love 43k

0

Entering edit mode

the order is ok, as I said, the only gene with a good padj was ok.

> head(res)
log2 fold change (MAP): grouping1 groupA vs groupB 
Wald test p-value: grouping1 groupA vs groupB
DataFrame with 6 rows and 7 columns

ADD REPLY • link 4.5 years ago ligia.mateiu • 0

0

Entering edit mode

Oh, for the first time I noticed that you have a lot of covariates.

If you have a confounded design (where covariates are correlated with your grouping), then it can easily happen that the sign of LFC changes compared to what you would guess by looking at marginal distributions.

E.g. replace your design with ~grouping and you will get the expected LFC. However, the problem is the confounding. The best you can do is include the covariates, and you can no longer set your expectations according to what the boxplots look like.

ADD REPLY • link 4.5 years ago Michael Love 43k

0

Entering edit mode

To be honest, that was my guess as well!! I'm using DESeq2 for quite some years already and I did not encouter this situation before. For this particular dataset, I have tried to add some covariates to deal with tissue heterogeneity which indeed are correlated with the grouping. Thanks a lot for your help!!

ADD REPLY • link 4.5 years ago ligia.mateiu • 0