【在國際上洼冻，R軟件的應(yīng)用是數(shù)據(jù)分析的主流發(fā)展趨勢之一瞧毙，但我發(fā)現(xiàn)在國內(nèi)R軟件的使用遠不如SPSS、SAS等軟件那么流行秦爆。為推廣R軟件的使用序愚，本博客將陸續(xù)推出“R高級教程”系列專輯，希望對生命科學(xué)領(lǐng)域的科技工作者有少許幫助......】
通常來講等限，對于一般的統(tǒng)計分析爸吮，基于傻瓜式操作的SPSS（PASW）軟件已經(jīng)足夠，但在涉及個性化要求很高的復(fù)雜數(shù)據(jù)處理時望门，SPSS就開始顯得力不從心形娇，這時必須依賴功能更為強大的SAS等軟件。以前在自己的科研過程中分析數(shù)據(jù)多用SPSS筹误、SAS等埂软。在統(tǒng)計遺傳和基因組學(xué)領(lǐng)域，SAS可以處理很多問題，但與此同時勘畔，SAS實現(xiàn)復(fù)雜問題過于麻煩所灸，很多問題SAS也不是首選。后來開始運用R環(huán)境中的各種免費統(tǒng)計包炫七，特別是Bioconductor的系列分析包爬立，我發(fā)覺非常適合生命科學(xué)領(lǐng)域的研究者。R有很多優(yōu)點：
（1）免費万哪，不需要去尋找破解版侠驯，不用擔(dān)心版權(quán)問題，使用非常方便奕巍；
（2）功能非常強大吟策，單個包的功能比較有限，但多個包組合起來使用則功能無比強大的止，遠勝于SPSS檩坚、SAS等；
（3）源代碼開放诅福，稍作修改后就能滿足個性化的復(fù)雜統(tǒng)計分析匾委，滿足個性化需求是R的最大特點之一；
（4）程序閱讀容易氓润，再加上參考學(xué)習(xí)資料很多赂乐，上手比較容易，提高也不是很難咖气，根據(jù)個人經(jīng)驗挨措，要比SAS高級階段的進階容易許多；
（5）國際同行高度認同R崩溪，我發(fā)現(xiàn)很多專用軟件都開發(fā)了軟件的R版运嗜，今后R將是數(shù)據(jù)分析的主流發(fā)展方向。
R軟件的安裝悯舟、基本使用等初級教程就不談了，隨便在官方網(wǎng)站找個學(xué)習(xí)資料就搞定了砸民〉衷酰“R系列”專輯擬推出中級、高級分析教程岭参。今天推出基因表達譜芯片的聚類分析專題反惕。
本專題示例芯片數(shù)據(jù)來自GEO數(shù)據(jù)庫中檢索號為GSE11787(http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11787)的Affymetrix芯片的CEL文件，共6個CEL文件演侯，3個正常對照組姿染，3個HPS刺激組，為免疫器官脾臟的表達數(shù)據(jù)。
（一）原始數(shù)據(jù)的讀入悬赏、RNA降解評估和標(biāo)準(zhǔn)化
>?pd <- read.AnnotatedDataFrame("Target.txt",header=TRUE,row.names=1,as.is=TRUE)
>rawAffyData <- ReadAffy(filenames=pData(pd)$FileName, phenoData=pd)
> summary(exprs(rawAffyData))
> deg <- AffyRNAdeg(rawAffyData)
>plotAffyRNAdeg(deg, col=c(1,2,3,4,5,6))
> eset <- rma(rawAffyData)
> summary(exprs(eset))
> op <- par(mfrow=c(1,2))
>cols <- brewer.pal(6, "Set3")
>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ù)分析（如差異表達基因的檢測）的統(tǒng)計功效闽颇，往往需要進一步經(jīng)過Detection Call Filter和IQR filter等過濾（“基因芯片數(shù)據(jù)的特異性過濾與非特異性過濾”將在另一專題里專門討論）盾戴。
需要說明的是，常規(guī)做法是先篩選出差異表達基因兵多，然后只用差異表達基因進行聚類分析（本示例直接用了過濾后的數(shù)據(jù)集尖啡，聚類圖的效果稍差一點）。
（1）樣本聚類
>dd <- dist2(log2(exprs(eset2)))
>diag(dd) <- 0
>dd.row <- as.dendrogram(hclust(as.dist(dd)))
>row.ord <- order.dendrogram(dd.row)
>library("latticeExtra")
>legend <- list(top = list(fun = dendrogramGrob,
args = list(x = dd.row, side = "top")))
>lp <- levelplot(dd[row.ord, row.ord],
scales = list(x = list(rot = 90)),
xlab = "", ylab = "", legend = legend)
>plot(lp)
（2）二維聚類
>source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/my.colorFct.R")
>mydata<-exprs(eset2)
>mydatascale <- t(scale(t(mydata)))
>hr <- hclust(as.dist(1-cor(t(mydatascale), method="pearson")), method="complete")
>hc <- hclust(as.dist(1-cor(mydatascale, method="spearman")), method="complete")
>heatmap.2(mydata, Rowv=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)
上述聚類圖一般和論文里的聚類圖有點不同剩膘，聚類的模式不太直觀衅斩，你也可以用下面的語句進行更直觀的作圖：
>mycl <- cutree(hr, h=max(hr$height)/1.5);
>mycolhc <- sample(rainbow(256));mycolhc <- mycolhc[as.vector(mycl)]
>myc2 <- cutree(hc, h=max(hc$height)/1.5); mycolhr <- sample(rainbow(256)); mycolhr <- mycolhr[as.vector(myc2)]
>heatmap(mydatascale, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=my.colorFct(), scale="row", ColSideColors=mycolhr, RowSideColors=mycolhc)
（3）MantelCorrs聚類程序
>kmeans.result <- GetClusters(eset2, 500, 100)
>x=exprs(eset2)
>DistMatrices.result <- DistMatrices(x, kmeans.result$clusters)
>MantelCorrs.result <- MantelCorrs(DistMatrices.result$Dfull,DistMatrices.result$Dsubsets)
>permuted.pval <- PermutationTest(DistMatrices.result$Dfull, DistMatrices.result$Dsubsets,100, 16, 0.05)
>ClusterLists <- ClusterList(permuted.pval, kmeans.result$cluster.sizes,MantelCorrs.result)
>ClusterGenes <- ClusterGeneList(kmeans.result$clusters, ClusterLists$SignificantClusters,eset2)
>h=hclust(dist(MantelCorrs.result))
>plot(h)
【注：除Bioconductor圖標(biāo)外，所有圖片均為軟件實際運行所得】
本文引用地址：http://blog.sciencenet.cn/blog-295006-394794.html此文來自科學(xué)網(wǎng)朱猛進博客怠褐，轉(zhuǎn)載請注明出處畏梆。
