Question

makeContrasts question

0

Entering edit mode

HS ▴ 10

@hakusen03-10599

Last seen 6 months ago

United States

Hi, I have samples collected from two treatments (A and B), and 5 time points (0,1,7,14,22 days). I'm interested in comparing gene expression between treatment B and A at each time point after accounting for expression at time 0 (my code is below). I don't understand why I'm getting a low number of differentially expressed genes for each contrasts. I was expecting a larger number of differentially expressed genes after looking at the MDS plot, especially at time 22. I still got a low number of DEGs after lowering the FC threshold to 1.1. I'm not sure if there is a problem in my code. Is there something wrong with my code? Have I set the design matrix, model and contrasts correctly?

x <- read.delim(file="counts.txt", header=T, com='', sep="\t", row.names=1, check.names=F)
groups <- factor(c(rep("A.0d",2),rep("A.1d",3),rep("A.7d",3),rep("A.14d",3),rep("A.22d",3),rep("B.0d",3),rep("B.1d",3),rep("B.7d",3),rep("B.14d",3),rep("B.22d",3)))
design <- model.matrix(~0+groups)
colnames(design) <- levels(groups)
y <- DGEList(counts=x,group=groups)
keep <- filterByExpr(y,design)
y <- y[keep,,keep.lib.sizes=FALSE]
y <- calcNormFactors(y, method="TMM")
y <- estimateDisp(y,design,robust=TRUE)
fit <- glmQLFit(y, design,robust=TRUE)
my.contrasts <- makeContrasts(BvsA.1d = (B.1d-B.0d)-(A.1d-A.0d),BvsA.7d = (B.7d-B.0d)-(A.7d-A.0d), BvsA.14d = (B.14d-B.0d)-(A.14d-A.0d),BvsA.22d = (B.22d-B.0d)-(A.22d-A.0d),levels=design)
BvsA.1d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.1d"],lfc=log2(1.5))
is.de.1d <- decideTestsDGE(BvsA.1d,adjust.method="fdr", p.value=0.05)
summary(is.de.1d)
BvsA.7d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.7d"],lfc=log2(1.5))
is.de.7d <- decideTestsDGE(BvsA.7d,adjust.method="fdr", p.value=0.05)
summary(is.de.7d)
BvsA.14d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.14d"],lfc=log2(1.5))
is.de.14d <- decideTestsDGE(BvsA.14d,adjust.method="fdr", p.value=0.05)
summary(is.de.14d)
BvsA.22d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.22d"],lfc=log2(1.5))
is.de.22d <- decideTestsDGE(BvsA.22d,adjust.method="fdr", p.value=0.05)
summary(is.de.22d)

       1*A.0d -1*A.1d -1*B.0d 1*B.1d
Down                                                   0
NotSig                                             25915
Up                                                     0

       1*A.0d -1*A.7d -1*B.0d 1*B.7d
Down                                                   0
NotSig                                             25893
Up                                                    22

       1*A.0d -1*A.14d -1*B.0d 1*B.14d
Down                                                     0
NotSig                                               25914
Up                                                       1

       1*A.0d -1*A.22d -1*B.0d 1*B.22d
Down                                                     0
NotSig                                               25912
Up                                                       3



Code should be placed in three backticks as shown below

```r

# include your problematic code here with any corresponding output 
# please also include the results of running the following in an R session 

sessionInfo( )

makeContrasts edgeR • 1.4k views

ADD COMMENT • link updated 3.3 years ago by Steve Lianoglou ★ 13k • written 3.3 years ago by HS ▴ 10

score 0 · Answer 1 · 2021-08-11

0

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 57 minutes ago

United States

Those are pretty large fold changes. You might consider using the default for glmTreat. Put another way, when you use treat you are testing for evidence that the underlying population differences are larger than a particular fold change, so your observed fold changes will end up being much larger than 1.5-fold. Using something like 1.2 fold (the default) is more reasonable, or maybe don't use treat at all, and see what you get without a fold change criterion.

ADD COMMENT • link 3.3 years ago James W. MacDonald 67k

0

Entering edit mode

The number of differentially expressed didn't change much after running a lower fold change threshold (1.1) or glmQLFTest. Should the contrast coefficients be -1A.0d 1A.1d -1B.0d 1B.1d instead of 1A.0d -1A.1d -1B.0d 1B.1d? (We are interested in comparing trt B vs trt A at each time point after accounting for expression at 0d). I think, 1A.0d -1A.1d -1B.0d 1B.1d is correct, but I'm not sure.

HS.

``` BvsA.1d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.1d"],lfc=log2(1.1)) is.de.1d <- decideTestsDGE(BvsA.1d,adjust.method="fdr",p.value=0.05) summary(is.de.1d) 1A.0d -1A.1d -1B.0d 1B.1d Down 0 NotSig 25915 Up 0 BvsA.7d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.7d"],lfc=log2(1.1)) is.de.7d <- decideTestsDGE(BvsA.7d,adjust.method="fdr",p.value=0.05) summary(is.de.7d) 1A.0d -1A.7d -1B.0d 1B.7d Down 1 NotSig 25839 Up 75 BvsA.14d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.14d"],lfc=log2(1.1)) is.de.14d <- decideTestsDGE(BvsA.14d,adjust.method="fdr",p.value=0.05) summary(is.de.14d) 1A.0d -1A.14d -1B.0d 1B.14d Down 0 NotSig 25914 Up 1 BvsA.22d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.22d"],lfc=log2(1.1)) is.de.22d <- decideTestsDGE(BvsA.22d,adjust.method="fdr",p.value=0.05) summary(is.de.22d) 1A.0d -1A.22d -1B.0d 1B.22d Down 1 NotSig 25908 Up 6

#

BvsA.1d <- glmQLFTest(fit, contrast=my.contrasts[,"BvsA.1d"]) degs.1d <- topTags(BvsA.1d) is.de.1d <- decideTestsDGE(BvsA.1d) summary(is.de.1d) 1A.0d -1A.1d -1B.0d 1B.1d Down 0 NotSig 25915 Up 0

BvsA.7d <- glmQLFTest(fit, contrast=my.contrasts[,"BvsA.7d"]) is.de.7d <- decideTestsDGE(BvsA.7d) degs.7d <- topTags(BvsA.7d) summary(is.de.7d) 1A.0d -1A.7d -1B.0d 1B.7d Down 3 NotSig 25817 Up 95

BvsA.14d <- glmQLFTest(fit, contrast=my.contrasts[,"BvsA.14d"]) degs.14d <- topTags(BvsA.14d) is.de.14d <- decideTestsDGE(BvsA.14d) summary(is.de.14d) 1A.0d -1A.14d -1B.0d 1B.14d Down 0 NotSig 25914 Up 1

BvsA.22d <- glmQLFTest(fit, contrast=my.contrasts[,"BvsA.22d"]) degs.22d <- topTags(BvsA.22d) is.de.22d <- decideTestsDGE(BvsA.22d) summary(is.de.22d) 1A.0d -1A.22d -1B.0d 1B.22d Down 1 NotSig 25908 Up 6

ADD REPLY • link 3.3 years ago HS ▴ 10

score 0 · Answer 2 · 2021-08-12

0

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 5 hours ago

WEHI, Melbourne, Australia

You're not actually testing for differential expression but instead testing interactions. Interactions usually give fewer significant results than straight comparisons.

ADD COMMENT • link 3.3 years ago Gordon Smyth 51k

0

Entering edit mode

Hi Gordon,

I'm having a hard time understanding why I'm getting few significantly expressed genes after testing for the interaction in this dataset. For example:

``` A.22vs0d <- glmTreat(fit, contrast=my.time.contrasts[,"A.22vs0d"],lfc=log2(1.5)) is.de.22d <- decideTestsDGE(A.22vs0d,adjust.method="fdr",p.value=0.05) summary(is.de.22d) -1A.0d 1A.22d Down 2349 NotSig 22966 Up 600

B.22vs0d <- glmTreat(fit, contrast=my.time.contrasts[,"B.22vs0d"],lfc=log2(1.5)) is.de.22d <- decideTestsDGE(B.22vs0d,adjust.method="fdr",p.value=0.05) summary(is.de.22d) -1B.0d 1B.22d Down 2066 NotSig 22387 Up 1462

BvsA.22d <- glmTreat(fit, contrast=my.contrasts[,"BvsA.22d"],lfc=log2(1.5)) is.de.22d <- decideTestsDGE(BvsA.22d,adjust.method="fdr",p.value=0.05) summary(is.de.22d) 1A.0d -1A.22d -1B.0d 1B.22d Down 0 NotSig 25912 Up 3

Is it related to testing the interaction with few replicates per treatment combo? Is the variation too small to find anything significant when testing the interaction? Can you please further explain why is harder to find significant genes at the interaction level?

ADD REPLY • link 3.3 years ago HS ▴ 10

1

Entering edit mode

Perhaps Gordon is referring to the same topic Andrew Gelman discusses in this post about testing interactions:

You need 16 times the sample size to estimate an interaction than to estimate a main effect

I’m linking you to a specific reply in the comments that writes out explicitly the inflation of the variance when testing a ratio of ratios (an interaction), but the whole post is worth a read.

ADD REPLY • link 3.3 years ago Steve Lianoglou ★ 13k