Comparing three groups of paired data with limma. How to make the paired design

Question

So I am using r, with the packages Bioconductor (oligo), (limma) to analyze some microarray data.

I am having trouble in the paired analysis.

So this is my phenotype data

ph@data ph@data index filename group WT1 WT WT1 WT WT2 WT WT2 WT WT3 WT WT3 WT WT4 WT WT4 WT LT1 LT LT1 LT LT2 LT LT2 LT LT3 LT LT3 LT LT4 LT LT4 LT TG1 TG TG1 TG TG2 TG TG2 TG TG3 TG TG3 TG TG4 TG TG4 TG

So to analyze I made this code:

colnames(ph@data)[2]="source"
ph@data
groups = ph@data$source
f = factor(groups,levels = c("WT","LT","TG"))
design = model.matrix(~ 0 + f) 
colnames(design)=c("WT","LT","TG")
design
data.fit = lmFit(normData,design)

> design
   WT LT TG
1   1  0  0
2   1  0  0
3   1  0  0
4   1  0  0
5   0  1  0
6   0  1  0
7   0  1  0
8   0  1  0
9   0  0  1
10  0  0  1
11  0  0  1
12  0  0  1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$f
[1] "contr.treatment"

Then I sue this to compare between my groups:

contrast.matrix = makeContrasts(LT-WT,TG-WT,LT-TG,levels=design) data.fit.con = contrasts.fit(data.fit,contrast.matrix) data.fit.eb = eBayes(data.fit.con)

So after all this I want to compare my groups:

ph@data[ ,2] = c("control","control","control","control","littermate","littermate","littermate","littermate","transgenic","transgenic","transgenic","transgenic")
colnames(ph@data)[2]="Treatment"
ph@data[ ,3] = c("B1","B2","B3","B4","B1","B2","B3","B4","B1","B2","B3","B4")
colnames(ph@data)[3]="BioRep"
ph@data```

groupsB = ph@data$BioRep 
groupsT = ph@data$Treatment
fb = factor(groupsB,levels=c("B1","B2","B3","B4"))
ft = factor(groupsT,levels=c("control","littermate","transgenic"))```

paired.design = (~ fb + ft)

So now my phenoData looks like this:

    index  Treatment BioRep
WT1    WT    control     B1
WT2    WT    control     B2
WT3    WT    control     B3
WT4    WT    control     B4
LT1    LT littermate     B1
LT2    LT littermate     B2
LT3    LT littermate     B3
LT4    LT littermate     B4
TG1    TG transgenic     B1
TG2    TG transgenic     B2
TG3    TG transgenic     B3
TG4    TG transgenic     B4```

B1 is the bioreplicates, then I have the groups that are wild type, littermate and transgenic

to compare my samples I am trying this

colnames(paired.design)=c("Intercept","B4vsB1","B3vsB1","B2vsB1","B4vsB2","B3vsB2","B4vsB3","littermatevscontrol","transgenicvscontrol")

But then I got this error: Error in `colnames<-`(`*tmp*`, value = c("Intercept", "WTvsLT", "WTvsTG", : attempt to set 'colnames' on an object with less than two dimensions

What I am doing wrong, and this is the right way to compare my data?

data.fit = lmFit(data.rma,paired.design)
data.fit$coefficients

score 1 · Answer 1 · answered Nov 28 '19 at 22:41

I made ft and fb since I do not have your data:

fb <- structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("B1", 
"B2", "B3", "B4"), class = "factor")

ft <- structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L), .Label = c("control", "littermate", "transgenic"), class = "factor")

There's a typo in this line:

paired.design = (~ fb + ft)
dim(paired.design)
NULL

It should be:

paired.design = model.matrix(~ fb + ft)
head(paired.design)
  (Intercept) fbB2 fbB3 fbB4 ftlittermate fttransgenic
1           1    0    0    0            0            0
2           1    1    0    0            0            0

If you look at your model.matrix, it has 6 columns which is not what you were trying to assign. So when you specify ~ fb + ft, the first level of the factor is chosen as the reference, and you find effect of other levels with respect to this. For example, for fb, B1 is the reference and you estimate the effect of B2,B3,B4.

To make the pairwise comparisons you want, for everything vs reference, you just need to specify the column itself, for example, for B4 vs B1, it will just be fbB4, because B4 is estimated with respect to B1. For others, you specify the "-" like before:

contrast.matrix = makeContrasts(fbB4,fbB3,fbB2,fbB4-fbB2,
fbB3-fbB2,fbB4-fbB3,ftlittermate,fttransgenic,levels=paired.design)

Now we can assign the col names as you wanted:

colnames(contrast.matrix) <-c("B4vsB1","B3vsB1","B2vsB1","B4vsB2","B3vsB2","B4vsB3",
"littermatevscontrol","transgenicvscontrol")

I simulate some data for normData to fit the model and contrasts:

set.seed(100)
normData = matrix(rpois(100*12,100),100,12)
data.fit = lmFit(normData,paired.design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)

Note, there is a warning:

Warning message:
In contrasts.fit(data.fit, contrast.matrix) :
  row names of contrasts don't match col names of coefficients

Because "(Intercept)" was renamed as Intercept.

And you can look at the coefficients under data.fit.con to see the fold changes of the comparisons you made

  Contrasts
           B4vsB1    B3vsB1      B2vsB1     B4vsB2     B3vsB2     B4vsB3
  [1,]  14.333333 11.000000   5.0000000   9.333333   6.000000   3.333333
  [2,]  -3.333333  2.000000   7.0000000 -10.333333  -5.000000  -5.333333
  [3,]   3.666667 -2.666667   6.3333333  -2.666667  -9.000000   6.333333
  [4,] -22.666667 -1.666667 -10.3333333 -12.333333   8.666667 -21.000000
  [5,]  -8.333333 -3.666667   8.3333333 -16.666667 -12.000000  -4.666667
  [6,]  15.333333 -5.666667  -0.3333333  15.666667  -5.333333  21.000000
      Contrasts
       littermatevscontrol transgenicvscontrol
  [1,]                1.25               -7.50
  [2,]                3.50               10.50
  [3,]               -3.75                2.75
  [4,]                7.75               14.00
  [5,]               16.75               -2.50
  [6,]                2.00                9.50

You can crosscheck with the original fitted coefficients to see that the correct contrast is made:

head(data.fit$coefficients)
     (Intercept)        fbB2      fbB3       fbB4 ftlittermate fttransgenic
[1,]    95.41667   5.0000000 11.000000  14.333333         1.25        -7.50
[2,]    94.33333   7.0000000  2.000000  -3.333333         3.50        10.50
[3,]    97.66667   6.3333333 -2.666667   3.666667        -3.75         2.75
[4,]    93.41667 -10.3333333 -1.666667 -22.666667         7.75        14.00
[5,]    98.91667   8.3333333 -3.666667  -8.333333        16.75        -2.50
[6,]    97.16667  -0.3333333 -5.666667  15.333333         2.00         9.50

Comparing three groups of paired data with limma. How to make the paired design

So now my phenoData looks like this:

to compare my samples I am trying this

1 Answers1