首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >用limma比较三组配对数据。如何进行配对设计

用limma比较三组配对数据。如何进行配对设计
EN

Stack Overflow用户
提问于 2019-11-29 04:29:21
回答 1查看 862关注 0票数 2

因此,我使用r和软件包Bioconductor (oligo),(limma)来分析一些微阵列数据。

我在配对分析中遇到了麻烦。

这是我的表型数据

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

所以为了分析,我写了这段代码:

代码语言:javascript
复制
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"

然后我在我的组之间进行比较:

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)

因此,在这之后,我想比较一下我的团队:

代码语言:javascript
复制
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)

所以现在我的phenoData看起来像这样:

代码语言:javascript
复制
    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是生物复制体,然后我有野生型,山羊型和转基因的群体

为了比较我的样本,我试着这样做

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

但是我得到了这个错误:Error in `colnames<-`(`*tmp*`, value = c("Intercept", "WTvsLT", "WTvsTG", : attempt to set 'colnames' on an object with less than two dimensions

我做错了什么,这是比较我的数据的正确方式吗?

代码语言:javascript
复制
data.fit = lmFit(data.rma,paired.design)
data.fit$coefficients
EN

回答 1

Stack Overflow用户

发布于 2019-11-29 06:41:01

我做了ft和fb,因为我没有你的数据:

代码语言:javascript
复制
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")

这一行有一个打字错误:

代码语言:javascript
复制
paired.design = (~ fb + ft)
dim(paired.design)
NULL

它应该是:

代码语言:javascript
复制
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

如果您查看您的model.matrix,它有6列,这不是您试图分配的列。因此,当您指定~ fb + ft时,将选择因子的第一个级别作为参考,您会发现与此相关的其他级别的影响。例如,对于fb,B1是参考,您可以估计B2、B3、B4的影响。

要进行您想要的成对比较,对于所有内容与引用,您只需指定列本身,例如,对于B4与B1,它将只是fbB4,因为B4是相对于B1进行估计的。对于其他情况,您可以像前面一样指定"-“:

代码语言:javascript
复制
contrast.matrix = makeContrasts(fbB4,fbB3,fbB2,fbB4-fbB2,
fbB3-fbB2,fbB4-fbB3,ftlittermate,fttransgenic,levels=paired.design)

现在,我们可以根据您的需要分配col名称:

代码语言:javascript
复制
colnames(contrast.matrix) <-c("B4vsB1","B3vsB1","B2vsB1","B4vsB2","B3vsB2","B4vsB3",
"littermatevscontrol","transgenicvscontrol")

我为normData模拟了一些数据,以符合模型并进行对比:

代码语言:javascript
复制
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)

注意,这里有一个警告:

代码语言:javascript
复制
Warning message:
In contrasts.fit(data.fit, contrast.matrix) :
  row names of contrasts don't match col names of coefficients

因为"( Intercept )“改名为Intercept。

您可以查看data.fit.con下的系数,以查看所做比较的折叠变化

代码语言:javascript
复制
  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

您可以与原始拟合系数交叉检查,以查看是否进行了正确的对比度:

代码语言:javascript
复制
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
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/59095788

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档