【R高級(jí)教程】專題一：表達(dá)譜芯片的聚類分析

1、【在國(guó)際上，R軟件的應(yīng)用是數(shù)據(jù)分析的主流發(fā)展趨勢(shì)之一，但我發(fā)現(xiàn)在國(guó)內(nèi)R軟件的使用遠(yuǎn)不如SPSS、SAS等軟件那么流行。為推廣R軟件的使用，本博客將陸續(xù)推出“R高級(jí)教程”系列專輯，希望對(duì)生命科學(xué)領(lǐng)域的科技工作者有少許幫助.】 通常來(lái)講，對(duì)于一般的統(tǒng)計(jì)分析，基于傻瓜式操作的SPSS（PASW）軟件已經(jīng)足夠，但在涉及個(gè)性化要求很高的復(fù)雜數(shù)據(jù)處理時(shí)，SPSS就開(kāi)始顯得力不從心，這時(shí)必須依賴功能更為強(qiáng)大的SAS等軟件。以前在自己的科研過(guò)程中分析數(shù)據(jù)多用SPSS、SAS等。在統(tǒng)計(jì)遺傳和基因組學(xué)領(lǐng)域，SAS可以處理很多問(wèn)題，但與此同時(shí)，SAS實(shí)現(xiàn)復(fù)雜問(wèn)題過(guò)于麻煩，很多問(wèn)題SAS也不是首選。后來(lái)開(kāi)始運(yùn)用R環(huán)境

2、中的各種免費(fèi)統(tǒng)計(jì)包，特別是Bioconductor的系列分析包，我發(fā)覺(jué)非常適合生命科學(xué)領(lǐng)域的研究者。R有很多優(yōu)點(diǎn)：（1）免費(fèi)，不需要去尋找破解版，不用擔(dān)心版權(quán)問(wèn)題，使用非常方便； （2）功能非常強(qiáng)大，單個(gè)包的功能比較有限，但多個(gè)包組合起來(lái)使用則功能無(wú)比強(qiáng)大，遠(yuǎn)勝于SPSS、SAS等； （3）源代碼開(kāi)放，稍作修改后就能滿足個(gè)性化的復(fù)雜統(tǒng)計(jì)分析，滿足個(gè)性化需求是R的最大特點(diǎn)之一； （4）程序閱讀容易，再加上參考學(xué)習(xí)資料很多，上手比較容易，提高也不是很難，根據(jù)個(gè)人經(jīng)驗(yàn)，要比SAS高級(jí)階段的進(jìn)階容易許多； （5）國(guó)際同行高度認(rèn)同R，我發(fā)現(xiàn)很多專用軟件都開(kāi)發(fā)了軟件的R版，今后R將是數(shù)據(jù)分析的主流發(fā)展方向

3、。 R軟件的安裝、基本使用等初級(jí)教程就不談了，隨便在官方網(wǎng)站找個(gè)學(xué)習(xí)資料就搞定了?！癛系列”專輯擬推出中級(jí)、高級(jí)分析教程。今天推出基因表達(dá)譜芯片的聚類分析專題。 本專題示例芯片數(shù)據(jù)來(lái)自GEO數(shù)據(jù)庫(kù)中檢索號(hào)為GSE11787的Affymetrix芯片的CEL文件，共6個(gè)CEL文件，3個(gè)正常對(duì)照組，3個(gè)HPS刺激組，為免疫器官脾臟的表達(dá)數(shù)據(jù)。（一）原始數(shù)據(jù)的讀入、RNA降解評(píng)估和標(biāo)準(zhǔn)化pd rawAffyData summary(exprs(rawAffyData) deg plotAffyRNAdeg(deg, col=c(1,2,3,4,5,6) eset summary(exprs(eset

4、) op cols boxplot(rawAffyData,col=cols,names=1:6, main = unnormalized.data)boxplot(data.frame(exprs(eset) ,names=1:6, main = normalization.data, col=blue, border=brown)par(op)（二）聚類分析 原始數(shù)據(jù)讀入，經(jīng)AffyBatch目標(biāo)轉(zhuǎn)成ExpressionSet目標(biāo)后，為提高后續(xù)分析（如差異表達(dá)基因的檢測(cè)）的統(tǒng)計(jì)功效，往往需要進(jìn)一步經(jīng)過(guò)Detection Call Filter和IQR filter等過(guò)濾（“基因芯片數(shù)據(jù)的特

5、異性過(guò)濾與非特異性過(guò)濾”將在另一專題里專門(mén)討論）。 需要說(shuō)明的是，常規(guī)做法是先篩選出差異表達(dá)基因，然后只用差異表達(dá)基因進(jìn)行聚類分析（本示例直接用了過(guò)濾后的數(shù)據(jù)集，聚類圖的效果稍差一點(diǎn)）。（1）樣本聚類dd diag(dd) dd.row row.ord library(latticeExtra)legend lp plot(lp)（2）二維聚類source(/tgirke/Documents/R_BioCond/My_R_Scripts/my.colorFct.R)mydatamydatascale hr hc heatmap.2(mydata, Ro

6、wv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=redgreen(75), scale=row, ColSideColors=heat.colors(length(hc$labels), RowSideColors=heat.colors(length(hr$labels), trace=none, key=T) 上述聚類圖一般和論文里的聚類圖有點(diǎn)不同，聚類的模式不太直觀，你也可以用下面的語(yǔ)句進(jìn)行更直觀的作圖：mycl mycolhc - sample(rainbow(256);mycolhc myc2 - cutree(hc, h=max(

7、hc$height)/1.5); mycolhr - sample(rainbow(256); mycolhr heatmap(mydatascale, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=my.colorFct(), scale=row, ColSideColors=mycolhr, RowSideColors=mycolhc)（3）MantelCorrs聚類程序kmeans.result x=exprs(eset2)DistMatrices.result MantelCorrs.result permuted.pval C

8、lusterLists ClusterGenes h=hclust(dist(MantelCorrs.result)plot(h)【注：除Bioconductor圖標(biāo)外，所有圖片均為軟件實(shí)際運(yùn)行所得】應(yīng)學(xué)生及個(gè)別博友的要求，盡管專業(yè)博文點(diǎn)擊率和反應(yīng)均很差，但在去San Diego參加PAG會(huì)議之前，還是抽時(shí)間給出【R高級(jí)教程】的第二專題。專題一給出了聚類分析的示例，本專題主要談在表達(dá)譜芯片分析中如何利用Bioconductor鑒定差異表達(dá)基因。 鑒定差異表達(dá)基因是表達(dá)譜芯片分析pipeline中必須的分析步驟。差異表達(dá)基因分析是根據(jù)表型協(xié)變量（分類變量）鑒定組間差異表達(dá)，它屬于監(jiān)督性分類的一種

9、。在鑒定差異表達(dá)基因以前，一般需要對(duì)表達(dá)值實(shí)施非特異性過(guò)濾（在機(jī)器學(xué)習(xí)框架下屬于非監(jiān)督性分類），因?yàn)檫m當(dāng)?shù)姆翘禺愋赃^(guò)濾可以提高差異表達(dá)基因的檢出率、甚至是功效。R分析差異表達(dá)基因的library有很多，但目前運(yùn)用最廣泛的Bioconductor包是limma。 本專題示例依然來(lái)自GEO數(shù)據(jù)庫(kù)中檢索號(hào)為GSE11787 的Affymetrix芯片的數(shù)據(jù)，數(shù)據(jù)介紹參閱專題一。library(limma)design - model.matrix( -1+factor(c(1,1,1, 2,2,2) 這個(gè)是根據(jù)芯片試驗(yàn)設(shè)計(jì)，對(duì)表型協(xié)變量的水平進(jìn)行design，比如本例中共有6張芯片，前3張為cont

10、rol對(duì)照組，后3張芯片為實(shí)驗(yàn)處理組，用1表示對(duì)照組，用2表示處理組。其他試驗(yàn)設(shè)計(jì)同理，比如2*2的因子設(shè)計(jì)試驗(yàn)，如果每個(gè)水平技術(shù)重復(fù)3次，那么可以表示為：design colnames(design) fit contrast.matrix fit fit2 fit2 resultssummary(results)vennCounts(results)vennDiagram(results)比較遺憾的是，目前l(fā)imma自帶的venn作圖函數(shù)不能做超過(guò)3維的高維venn圖，只能畫(huà)出3個(gè)圓圈的venn圖，即只能同時(shí)對(duì)三個(gè)coef進(jìn)行venn作圖。上面的venn圖只有一個(gè)coef，太簡(jiǎn)單了。下面是

11、一個(gè)由本實(shí)驗(yàn)室芯片數(shù)據(jù)得出的三個(gè)coef的venn圖例：heatDiagram(results,fit2$coef) 紅色為control中（與LPS相比）的高表達(dá)基因，綠色為control中（與LPS相比）的低表達(dá)基因，x軸的數(shù)字表示差異表達(dá)基因在eset2中所處的位置。xwrite.table(x, file=limma.xls, s=F, sep=t) 將結(jié)果寫(xiě)入limma.xls文件中，內(nèi)容包括AveExpr值（比較組間絕對(duì)值的平均差異值）、logFC值（差異倍數(shù)）、t值、P值、q值（即adj.P.Val值）和B值。一般logFC值、P值、q值和AveExpr值用來(lái)作為差異表達(dá)的判斷標(biāo)準(zhǔn)，比如差異倍數(shù)在2倍以上、絕對(duì)差異表達(dá)在10以上、P值小于0.01等。在Excel文件中，根據(jù)各項(xiàng)判斷標(biāo)準(zhǔn)排序，可以很容易地得到差異表達(dá)基因列表，這個(gè)列表可以用來(lái)進(jìn)行后續(xù)的分析，如GO注釋、基因網(wǎng)絡(luò)繪制等。 專題一中提到實(shí)際研究中，一般只用差異表達(dá)基因進(jìn)行聚類分析，在R中，根據(jù)差異表達(dá)結(jié)果過(guò)濾表達(dá)值很簡(jiǎn)單（具體的值可以依據(jù)芯片數(shù)據(jù)的實(shí)際情況設(shè)定，比如P值可以設(shè)寬松點(diǎn)0.05、logFC的絕對(duì)值