Limma: difference between groups of paired data
1
0
Entering edit mode
johnmcma ▴ 10
@johnmcma-12132
Last seen 2.5 years ago
United States

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 paired analysis limma • 2.5k views
ADD COMMENT
2
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 16 hours ago
The city by the bay

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 COMMENT

Login before adding your answer.

Traffic: 585 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6