[R]無監督聚類-ConsensusCluster & NMF

針對基因表達矩陣或者晶片資料，可以基於表達水平對資料進行無監督聚類分析。 完成這個功能的演算法有很多種，如非負矩陣分解（nonnegative matrix factorization），一致性聚類（consensus clustering）， 自組織對映（SOM clustering）等，文章

Metagenes and molecular pattern discovery using matrix factorization

介紹和比較了幾種演算法，文章總結NMF演算法更優，但其他幾種演算法針對不同特徵的資料也各有偏好，可以多瞭解一下演算法再自行選擇。

之前看一篇文章內容中涉及NMF演算法而程式碼用的是ConsensusClusterPlus包，未確認

演算法文章

就認為是用了同一個演算法，混淆了演算法和對應R包，謝謝評論指正。以此為鑑，演算法千千萬，還是要再仔細謹慎一點。

輸入資料格式

每列一個樣本，每行一個基因

- nonnegative matrix factorization

選擇cophenetic下降最快的前一點做較優rank（第二個引數）。

library（NMF）

library（doMPI） #rank設定4以上需要用這個包

#選定rank

# compute NMF for each method， rank設定為2：7 （2，3，4，5，6，7都試一遍一起比較）

res <- nmf（dat，2：7，nrun=10）

plot（res）

結果選rank=5

#選用較優rank再執行一次NMF

res <- nmf（dat，5，nrun=10）

coefmap（res）

consensusmap（res）

這裡的選擇用預設設定的熱圖來展示結果，詳細的還有很多引數可以最佳化結果輸出形式，這一步還沒有太理解，僅放圖做參考哈。

coefmap（res）

consensusmap（res）

NMF包中包含不同的演算法method（第三個引數），預設使用brunet，可以比較演算法擇優錄用

res。multi。method <- nmf（dat， 3， list（‘brunet’， ‘lee’， ‘ns’）， seed=222）

compare（res。multi。method）

- consensus clustering

根據consensus10的圖，選擇轉向平緩的拐點對應的x作為最優k（分組數）。所有直接輸出到 output/ （step4）資料夾下。 要注意調整step 2 中運用的基因數目，如果運用基因過多可能效果適得其反。

### step 1： ConsensusClusterPlus。

info <- as。data。frame（read。table（“data/exprSet_symbol。csv”，header = TRUE，sep = “，”， dec = “。”，na。strings = “NA”，stringsAsFactors=FALSE，check。names = FALSE））

#row。names（info） <- info［，1］

#info<-as。data。frame（info［，-1］）

#remove outlier if necessary

info<-info［，！colnames（info） %in% c（“sample20”）］

dat=info

### step 2： ConsensusClusterPlus。

d=dat

gene_e=8500 # 要提取的基因數，一般是70%左右

#按需選出決定分組的基因集

mads=apply（d，1，mad）

d=d［rev（order（mads））［1：gene_e］，］

### step 3： ConsensusClusterPlus。

d = sweep（d，1， apply（d，1，median，na。rm=T））

### step 4： ConsensusClusterPlus。

library（ConsensusClusterPlus）

#title=tempdir（）

#最大設定6個分組，重複1000次，結果輸出到資料夾output中

results = ConsensusClusterPlus（as。matrix（d），maxK=6，reps=1000，pItem=0。8，pFeature=1，

title=“output”，

clusterAlg=“hc”，distance=“pearson”，seed=12621，plot=“png”，writeTable=F）

### step 5： ConsensusClusterPlus。

#cat（sprintf（“\\graphicspath{{%s}}”， paste（gsub（“［\\］”，“/”，title），“/”，sep=“”）））

#cat（“\n”）

### step 6： ConsensusClusterPlus。

#consensusMatrix - the consensus matrix。

#For 。example， the top five rows and columns of results for k=2：

results［［3］］［［“consensusMatrix”］］［1：5，1：5］

#consensusTree - hclust object

results［［3］］［［“consensusTree”］］

#consensusClass - the sample classifications

clu3=as。data。frame（results［［3］］［［“consensusClass”］］）

colnames（clu3）=“cluster3”

#

clu4=as。data。frame（results［［4］］［［“consensusClass”］］）

colnames（clu4）=“cluster4”

NMF螢幕輸出

圖片輸出，分別是k=1到k=6的分組情況

落石圖表明按分組來說整體分組情況的穩定性，一般選取拐點（開始平緩點）來作為可行的k。 後續需要結合樣本的臨床資料，看這些分組是否有實際意義。

以這個圖為例，有一類是隻有一個數字的，可以檢視一下是哪一個樣本，看移除後再分析效果是否會好一點

#檢視每個分組的情況

table（clu4）

clu4$ID=row。names（clu4）

clu4［which（clu4$cluster4==4），］

關於分組結果驗證，一般來說，可以選取另外一個相同型別的資料集，然後選取同樣的基因集（即上文中根據資料集選出來的10k個基因），再次進行無監督聚類分組，看結果是否穩定，臨床資訊是否也具有同樣的規律。

先照常計算實驗資料集dat，得到基因集d，再計算驗證資料集dat2

需要把 chunk number 2 的部分改為

d=dat2［row。names（d），］

即可。