24小时热门版块排行榜

返回列表

xpb2005

银虫 (小有名气)

应助: 0 (幼儿园)
金币: 548
红花: 1
帖子: 254
在线: 68小时
虫号: 985644
注册: 2010-03-30
专业: 植物生理与生化

[求助] bioconductor中limma包，affylmGUI的使用已有1人参与

想分析芯片的数据，在R语言中用LIMMA包分析的，找了资料看了看，发现用LIMMA分析时要做一个矩阵设计，对此不太理解，希望前辈们给指导一下。下面
estrogen是一个2X2因素的基因芯片实验数据。细胞为MCF7癌细胞，芯片为HGU95av2。实验因子为estrogen存在与否，及其暴露时长（10或者48小时）。这个研究的目的就是为了调查哪些基因对estrogen做出应答，哪些是早期应答，哪些是晚期应答。期待早期应答与晚期应答的基因分布在一个代谢途径上，只是所处的位置不同。
使用limma来分析差异表达的基因，主要分几步走：
1.    读取数据
2.    预处理数据
3.    构建实验设计矩阵
4.    使用线性模型估计差异表达的倍数
5.    使用贝叶斯平滑标准差
6.    试用不同的参数来输出差异表达基因结果。
下面我们就开始示例：
##首先读入limma库
> library(limma)
>library(affy)
Loading required package: Biobase
Welcome to Bioconductor
  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation("Biobase"

' and for packages 'citation(pkgname)'.
##载入CDF文件库
> library(hgu95av2cdf)
> workingdir<-getwd()
> workingdir<-paste(workingdir,"/4R/estrogen",sep=""

> setwd(workingdir)
> targets<-readTargets("estrogen.txt",sep="",row.names="filename"

> targets
               filename estrogen time.h
low10-1.cel low10-1.cel absent    10
low10-2.cel low10-2.cel absent    10
high10-1.cel high10-1.cel  present    10
high10-2.cel high10-2.cel  present    10
low48-1.cel low48-1.cel absent    48
low48-2.cel low48-2.cel absent    48
high48-1.cel high48-1.cel  present    48
high48-2.cel high48-2.cel  present    48
> ab <- ReadAffy(filenames=targets$filename)
##使用rma算法来预处理数据
> eset <- rma(ab)
Background correcting
Normalizing
Calculating Expression
##构建实验设计矩阵，这一步非常关键。
##有很多办法都可以构建这个矩阵。
##在这组实验数据中，我们关心的是早期和晚期细胞对estrogen的应答，
  ##所以我们需要构建两组对比。
> treatments <- factor(c(rep(1:4,each=2)),labels=c("e10","E10","e48","E48"

)
> treatments
[1] e10 e10 E10 E10 e48 e48 E48 E48
Levels: e10 E10 e48 E48
> contrasts(treatments) <- cbind(Time=c(0,0,1,1),E10=c(0,1,0,0),E48=c(0,0,0,1))
> treatments
[1] e10 e10 E10 E10 e48 e48 E48 E48
attr(,"contrasts"

Time E10 E48
e10 0 0 0
E10 0 1 0
e48 1 0 0
E48 1 0 1
Levels: e10 E10 e48 E48
> design<-model.matrix(~treatments)
> design
  (Intercept) treatmentsTime treatmentsE10 treatmentsE48
1          1             0          0          0
2          1             0          0          0
3          1             0          1          0
4          1             0          1          0
5          1             1          0          0
6          1             1          0          0
7          1             1          0          1
8          1             1          0          1
attr(,"assign"

[1] 0 1 1 1
attr(,"contrasts"

attr(,"contrasts"

$treatments
Time E10 E48
e10 0 0 0
E10 0 1 0
e48 1 0 0
E48 1 0 1
> colnames(design) <- c("Intercept","Time","E10","E48"

> design
  Intercept Time E10 E48
1       1 0 0 0
2       1 0 0 0
3       1 0 1 0
4       1 0 1 0
5       1 1 0 0
6       1 1 0 0
7       1 1 0 1
8       1 1 0 1
attr(,"assign"

[1] 0 1 1 1
attr(,"contrasts"

attr(,"contrasts"

$treatments
Time E10 E48
e10 0 0 0
E10 0 1 0
e48 1 0 0
E48 1 0 1
##拟合线性模型
> fit <- lmFit(eset, design)
> dim(fit)
[1] 12625    4
> cont.matrix <- cbind(E10=c(0,0,1,0),E48=c(0,0,0,1))
> cont.matrix
   E10 E48
[1,] 0 0
[2,] 0 0
[3,] 1 0
[4,] 0 1
> fit2 <- contrasts.fit(fit,cont.matrix)
> dim(fit2)
[1] 12625    2
我的问题是：他做的矩阵设计，这一步我不太理解，矩阵中的1和0都是什么意思？貌似构建矩阵时根据要比较的对象不一样矩阵是有变化的，本例中关心的是早期细胞和晚期细胞对所添加物质的应答情况。我的实验是对野生型（wt）和突变体（hm）分别用物质A处理（用大写的A表示）和不处理（用小写的a表示）做的芯片，如果按照上面的程序走下来，那结果是野生型和突变体对物质A的应答情况。我用affylmGUI包做了三组对比，分别是awt-ahm;Awt-Ahm;awt-Ahm,我理解的是第一种比较是在不加物质A处理时两种材料之间的差异，第二种比较是在用物质A处理后两种材料的差异，第三种比较得出的差异基因是受材料和物质A共同影响的，不知道我理解的对不对，
恳请高手给予矩阵设计方面的指点。
下面接着说例子（接上面的）：
##贝叶斯平滑
> fit2 <- eBayes(fit2)
> dim(fit2)
[1] 12625    2
##选择测试方法，global或者nestedF。
> results.global <- decideTests(fit2, method="global"

> results.nestedF<- decideTests(fit2, method="nestedF"

##去除嵌入探针
> i <- grep("AFFX",featureNames(eset))
##获取p值的总结，用于p值背景去除。
  > summary(fit2$F.p.value)
   Min. 1st Qu. Median    Mean 3rd Qu.    Max.
0.0001391 0.1726000 0.3562000 0.4207000 0.6825000 0.9925000
> results <- classifyTestsF(fit2, p.value=0.0001)
我的问题是p值总结中的各项应该是最小，四分之一分位，中位数，均值，四分之三分位，最大值，列出这些的目的是干什么的，意义何在，还有就是> results <- classifyTestsF(fit2, p.value=0.0001)中的p.value=0.0001是代表设定的p值的界限吗？
接着上面的例子继续：
##获取差异表达基因的总结
> summary(results) #  问题是这个表格看不太懂？
   E10 E48
-1 40 76
0  12469 12410
1 116 139
> table(E10=results[,1],E48=results[,2]) #表格看不太懂？
E48
E10    -1    0    1
  -1 29 11    0
  0    47 12370 52
  1    0 29 87
##绘制文氏图
> vennDiagram(results, include="up"
> x11()
> vennDiagram(results, include="down"
##生成差异表达基因的结果报表
> options(digits=3) ？？？？digits=3是什么意思？
> topTable(fit2,coef="E10",n=20) ？？？？不理解coef=E10是什么意思，是代表得到的差异基因是e10和E10的比较吗？
         ID logFC AveExpr    t  P.Value adj.P.Val B
9735 39642_at  2.94 7.88  23.7 4.74e-09  3.13e-05 9.97
12472 910_at  3.11 9.66  23.6 4.96e-09  3.13e-05 9.94
1814 31798_at  2.80 12.12  16.4 1.03e-07  3.51e-04 7.98
11509  41400_at  2.38 10.04  16.2 1.11e-07  3.51e-04 7.92
10214  40117_at  2.56 9.68  15.7 1.47e-07  3.58e-04 7.71
953    1854_at  2.51 8.53  15.2 1.95e-07  3.58e-04 7.49
9848 39755_at  1.68 12.13  15.1 2.05e-07  3.58e-04 7.45
922 1824_s_at  1.91 9.24  14.9 2.27e-07  3.58e-04 7.37
140 1126_s_at  1.78 6.88  13.8 4.12e-07  5.78e-04 6.89
580    1536_at  2.66 5.94  13.3 5.80e-07  7.32e-04 6.61
12542 981_at  1.82 7.78  13.1 6.46e-07  7.42e-04 6.52
3283 33252_at  1.74 8.00  12.6 8.86e-07  9.20e-04 6.25
546    1505_at  2.40 8.76  12.5 9.48e-07  9.20e-04 6.19
4405 34363_at -1.75 5.55 -12.2 1.14e-06  1.03e-03 6.03
985 1884_s_at  2.80 9.03  12.1 1.26e-06  1.06e-03 5.95
6194 36134_at  2.49 8.28  11.8 1.50e-06  1.19e-03 5.79
7557 37485_at  1.61 6.67  11.4 1.99e-06  1.48e-03 5.55
1244    239_at  1.57 11.25  10.4 4.07e-06  2.66e-03 4.90
8195 38116_at  2.32 9.51  10.4 4.09e-06  2.66e-03 4.90
10634  40533_at  1.26 8.47  10.4 4.21e-06  2.66e-03 4.87
> topTable(fit2,coef="E48",n=20) ？？？这里的coef=E48是代表得到的差异基因是e48和E48比较的结果吗？
         ID logFC AveExpr    t  P.Value adj.P.Val    B
12472 910_at  3.86 9.66  29.2 8.27e-10  1.04e-05 11.61
1814 31798_at  3.60 12.12  21.0 1.28e-08  7.63e-05  9.89
953    1854_at  3.34 8.53  20.2 1.81e-08  7.63e-05  9.64
8195 38116_at  3.76 9.51  16.9 8.12e-08  2.51e-04  8.48
8143 38065_at  2.99 9.10  16.2 1.12e-07  2.51e-04  8.21
9848 39755_at  1.77 12.13  15.8 1.36e-07  2.51e-04  8.05
642    1592_at  2.30 8.31  15.8 1.39e-07  2.51e-04  8.03
11509  41400_at  2.24 10.04  15.3 1.81e-07  2.75e-04  7.81
366 33730_at -2.04 8.57 -15.1 1.96e-07  2.75e-04  7.74
732    1651_at  2.97 10.50  14.8 2.39e-07  3.02e-04  7.57
8495 38414_at  2.02 9.46  14.6 2.66e-07  3.05e-04  7.48
1049 1943_at  2.19 7.60  14.0 3.72e-07  3.69e-04  7.18
10214  40117_at  2.28 9.68  14.0 3.80e-07  3.69e-04  7.17
10634  40533_at  1.64 8.47  13.5 4.94e-07  4.45e-04  6.93
9735 39642_at  1.61 7.88  13.0 6.71e-07  5.18e-04  6.65
4898 34851_at  1.96 9.96  12.8 7.51e-07  5.18e-04  6.55
922 1824_s_at  1.64 9.24  12.8 7.95e-07  5.18e-04  6.50
6053 35995_at  2.76 8.87  12.7 8.32e-07  5.18e-04  6.46
12455 893_at  1.54 10.95  12.7 8.43e-07  5.18e-04  6.45
10175  40079_at -2.41 8.23 -12.6 8.62e-07  5.18e-04  6.42

我用上面的程序分析了我的数据，也用affylmGUI分析了我的数据，见下图
图中的第一行的M 和A  分别代表什么，我理解M带表倍数变化？最后的B代表什么，是贝叶斯的什么吗？
另外用affylmGUI包分析时，P-value的值怎么设定，或者说怎么知道默认的alfa值是0.05还是0.01
我自己统计方面知识薄弱，也不懂编程，还望前辈们指导指导。不胜感激。

[ Last edited by xpb2005 on 2012-7-22 at 11:23 ]