Question: Limma analysis, controlling for Donors.
0
9 months ago by
reubenmcgregor880 wrote:

I am using Limma to analyse some microarray data and have managed a simple workflow so far, but am struggling with a slightly more complex design matrix (despite trying to adapt some other answers on here to my question, sorry).

I have phenotypic data of an experiment that looks like this:

    tissue donor anat_areas heart_area    valves MV_only
Ao3     Ao     3    vessels         na        na   other
Ao8     Ao     8    vessels         na        na   other
Ao4     Ao     4    vessels         na        na   other
AV3     AV     3     valves         na sl_valves   other
AV4     AV     4     valves         na sl_valves   other
AV8     AV     8     valves         na sl_valves   other
LA3     LA     3   cavities     atrium        na   other
LA4     LA     4   cavities     atrium        na   other
LA8     LA     8   cavities     atrium        na   other
LV3     LV     3   cavities  ventricle        na   other
LV4     LV     4   cavities  ventricle        na   other
LV8     LV     8   cavities  ventricle        na   other
MV3     MV     3     valves         na av_valves      MV
MV4     MV     4     valves         na av_valves      MV
MV8     MV     8     valves         na av_valves      MV


I am currently comparing the "MV" vs the "other" group as follows:

group <- factor(pData$MV_only) design <- model.matrix(~0+group) colnames(design) <- levels(group) contrast.matrix <- makeContrasts( mvvsother = MV-other,levels=design) contrast.matrix  Contrasts Levels mv_vs_other MV 1 other -1  fit <- lmFit(y, design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2, trend=T) As you can see in the pData there are 3 donors (3, 8 and 4) factor(pData$donor) [1] 3 8 4 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 3 4 8 Levels: 3 4 8

I would like to carry out the same analysis taking into consideration this "donor effect".

In my lay language I guess I would say how do I construct a contrast matrix to "control for Donor variability" or "do a paired analysis"?

Thank you and sorry for a question that has seemingly been answered before. I am just going in circles.

microarray limma R • 324 views
modified 9 months ago by Gordon Smyth38k • written 9 months ago by reubenmcgregor880
Answer: Limma analysis, controlling for Donors.
2
9 months ago by
Gordon Smyth38k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth38k wrote:

All you need is

design <- model.matrix( ~ Donor + MV_only, data = pData)


You don't even need a contrast.

Thanks,

I did do that originally and tried to carry on but was not sure about how to then setup the contrast matrix.

currently the design looks like:

> head(design, 10)
(Intercept) donor4 donor8 MV_onlyother
Ao3           1      0      0            1
Ao8           1      0      1            1
Ao4           1      1      0            1
AV3           1      0      0            1
AV4           1      1      0            1
AV8           1      0      1            1
LA3           1      0      0            1
LA4           1      1      0            1
LA8           1      0      1            1
LV3           1      0      0            1


so how would I set up the right contrast matrix as now obviously the contrast matrix I had is no longer applicable:

contrast.matrix <- makeContrasts(
mv_vs_other = MV-other,levels=design)


I would still like to compare MV vs other samples in "MV_only" column. My guess was

contrast.matrix <- makeContrasts(
mv_vs_other = MV_onlyother,levels=design)


But that does not seem to be right

1

As I said in my answer above, you don't need a contrast. Why do you think you need one? Just run

fit <- eBayes(fit)
topTable(fit, coef="MY_onlyother")


If you insist on computing a contrast explicitly, then your last contrast matrix above is correct and would give the same result.

I was getting confused with other analyses I had done with the same dataset where I needed a contrast and wanted to use the same structure. Thank you for your help yet again Gordon!

Gordon, a related question I have (though happy to put it on a seperate thread if the question is too unrelated) is how to remove the effect of donors with this design?

My attempt was:

donor <- pData\$donor
design_donor_effect <- model.matrix(~MV_only, data=pData)
y_batchremoved <- removeBatchEffect(y, batch=donor, design=design_donor_effect)


It seems to make sense to me but may actually be too long winded?

Yes, that's how it's done. Should only be using for plotting purposes though.