Question: Limma analysis, controlling for Donors.
0
gravatar for reubenmcgregor88
6 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:

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.

microarray limma R • 241 views
ADD COMMENTlink modified 6 months ago by Gordon Smyth37k • written 6 months ago by reubenmcgregor880
Answer: Limma analysis, controlling for Donors.
2
gravatar for Gordon Smyth
6 months ago by
Gordon Smyth37k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth37k wrote:

All you need is

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

You don't even need a contrast.

ADD COMMENTlink modified 6 months ago • written 6 months ago by Gordon Smyth37k

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 REPLYlink modified 6 months ago • written 6 months ago by reubenmcgregor880
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.

ADD REPLYlink modified 6 months ago • written 6 months ago by Gordon Smyth37k

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 REPLYlink written 6 months ago by reubenmcgregor880

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 REPLYlink written 6 months ago by reubenmcgregor880

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

ADD REPLYlink written 6 months ago by Gordon Smyth37k

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

ADD REPLYlink written 6 months ago by reubenmcgregor880
Please log in to add an answer.

Help
Access

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