版塊導(dǎo)航: 正在加載中...

登錄注冊

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號能夠正常使用，請盡快對帳號進(jìn)行手機(jī)號驗證，感謝您的理解與支持！

24小時熱門版塊排行榜

返回列表

xpb2005

銀蟲 (小有名氣)

應(yīng)助: 0 (幼兒園)
金幣: 548
紅花: 1
帖子: 254
在線: 68小時
蟲號: 985644
注冊: 2010-03-30
專業(yè): 植物生理與生化

[求助] bioconductor中l(wèi)imma包，affylmGUI的使用已有1人參與

想分析芯片的數(shù)據(jù)，在R語言中用LIMMA包分析的，找了資料看了看，發(fā)現(xiàn)用LIMMA分析時要做一個矩陣設(shè)計，對此不太理解，希望前輩們給指導(dǎo)一下。下面
estrogen是一個2X2因素的基因芯片實驗數(shù)據(jù)。細(xì)胞為MCF7癌細(xì)胞，芯片為HGU95av2。實驗因子為estrogen存在與否，及其暴露時長（10或者48小時）。這個研究的目的就是為了調(diào)查哪些基因?qū)strogen做出應(yīng)答，哪些是早期應(yīng)答，哪些是晚期應(yīng)答。期待早期應(yīng)答與晚期應(yīng)答的基因分布在一個代謝途徑上，只是所處的位置不同。
使用limma來分析差異表達(dá)的基因，主要分幾步走：
1.    讀取數(shù)據(jù)
2.    預(yù)處理數(shù)據(jù)
3.    構(gòu)建實驗設(shè)計矩陣
4.    使用線性模型估計差異表達(dá)的倍數(shù)
5.    使用貝葉斯平滑標(biāo)準(zhǔn)差
6.    試用不同的參數(shù)來輸出差異表達(dá)基因結(jié)果。
下面我們就開始示例：
##首先讀入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算法來預(yù)處理數(shù)據(jù)
> eset <- rma(ab)
Background correcting
Normalizing
Calculating Expression
##構(gòu)建實驗設(shè)計矩陣，這一步非常關(guān)鍵。
##有很多辦法都可以構(gòu)建這個矩陣。
##在這組實驗數(shù)據(jù)中，我們關(guān)心的是早期和晚期細(xì)胞對estrogen的應(yīng)答，
  ##所以我們需要構(gòu)建兩組對比。
> 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
我的問題是：他做的矩陣設(shè)計，這一步我不太理解，矩陣中的1和0都是什么意思？貌似構(gòu)建矩陣時根據(jù)要比較的對象不一樣矩陣是有變化的，本例中關(guān)心的是早期細(xì)胞和晚期細(xì)胞對所添加物質(zhì)的應(yīng)答情況。我的實驗是對野生型（wt）和突變體（hm）分別用物質(zhì)A處理（用大寫的A表示）和不處理（用小寫的a表示）做的芯片，如果按照上面的程序走下來，那結(jié)果是野生型和突變體對物質(zhì)A的應(yīng)答情況。我用affylmGUI包做了三組對比，分別是awt-ahm;Awt-Ahm;awt-Ahm,我理解的是第一種比較是在不加物質(zhì)A處理時兩種材料之間的差異，第二種比較是在用物質(zhì)A處理后兩種材料的差異，第三種比較得出的差異基因是受材料和物質(zhì)A共同影響的，不知道我理解的對不對，
懇請高手給予矩陣設(shè)計方面的指點(diǎn)。
下面接著說例子（接上面的）：
##貝葉斯平滑
> 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值的總結(jié)，用于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值總結(jié)中的各項應(yīng)該是最小，四分之一分位，中位數(shù)，均值，四分之三分位，最大值，列出這些的目的是干什么的，意義何在，還有就是> results <- classifyTestsF(fit2, p.value=0.0001)中的p.value=0.0001是代表設(shè)定的p值的界限嗎？
接著上面的例子繼續(xù)：
##獲取差異表達(dá)基因的總結(jié)
> 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"
##生成差異表達(dá)基因的結(jié)果報表
> 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比較的結(jié)果嗎？
         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

我用上面的程序分析了我的數(shù)據(jù)，也用affylmGUI分析了我的數(shù)據(jù)，見下圖
圖中的第一行的M 和A  分別代表什么，我理解M帶表倍數(shù)變化？最后的B代表什么，是貝葉斯的什么嗎？
另外用affylmGUI包分析時，P-value的值怎么設(shè)定，或者說怎么知道默認(rèn)的alfa值是0.05還是0.01
我自己統(tǒng)計方面知識薄弱，也不懂編程，還望前輩們指導(dǎo)指導(dǎo)。不勝感激。

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

回復(fù)此樓

» 猜你喜歡

340求調(diào)劑已經(jīng)有3人回復(fù)
材料專碩找調(diào)劑已經(jīng)有4人回復(fù)
0854電子信息求調(diào)劑已經(jīng)有6人回復(fù)
0805 316求調(diào)劑已經(jīng)有4人回復(fù)
0854 考研調(diào)劑招生了！AI 方向已經(jīng)有17人回復(fù)
【考研調(diào)劑】化學(xué)專業(yè) 281分，一志愿四川大學(xué)，誠心求調(diào)劑已經(jīng)有15人回復(fù)
資源與環(huán)境調(diào)劑申請(333分) 已經(jīng)有7人回復(fù)
306求0703調(diào)劑一志愿華中師范已經(jīng)有11人回復(fù)
0703化學(xué)調(diào)劑，求導(dǎo)師收已經(jīng)有6人回復(fù)
求調(diào)劑已經(jīng)有5人回復(fù)

» 本主題相關(guān)商家推薦: (我也要在這里推廣)

1樓 2012-07-22 11:12:38

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

IsaacGK

鐵蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 99
帖子: 16
在線: 21.2小時
蟲號: 3921193
注冊: 2015-06-12
性別: GG
專業(yè): 生物信息學(xué)

【答案】應(yīng)助回帖

老師，請問您是如何用limma包去篩選差異表達(dá)基因的，現(xiàn)在真的很需要您的幫助

贊一下

回復(fù)此樓

2樓2015-06-13 21:18:32

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 xpb2005 的主題更新