Question: Limma: difference between groups of paired data
gravatar for johnmcma
2.5 years ago by
johnmcma10 wrote:

So let us assume in addition to the setup in Part 9.4.1 of the User Guide, I'm also adding a variable called Schedule, referring to two types of scheduling for the tested treatment. So the targets frame would look like this:

FileName SibShip Treatment Schedule
File1 1 C A
File2 1 T A
File3 2 C B
File4 2 T B
File5 3 C A
File6 3 T A
File7 4 C B
File8 4 T B
File9 5 C A
File10 5 T A
File11 6 C B
File12 6 T B

What we wanted to identify are genes whose in-sibling C/T differences are different from A and B. How do I handle contrast, and should I use blocking or `duplicateCorrelation`?

Thanks a lot for the answer.

limma limma paired analysis • 1.3k views
ADD COMMENTlink modified 2.5 years ago by Aaron Lun25k • written 2.5 years ago by johnmcma10
Answer: Limma: difference between groups of paired data
gravatar for Aaron Lun
2.5 years ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

For your specific contrast of interest, I would set up the design matrix to block on the sibling factor:

data$SibShip <- factor(data$SibShip) # Make sure it is a factor
design <- model.matrix(~0 + SibShip + Treatment:Schedule, data)
design <- design[,-grep("TreatmentC", colnames(design))] # Get to full rank
colnames(design) <- sub(":", "_", colnames(design)) # Rename for convenience

Now, if we have a look at the column names of the design matrix, we get:

[1] "SibShip1"             "SibShip2"             "SibShip3"            
[4] "SibShip4"             "SibShip5"             "SibShip6"            
[7] "TreatmentT_ScheduleA" "TreatmentT_ScheduleB"

The first 6 coefficients are the sibling-specific effects, which are not of interest to us. The last two coefficients represent the effect (i.e., log-fold change in expression) of treatment over control in schedules A and B, respectively. If you want to find differences in the treatment effect between schedules, you can set up a contrast like this:

con <- makeContrasts(TreatmentT_ScheduleA - TreatmentT_ScheduleB, levels=design)

Obviously, you can also drop either coefficient separately to determine if there is a non-zero treatment effect for the corresponding schedule. This is often a good idea to help in interpreting the log-fold change from the above contrast, e.g., is a positive log-fold change due to an increase in the treatment effect for schedule A or a decrease in the effect for schedule B?

You could also use duplicateCorrelation here. However, this is unnecessary as the calculation of the treatment effect is done within the blocking factor, such that the sibling effect cancels out naturally. If you wanted to compare, say, treatment samples directly between schedules, then you would have to block on sibling in duplicateCorrelation instead of in the design matrix.

ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by Aaron Lun25k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 299 users visited in the last hour