Question: Paired analysis setup using edgeR - follow up question
0
10 days ago by
candida.vaz30
Singapore
candida.vaz30 wrote:

Dear Bioconductor support team,

This question is a follow up of the same I asked before, and then I carried out some analysis, but need help:-

I have samples from two time points. The samples belong to two categories: Treatment and Placebo as follows: 005v1 Placebo 005v4 Placebo 008v1 Treatment 008v4 Treatment

The ID number 005 etc represent the same sample at two time points. At V1 time-point there is no treatment given, its just that the samples are assigned to be given a treatment (treatment) or not (placebo). At V4 time-point there would be difference in placebo and treatment sample.

Now, if the question is to identify the changes caused due to the treatment, how should I define the design? Should V1 time-point be taken into consideration or not, given that the treatment intervention was given only after the assignment of samples at V1.

It seems like a paired analysis, where the blocking factor is the IDs of the samples, so using the edgeR guide (3.4), I tried to block on the ID (called as Cid in the Rcode) in the following way:-

Creating DGElist

raw <- DGEList(counts=counts, group=Group) raw$samples <- cbind(raw$samples,Cid,Visit) data.frame(Sample=colnames(raw),Group,Cid,Visit)

Filter

plot (table (rowSums(cpm(raw) >= 1) ) ) keep <- (rowSums(cpm(raw) >= 1) >= 17) raw <- raw[keep, , keep.lib.sizes=FALSE]

Normalization

norm <- calcNormFactors(raw)

design <- model.matrix( ~0 + Group + Visit, data=norm$samples) v <- voomWithQualityWeights(norm, design, plot=T) dupCor <- duplicateCorrelation(v, design, block=norm$samples$Cid) v <- voomWithQualityWeights(norm, design, block=norm$samples$Cid, correlation=dupCor$consensus, plot=T) fit <- lmFit(v, design, block=norm$samples$Cid, correlation=dupCor$consensus) cm <- makeContrasts(treat = Grouptreatment-Groupplacebo, levels=design ) fit2 <-contrasts.fit(fit, cm) efit <-eBayes(fit2) Thanks a ton! Appreciate you kind support ADD REPLYlink modified 9 days ago • written 10 days ago by candida.vaz30 1 1. I would fuse Group and Visit into a single variable and use a one-way layout. Your additive model assumes that the treatment and time are additive, but this won't be true for the most interesting genes. 2. It wouldn't hurt to run duplicateCorrelation again after the second voomWithQualityWeights call. 3. Your filtering looks odd to me. Why keep 17? What happens to a gene that is not expressed in any group except the treated condition at time point 4 (which presumably would only have non-zero CPMs in a subset of the 17)? ADD REPLYlink modified 10 days ago • written 10 days ago by Aaron Lun22k Thanks a ton Aaron! 1. To fuse into a single variable can I use: Group <- factor(paste(Categ, Visit, sep=".")) raw <- DGEList(counts=counts, group=Group) raw$samples <- cbind(raw$samples,Cid,Visit) (Where Categ =placebo/treatment & Visit=two time points) 2. For filtering (you are right, i'm now using 8 which is the number of samples in treated category)...is this ok? keep <- (rowSums(cpm(raw) >= 1) >= 8) raw <- raw[keep, , keep.lib.sizes=FALSE] 3. Design and DGE calculation design <- model.matrix( ~0 + Group, data=norm$samples) # one way layout

v <- voomWithQualityWeights(norm, design, plot=T)

dupCor <- duplicateCorrelation(v, design, block=norm$samples$Cid)

v <- voomWithQualityWeights(norm, design, block=norm$samples$Cid, correlation=dupCor$consensus, plot=T) dupCor <- duplicateCorrelation(v, design, block=norm$samples$Cid) # second duplicateCorrelation fit <- lmFit(v, design, block=norm$samples$Cid, correlation=dupCor$consensus)

cm <- makeContrasts(treat = Grouptreatment.MV4-Groupplacebo.MV4, levels=design ) # contrast i'm interested in

fit2 <-contrasts.fit(fit, cm)

efit <-eBayes(fit2)

One important fact I would like to stress on, in the experimental design is that some samples remain placebo (9X2) at both time points and some (8X2) get treatment at 2nd time-point. Hope the model is sufficient to get DE genes (Treatment vs Placebo).

Appreciate your time and help on this, getting to learn a lot!

1

If you want to identify differences between treatment and placebo, the most relevant contrast seems to be:

makeContrasts((Grouptreatment.MV4 - Grouptreatment.MV1)
- (Groupplacebo.MV4 - Groupplacebo.MV1), levels=design)


This corresponds to the contrast for edgeR that I proposed in my original answer. It should be better than your cm, as a within-patient comparison avoids the penalty that duplicateCorrelation must impose (i.e., an increase in the standard error, and thus a decrease in power) when comparing across patients.

Dear Aaron,

Thanks a ton for your suggestions!