Question

Limma analysis, controlling for Donors.

0

Entering edit mode

reubenmcgregor88 • 0

@reubenmcgregor88-13722

Last seen 3.8 years ago

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:

head(pData, n=15)

    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.

limma r microarray • 1.4k views

ADD COMMENT • link updated 5.9 years ago by Gordon Smyth 51k • written 5.9 years ago by reubenmcgregor88 • 0

score 2 · Accepted Answer · 2019-01-15

2

Entering edit mode

Gordon Smyth 51k

@gordon-smyth

Last seen 59 minutes ago

WEHI, Melbourne, Australia

All you need is

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

You don't even need a contrast.

ADD COMMENT • link 5.9 years ago Gordon Smyth 51k

0

Entering edit mode

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

ADD REPLY • link 5.9 years ago reubenmcgregor88 • 0

1

Entering edit mode

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.

ADD REPLY • link 5.9 years ago Gordon Smyth 51k

0

Entering edit mode

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!

ADD REPLY • link 5.9 years ago reubenmcgregor88 • 0

0

Entering edit mode

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?

ADD REPLY • link 5.9 years ago reubenmcgregor88 • 0

0

Entering edit mode

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

ADD REPLY • link 5.9 years ago Gordon Smyth 51k

0

Entering edit mode

Thanks, yes, I want to plot a few heat maps with DE genes.

ADD REPLY • link 5.9 years ago reubenmcgregor88 • 0