A. 層次聚類分析案例(三)
之前的筆記:
聚類介紹: 點這里
層次聚類分析案例(一)
層次聚類分析案例(二)
獲取全基因組表達數據的能力是一項計算復雜度非常高的任務。由於人腦的局限性,是無法解決這個問題。但是,通過將基因分類進數量較少的類別後再進行分析,就能將基因數據加工到更易理解的水平。
聚類的目標是將一組基因進行劃分,使相似的基因落入同一個簇,同時不相似的基因落入不同的簇。這里需要考慮的關鍵問題是如何定義相似性,以及處理已分類基因。這里我們使用兩種基因類型的感光性來探索基因聚類問題。
准備工作
為了進行層次聚類,我們使用從實驗鼠身上採集的數據集。
第1步:收集和描述數據
該任務使用名為GSE4051_data和GSE4051_design的數據集。該數據集以標准格式存儲在名為GSE4051_data.csv和GSE4051_design.csv的CSV格式的文件中。數據獲取路徑: 在這里
GSE4051_data數據集包含29949行數據和39個變數。數值型變數如下:
GSE4051_design數據集包含39行數據和4個變數。數值型變數是:sidNum
非數值型變數是:sidChar;devStage;gType;
具體實施步驟以下為實現細節。
第2步:探索數據
RColorBrewer包是一個R包,可從 http://colorbrewer2.org 獲取,它提供地圖和其他圖形的彩色模板。
pvclust包用來實現非確定性的層次聚類分析。在層次聚類中,每個簇通過多尺度有放回抽樣計算p值。一個簇的p值在0~1之間。p值有兩種類型:近似無偏(approximately unbiased,AU)和有放回概率(bootstrap probability,BP)值。AU p值通過多尺度有放回採樣方法計算,經典的有放回採樣方法用來計算BP p值。AU p值相比BP p值存在優效性偏見。
xtable包可以生成LaTeX格式的表格。使用xtable可以將特定的R對象轉換成xtables。這些xtables能夠以LaTeX或HTML的格式輸出。
plyr包被用來進行分置合並(split-apply-combine,SAC)過程。它將一個大的問題切分成易處理的小塊,在每個小塊上進行操作,然後將所有小塊合並起來。
載入以下包:
讓我們探索並理解變數間的關系。從導入名為GSE4051_data.csv的CSV文件開始。我們將該文件數據存儲到GSE4051_data數據框中:
接下來,輸出GSE4051_data數據框的信息。str()函數返回GSE4051_data的結構信息。它簡略顯示了GSE4051_data數據框的內部結構。max.level指明了為了顯示網狀結構的最大等級。
結果如下:
下面,我們導入名為GSE4051_design.csv的CSV文件,將其數據保存到GSE4051_design數據框中:
輸出GSE4051_design數據框的內部結構。
結果如下:
第3步:轉換數據
為了便於後續的可視化階段,需要對每一行數據進行拉伸操作。這是由於在目前的要求下,不同基因表達之間存在絕對值的差距,因此需要對每一行數據進行拉伸。
中心化變數和創建z值是兩個常見的數據分析方法。scale函數中心化並拉伸數值型矩陣的列。
變換矩陣。傳入GSE4051_data數據框用t()函數進行數據框變換。
接下來,我們輸出GSE4051_data數據框的信息。通過設置give.attr=FALSE,次級結構的屬性不會被顯示。
結果如下:
round()函數用於舍入到最接近的整數。語法形式只有1種:Y = round(X),這里的X可以是數,向量,矩陣,輸出對應。
head()函數返回一個向量、矩陣、表、數據框或函數的頭部。GSE4051_data和trans_GSE4051_data數據框被當作對象傳入。rowMeans()函數計算每列的平均值。data.frame()函數創建數據框耦合變數集合,並且共享許多指標的性質:
結果如下:
第4步:訓練模型
接下來是訓練模型。第一步是計算距離矩陣。dist()函數用來計算並返回距離矩陣,可以使用特定的距離度量方法來計算數據矩陣中各行間的距離。這里可使用的距離度量方法有歐式距離、最大距離、曼哈頓距離、堪培拉距離、二進制距離,或閔可夫斯基距離。這里使用歐式距離。歐式距離計算兩個向量間的距離公式為sqrt(sum((x_i-y_i)^2))。轉換後的trans_GSE4051_data數據框被用來計算距離。結果存儲在pair_dist_GSE4051_data數據框中。
接下來,使用interaction()函數計算並返回gType、devStage變數間相互作用的無序因子。無序因子的結果連同GSE4051_design數據框一同被傳入with()函數。該函數計算產生一個新的因子代表gType、devStage變數的相互作用:
summary()函數用來生成GSE4051_design$group數據框的結果總結:
結果如下:
下面,使用多種不同的聯合類型計算層次聚類。
使用hclust()函數對n個不同對象進行聚類分析。第一個階段,每個對象被指派給自己的簇。演算法在每個階段迭代聚合兩個最相似的簇。持續該過程直到只剩一個單獨的簇。hclust()函數要求我們以距離矩陣的形式提供數據。pair_dist_GSE4051_data數據框被傳入。
在第一個例子中使用single聚類方法:
結果如下:
在第二個例子中使用complete聚集方法。
調用pr.hc.complete的結果是顯示所使用的聚集方法、距離計算方法和對象數量:
結果如下:
在第三個例子中使用average聚類方法:
調用pr.hc.complete的結果是顯示所使用的聚集方法、距離計算方法和對象數量:
結果如下:
在第四個例子中使用ward聚類方法:
pr.hc.ward的調用結果是顯示所使用的聚集方法、距離計算方法和對象數量:
結果如下:
plot()函數是繪制R對象的通用函數。
第一次調用plot()函數,傳遞pr.hc.single數據框作為輸入對象:
結果如下:
第二次調用plot()函數,傳入pr.hc.complete數據框作為輸入對象:
結果如下:
第三次調用plot()函數,傳入pr.hc.average數據框作為輸入對象:
結果如下:
第四次調用plot()函數,傳入pr.hc.ward數據框作為輸入對象:
結果如下:
第5步:繪制模型
plot()函數是繪制R對象的通用函數。這里,plot()函數用來繪制系統樹圖。
rect.hclust()函數強調不同的簇,並在系統樹圖的枝幹上繪制長方形。系統樹圖首先在某個等級上被剪切,之後在選定的枝幹上繪制長方形。
RColorBrewer使用從 http://colorbrewer2.org 獲得的包來選擇繪制R圖像的顏色模板。
顏色分為三組:
最重要的一個RColorBrewer函數是brewer.pal()。通過向該函數傳入顏色的數量和配色的名字,可以從display.brewer.all()函數中選擇一個配色方案。
在第一個例子中,pr.hc.single作為一個對象傳入plot()函數:
結果如下:
下面創建熱度圖,使用single聚集方法。heatmap()函數默認使用euclidean聚集方法:
結果如下:
在第二例子中,pr.hc.complete作為對象傳入plot()函數:
結果如下:
下面使用complete聚集方法創建熱度圖:
結果如下:
在第三個例子中,pr.hc.average作為對象傳入plot()函數:
結果如下:
下面創建average聚集方法的熱度圖:
結果如下:
在第四個例子中,pr.hc.ward作為對象傳入plot()函數:
結果如下:
下面繪制ward聚集方法的熱度圖:
結果如下:
B. 在進行系統聚類分析時,不同的類間距離計算方法有何區別
聚類分析有兩種主要計算方法,分別是凝聚層次聚類(Agglomerative hierarchical method)和K均值聚類(K-Means)。
一、層次聚類
層次聚類又稱為系統聚類,首先要定義樣本之間的距離關系,距離較近的歸為一類,較遠的則屬於不同的類。可用於定義「距離」的統計量包括了歐氏距離 (euclidean)、馬氏距離(manhattan)、 兩項距離(binary)、明氏距離(minkowski)。還包括相關系數和夾角餘弦。
層次聚類首先將每個樣本單獨作為一類,然後將不同類之間距離最近的進行合並,合並後重新計算類間距離。這個過程一直持續到將所有樣本歸為一類為止。在計算類間距離時則有六種不同的方法,分別是最短距離法、最長距離法、類平均法、重心法、中間距離法、離差平方和法。
下面我們用iris數據集來進行聚類分析,在R語言中所用到的函數為hclust。首先提取iris數據中的4個數值變數,然後計算其歐氏距離矩陣。然後將矩陣繪制熱圖,從圖中可以看到顏色越深表示樣本間距離越近,大致上可以區分出三到四個區塊,其樣本之間比較接近。
data=iris[,-5]
dist.e=dist(data,method='euclidean')
heatmap(as.matrix(dist.e),labRow = F, labCol = F)
X
然後使用hclust函數建立聚類模型,結果存在model1變數中,其中ward參數是將類間距離計算方法設置為離差平方和法。使用plot(model1)可以繪制出聚類樹圖。如果我們希望將類別設為3類,可以使用cutree函數提取每個樣本所屬的類別。
model1=hclust(dist.e,method='ward')
result=cutree(model1,k=3) 為了顯示聚類的效果,我們可以結合多維標度和聚類的結果。先將數據用MDS進行降維,然後以不同的的形狀表示原本的分類,用不同的顏色來表示聚類的結果。可以看到setose品種聚類很成功,但有一些virginica品種的花被錯誤和virginica品種聚類到一起。
C. 層次聚類分析案例(一)
關於聚類分析的介紹,可參見本人之前的筆記: 聚類分析
案例一:世界銀行樣本數據集
創建世界銀行的一個主要目標是對抗和消除貧困。在這個不斷發展的世界中,世界銀行持續的發展並精細地調整它的政策,已經幫助這個機構逐漸實現了消除貧困的目標。消除貧困的成果以下指標的改進衡量,這些指標包括健康、教育、衛生、基礎設施以及其他需要用於改進窮人生活的服務。與此同時,發展成果必須保證以一種環保的、全社會的、經濟可持續的方式達成。
准備工作
為了進行層次聚類,我們需要使用從世界銀行收集的數據集。
第1步:收集和描述數據
該任務使用名為WBClust2013的數據集。該數據以標准格式存儲在名為WBClust2013.csv的CSV格式的文件中。其有80行數據和14個變數。 點我獲取數據
第一列Country為非數值型變數,其他列均為數值型變數。
第2步:探索數據
讓我們探索數據並理解變數間的關系。我們通過導入名為WBClust2013.csv的CSV文件開始。存儲數據到wbclust數據框中:
下一步輸出wbclust數據框,head()函數返回wbclust數據框。wbclust數據框作為一個輸入參數傳入:
結果如下:
第3步:轉換數據
中心化變數和創建z值是兩個常見的用於歸一化數據的數據分析手段。上面提到的數值型變數需要創建z值。scale()函數是一個通用的函數,其默認方法中心化並比例縮放一個數值化矩陣的列。數據框wbclust被傳給該比例函數。只有數據框中數值化的變數會被縮放。結果存儲在wbnorm數據框中。
結果如下:
所有的數據框都有rownames屬性。rownames()函數用來獲取或設置矩陣類變數的行名或列名。數據框wbclust以及第一列被傳遞給rownames()函數。
調用rownames(wbnorm)方法顯示第一列的數值。結果如下:
第4步:訓練並評估模型效果
下一步是訓練模型。首先使用dist()函數計算距離矩陣。使用特定的距離度量方法計算數據矩陣行間的距離。使用的距離度量可以是歐式距離、最大距離、曼哈頓距離、堪培拉距離、二進制距離,或閔可夫斯基距離。這里的距離度量使用歐式距離。使用歐式距離計算兩個向量間的距離為sqrt(sum((x_i-y_i)^2))。結果被存儲在一個新的數據框dist1中。
下一步是使用Ward方法進行聚類。hclust()函數對一組不同的n個對象進行聚類分析。第一階段,每個對象被指派給它自己的簇。之後每個階段,演算法迭代聚合兩個最相似的簇。這個過程不斷持續直到只剩一個簇。hclust()函數要求我們以距離矩陣的形式提供數據。dist1數據框被作為輸入傳入。默認使用全鏈接演算法。此外還可以使用不同的聚集方法,包括ward.D、ward.D2、single、complete和average。
輸入clust1命令可顯示所使用的聚集方法,計算距離的方法,以及數據對象的數量。結果如下:
第5步:繪制模型
plot()函數是一個通用的繪制R語言對象的函數。這里plot()函數用來繪制系統樹圖:
結果如下:
rect.hclust()函數強調不同的簇,並在系統樹圖的枝幹處繪制長方形。系統樹圖首先在某個等級上被剪切,之後在選定的枝幹上繪制長方形。
clust1對象以及需要形成的簇的數量作為輸入變數傳入函數。
結果如下:
cuts()函數基於期望的簇數量或者切割高度將樹中的元素切割到不同的簇中。這里,clust1對象以及需要形成的簇的數量作為輸入變數傳入函數。
結果如下:
得到每個簇的國家列表:
結果如下: