因此,我使用带有Bioconductor(oligo),(limma)软件包的r来分析一些微阵列数据。
我在配对分析中遇到麻烦。
这是我的表型数据
ph @ dataph@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
因此,我分析了这段代码:
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)
所以所有这些之后,我想比较我的组:
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)
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
我做错了,这是比较我的数据的正确方法吗?
data.fit = lmFit(data.rma,paired.design)
data.fit$coefficients
由于没有您的数据,所以我做了ft和fb:
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")
此行中有错别字:
paired.design = (~ fb + ft)
dim(paired.design)
NULL
应该是:
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估算的。对于其他用户,您可以像以前一样指定“-”:
contrast.matrix = makeContrasts(fbB4,fbB3,fbB2,fbB4-fbB2,
fbB3-fbB2,fbB4-fbB3,ftlittermate,fttransgenic,levels=paired.design)
现在我们可以根据需要分配列名称:
colnames(contrast.matrix) <-c("B4vsB1","B3vsB1","B2vsB1","B4vsB2","B3vsB2","B4vsB3",
"littermatevscontrol","transgenicvscontrol")
我为normData模拟一些数据以适合模型和对比:
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)
注意,有警告:
Warning message:
In contrasts.fit(data.fit, contrast.matrix) :
row names of contrasts don't match col names of coefficients
因为“(拦截)”被重命名为拦截。
并且您可以查看data.fit.con下的系数,以查看所做比较的倍数变化
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
您可以对原始拟合系数进行交叉检查,以确保形成正确的对比度:
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