Question

limma - Multi-level Experiments - time course data

0

Entering edit mode

SPbeginner • 0

@spbeginner-15170

Last seen 5.7 years ago

Dear Bioconductor community,

I have a two-way experiment in which time course profiles are to be compared for two many groups.

The two experimental dimensions or factors here are Treatment (A-D-F) and Time. For covenience, as explain in Limma user manual, I used a combinaison of Treatment/Time = Group in the target frame.

Animals received different treatment (A, D or F). But for some reasons some samples are missing in some groupes

Should I use a blocked model on animals or is it a problem ?

I want to identify :

1/ Which genes respond at either the Day1, Day3 or Day27 in the different group ( and D1, D5 , D7 for Group F)

2/ Which genes respond differently over time in the different groups (D3, D7, Day 27/28)

The results for gene set enrichment analysis I obtained for comparison over time in Group F perplex me. This is why I'm wondering if my analysis is right.

Thanks in advance for your help

limma limma voom rnaseq_analysis DEG • 2.2k views

ADD COMMENT • link 6.7 years ago • updated 5.7 years ago SPbeginner • 0

0

Entering edit mode

One question to consider is whether all the 0 time points should be considered a single group, since presumably these are all pre-treatment, so the identity of the treatment is irrelevant. However, it's possible that in this experiment the different treatments require different initial conditions, in which case it would not be correct to combine the 0 time points.

ADD REPLY • link 6.7 years ago Ryan C. Thompson ★ 7.9k

score 1 · Answer 1 · 2018-03-29

1

Entering edit mode

Aaron Lun ★ 28k

@alun

Last seen 11 hours ago

The city by the bay

The code looks fine to me, assuming that you consider F's timepoint D28 to be equivalent to D's D27.

Note that your comparisons are all "within animal" comparisons; you are either comparing each timepoint to zero, or you are comparing the differences from zero between treatments, neither of which requires a direct comparison between different animals. This means that you can directly block in the design matrix:

design <- model.matrix(~0 + Animal + Group, Design_clean_simple)
design <- design[,!grepl("\\.0$", colnames(design))]

... which gives you a bunch of GroupX.Y terms that directly represent the log-fold change in expression between timepoint Y and 0 for Treatment X. This provides an alternative to using duplicateCorrelation, which can be anti-conservative at times due to the use of a consensus correlation for all genes.

I can't give you any specific advice for your gene set analysis, though I will mention that using treat with lfc of 0.5-0.8 can focuses on more biologically relevant genes with strong log-fold changes. This may improve the interpretability of your results.

ADD COMMENT • link 6.7 years ago Aaron Lun ★ 28k

0

Entering edit mode

Given the sparseness of which animals received which treatment group, I'd be inclined to stick with duplicateCorrelation here, as in OPs code. If there are very big differences between the animals however (with one set of animals different to others), then I would include animal in the design matrix instead.

ADD REPLY • link 6.7 years ago Gordon Smyth 52k

0

Entering edit mode

The Group A and D received an attenuated form of a virus whereas Group F received the virus. So A and D are close whereas F is really different.

So i'll include animal in the model matrix as Aaron Lun suggest :

design <- model.matrix(~0 + Animal + Group, Design_clean_simple)
design <- design[,!grepl("\\.0$", colnames(design))]

dge <- DGEList(counts=Count)
keep <- filterByExpr(dge)
dge.filtered <- dge[keep,,keep.lib.sizes=FALSE]

# calculate norm factor for libraries
dge.norm <- calcNormFactors(dge.filtered)

v.wts.bl<- voomWithQualityWeights(dge.norm, design, plot=TRUE)
colnames(design)<-make.names(colnames(design))
fit.wts.bl <- lmFit(v.wts.bl, design)
fit.wts.bl <- eBayes(fit.wts.bl)

As Aaron mention, I have a bunch of GroupX.Y terms that represent the log-fold change in expression between time Y and 0 for treatment X.
Now, to have comparison between groups :

contrast.matrix<-makeContrasts(GpAvsGpF.D3=GroupA.3-GroupF.3,
                                       GpAvsGpF.D7=GroupA.7-GroupF.7,
                                       GpAvsGpF.D27=GroupA.27-GroupF.28,
                                       GpDvsGpF.D3=GroupD.3-GroupF.3,
                                       GpDvsGpF.D7=GroupD.7-GroupF.7,
                                       GpDvsGpF.D27=GroupD.27-GroupF.28, 
                                       levels=colnames(design))

fit2.wts.bl <- contrasts.fit(fit.wts.bl, contrast.matrix)
fit2.wts.bl < eBayes(fit2.wts.bl,robust=TRUE)

But I have the following error :

Error in fit2.wts.bl < eBayes(fit2.wts.bl, robust = TRUE) : 
  comparison of these types is not implemented

It seems that I cannot use voomWithQualityWeights() in that particular comparison ?

I tested with

voom(dge.norm, design, plot=TRUE)

and it worked. I used voomWithQualityWeights() because quality of some samples are really bad.
The weigth in voomWithQualityWeights() correlates well with sample QC !

ADD REPLY • link 6.7 years ago SPbeginner • 0

0

Entering edit mode

Feel free to include animal in the design matrix if you like, as per Aaron's suggestion, but I think you have misunderstood my advice.

The choice of duplicateCorrelation vs blocked design matrix depends on how different the animals are (supposing they received the same treatment), not on how similar the groups are (A vs F etc).

ADD REPLY • link 6.7 years ago Gordon Smyth 52k

0

Entering edit mode

Animals are supposed to be the same (same sex + same age, all coming from the same farm).

ADD REPLY • link 6.7 years ago SPbeginner • 0

0

Entering edit mode

If you read your error message, you will see that you used < instead of <-. So fix that.

And as Gordon said, the choice between blocking on the design and blocking in duplicateCorrelation depends on whether there are systematic differences between animals in a given treatment condition. For example, you mention that samples in F were collected across two studies; this could result in differences between animals from different studies. You can check this by creating a MDS plot using the F samples only.

ADD REPLY • link 6.7 years ago Aaron Lun ★ 28k