最近在進(jìn)行ATAC和RNA-Seq的聯(lián)合分析抛杨，由于處理的材料是時(shí)間相關(guān)的驼抹，所以time-course也是一個(gè)可以分析的點(diǎn)荸百。在一位劉潛老哥的幫助下咕幻，我找到了一篇靠譜熊 轉(zhuǎn)錄組時(shí)間序列數(shù)據(jù)處理 的文章，里面提到了mfuzz這個(gè)包糯俗。里面的軟聚類的思想非常符合我的預(yù)期尿褪，然后就決定拿這個(gè)包進(jìn)行我time-course的分析。為了更好的分析得湘，我決定先翻譯下這個(gè)包杖玲，了解下這個(gè)包的大致思想。

1 Overview

這部分是我偷懶的隨便寫的淘正。摆马。臼闻。。囤采。述呐。

感覺time-course的方法一般就是聚類。常見的聚類分為三種斑唬，分別是層次聚類（Hierarchical Clustering）市埋、硬聚類（hard clustering）、軟聚類（soft clustering）恕刘。層次聚類好像一般就是像熱圖那種缤谎。看了pheatmap的文檔褐着，感覺pheatmap就是層次聚類坷澡，當(dāng)然你可以設(shè)置k_means，變成硬聚類含蓉。硬聚類常見的就是k-menas频敛。軟聚類就是我們這會(huì)要用到的這個(gè)包的核心思路。

2 Installation requirements

見Bioconductor的安裝方法馅扣。

3 Data pre-processing

數(shù)據(jù)集是來(lái)源于酵母細(xì)胞循環(huán)表達(dá)數(shù)據(jù)斟赚。6178個(gè)基因，橫跨160分鐘的17個(gè)時(shí)間點(diǎn)差油。用的是芯片數(shù)據(jù)拗军。

> head(yeast@assayData$exprs)
          cdc28_0 cdc28_10 cdc28_20 cdc28_30 cdc28_40 cdc28_50 cdc28_60 cdc28_70 cdc28_80 cdc28_90 cdc28_100 cdc28_110 cdc28_120 cdc28_130 cdc28_140 cdc28_150 cdc28_160
YDR132C      0.19     0.30    -0.29     0.29    -0.31     0.23     0.20    -0.08     0.19    -0.51      0.00     -0.31      0.11     -0.02      0.20      0.36     -0.54
YMR012W     -0.15    -0.15    -0.04    -0.28    -0.39     0.03     0.22     0.04    -0.15     0.37      0.47     -0.10     -0.09        NA     -0.04      0.07      0.19
YLR214W      0.38     0.30    -0.68    -0.52    -0.43    -0.13    -0.17     0.26    -0.03    -0.34     -0.01     -0.20      0.10        NA      0.45      0.40      0.63
YLR116W      0.17     0.06    -0.21     0.19     0.33     0.44     0.46     0.38    -0.15    -0.03      0.04     -0.42     -0.15      0.02        NA     -0.51     -0.61
YDR203W      0.85    -0.10    -0.56    -0.31    -0.43     0.00    -0.34     0.17     0.40    -0.37      0.15      0.24      0.24      0.17     -0.12     -0.02      0.02
YEL059C-A    0.45     0.20     0.06     0.10    -0.21    -0.08    -0.27    -0.01    -0.29     0.41     -0.08     -0.22     -0.27        NA     -0.30      0.25      0.26

3.1 Missing value

第一步，去除那些有超過(guò)25%數(shù)據(jù)缺失的基因蓄喇。注意這些數(shù)據(jù)缺失值應(yīng)該是設(shè)為NA发侵。

yeast.r <- filter.NA(yeast, thres=0.25)
49 genes excluded.

這里就如上面的數(shù)據(jù)一樣，一行即一個(gè)基因有16個(gè)時(shí)間點(diǎn)的數(shù)據(jù)妆偏，如果16個(gè)時(shí)間點(diǎn)里面有25%刃鳄，即4個(gè)時(shí)間點(diǎn)都是NA，則剔除這個(gè)基因钱骂。

> nrow(yeast)
Features 
    3000 
> nrow(yeast.r)
Features 
    2951 

Fuzzy c-means就像其他聚類算法一樣叔锐，其并不允許有缺失值的存在。所以我們會(huì)對(duì)剩下那些缺失值（16個(gè)數(shù)據(jù)點(diǎn)里面就缺了1個(gè)2個(gè)那種）進(jìn)行填充罐柳。用的是對(duì)應(yīng)基因的平均表達(dá)值掌腰。

對(duì)于RNA-Seq來(lái)說(shuō)，你可以加上一些pseudocount张吉，比如0.01。

 yeast.f <- fill.NA(yeast.r,mode="mean")

當(dāng)然催植，你也可以用（weighted） k-nearest neighbour method肮蛹。（mode='knn'/'wknn'）勺择。這些方法相比較而言比上面這種簡(jiǎn)單的方法要好，但需要耗費(fèi)更多的算力伦忠。

3.2 Filtering

許多已經(jīng)出版的聚類分析包含過(guò)濾的步驟省核，從而來(lái)去除那些表達(dá)相對(duì)比較低的，或者表達(dá)不怎么變化的昆码。通常來(lái)說(shuō)气忠，比較受歡迎的就是樣本的標(biāo)準(zhǔn)差作為閾值。

tmp <- filter.std(yeast.f,min.std=0)

然而在基因低表達(dá)到高表達(dá)的過(guò)程中赋咽，變化是非常平緩的旧噪。所以給定閾值篩選并不一定是可靠的，可能是非常武斷脓匿。因?yàn)楝F(xiàn)在并沒有很多有說(shuō)服力的篩選手段淘钟，所以我們還是避免對(duì)基因數(shù)據(jù)做提前的篩選。這可以避免損失一些有生物學(xué)重大意義的基因陪毡。

比如1,2,4,10,12,13,15米母。看起來(lái)變化很大毡琉，但方差可能并不如你想象中的那么大铁瞒。

Standardisation

由于聚類是在歐幾里德空間中進(jìn)行的，因此基因的表達(dá)值被標(biāo)準(zhǔn)化為平均值為零桅滋，標(biāo)準(zhǔn)差為1慧耍。該步驟確保了在歐幾里得空間中具有相似表達(dá)模式的基因是相互接近的。

yeast.s <- standardise(yeast.f)

重要的是虱歪，Mfuzz認(rèn)為輸入的表達(dá)數(shù)據(jù)是完全經(jīng)過(guò)前期數(shù)據(jù)標(biāo)準(zhǔn)化的蜂绎。standardise 并不能代替標(biāo)準(zhǔn)化步驟。注意差異：標(biāo)準(zhǔn)化是為了讓不同的樣品間可以比較笋鄙，而Mufzz中standardisation則是讓轉(zhuǎn)錄本或者基因間可以比較师枣。

4 Soft clustering of gene expression data

聚類可以用來(lái)解釋基因表達(dá)的調(diào)控機(jī)制。眾所周知的萧落，基因的表達(dá)并不是開和關(guān)的践美，而是一個(gè)逐漸變化的過(guò)程。一個(gè)聚類算法應(yīng)該展現(xiàn)出一個(gè)基因有多么的符合dominant cluster pattern找岖。軟聚類應(yīng)該是一個(gè)非常好的方法陨倡，因?yàn)槠淇梢岳胢embership 衡量一個(gè)基因 i跟cluster j的關(guān)系。

其實(shí)就是說(shuō)基因A跟每個(gè)cluster都有關(guān)系许布，無(wú)非是membership score的值不一樣而已兴革。

軟聚類的mfuzz函數(shù)基于的是e1071包的fuzzy c-means算法。對(duì)于軟過(guò)濾而言，聚類中心點(diǎn)來(lái)源于所有聚類成員的權(quán)重值杂曲。在圖中的membership值可以用mfuzz.plot來(lái)展現(xiàn)庶艾。你也可以用mfuzz.plot2來(lái)看，其會(huì)有更多的選項(xiàng)擎勘。

值得注意的是咱揍，clustering只會(huì)基于表達(dá)矩陣，不會(huì)使用phenoData的任何信息棚饵。還有煤裙，在mfuzz中重復(fù)會(huì)被當(dāng)作是獨(dú)立的信息，所以他們應(yīng)該提前被算好平均值噪漾，或者放進(jìn)不同的ExpressionSet對(duì)象里面硼砰。

> cl <- mfuzz(yeast.s,c=16,m=1.25)
> mfuzz.plot(yeast.s,cl=cl,mfrow=c(4,4),time.labels=seq(0,160,10))

# center代表的應(yīng)該是你選擇的16個(gè)中心點(diǎn)的表達(dá)模式
## 感覺可以用來(lái)畫圖
> head(cl$centers,2)
     cdc28_0   cdc28_10   cdc28_20   cdc28_30    cdc28_40   cdc28_50    cdc28_60    cdc28_70   cdc28_80  cdc28_90 cdc28_100    cdc28_110  cdc28_120 cdc28_130 cdc28_140 cdc28_150 cdc28_160
1  0.1971169 -1.0925729 -1.6203551 -0.7961482 -0.33954720 -0.1567524 -0.05036767  0.08380756  0.5122518 0.3843354 0.4905732  0.437668149  0.4526805 0.3170533 0.2866267 0.2890568 0.6045731
2 -0.7393245 -0.5872038  0.2438611 -0.1883262  0.03321276 -1.0122666 -0.38203192 -0.47328266 -0.6289479 2.1494891 0.5371715 -0.001270464 -0.5672875 0.2288121 0.2331435 0.5896372 0.5646142

# size代表的是各個(gè)聚類的基因數(shù)目 
> head(cl$size,2)
[1] 175 244

# cluster代表的是基因所屬的membership score最高的那個(gè)簇
> head(cl$cluster,5)
YDR132C YMR012W YLR214W YLR116W YDR203W 
      4      11      16      13      16 

# membership代表的是每個(gè)基因?qū)?yīng)16個(gè)簇的membership值
> head(cl$membership,5)
                  1            2            3          4           5           6           7           8          9          10           11          12          13          14          15         16
YDR132C 0.014386854 0.0023657075 0.0478663840 0.47749014 0.007117060 0.070156647 0.010004348 0.051767933 0.01711759 0.028929516 0.0064926466 0.011933880 0.110039692 0.025007449 0.033967103 0.08535705
YMR012W 0.082025023 0.1807366237 0.0088059459 0.00371571 0.012467968 0.002590922 0.127711735 0.007948114 0.06412901 0.004626184 0.3408675360 0.105343965 0.008278282 0.011112051 0.012020122 0.02762081
YLR214W 0.130082565 0.0047176611 0.0041702207 0.06382311 0.001340277 0.004636745 0.013298969 0.043267912 0.21204674 0.016651050 0.0850329333 0.023956515 0.010243391 0.003459240 0.022271137 0.36100153
YLR116W 0.002047923 0.0008467741 0.0412749579 0.01409627 0.002758020 0.034171711 0.001234725 0.002968499 0.00083482 0.003580831 0.0006540104 0.003016454 0.846273635 0.023517377 0.012508494 0.01021550
YDR203W 0.083941355 0.0008482787 0.0008575124 0.02562416 0.001318053 0.004043468 0.001274160 0.032185727 0.12237271 0.002649867 0.0125075149 0.003545481 0.005389429 0.001504605 0.007198611 0.69473907

mfuzz.plot(eset,cl,mfrow=c(1,1),colo,min.mem=0,time.labels,new.window=TRUE)

colo可以設(shè)置顏色，min.mem可以設(shè)置membership的閾值

4.1 Setting of parameters for FCM clustering

對(duì)于fuzzy c-means來(lái)說(shuō)怪与，模糊值m和聚類數(shù)c必須提前設(shè)置好夺刑。對(duì)于m，我們應(yīng)該選擇一個(gè)可以防止隨機(jī)數(shù)據(jù)聚類的值分别。值得注意的是遍愿，fuzzy 聚類可以遵守這樣的準(zhǔn)則，隨機(jī)數(shù)據(jù)并不能被聚類耘斩。這相比于硬聚類（例如k-means）來(lái)說(shuō)沼填，是一個(gè)明顯的優(yōu)點(diǎn)。因?yàn)槠浼词乖陔S機(jī)數(shù)據(jù)中括授，也可以檢測(cè)到cluster坞笙。為了達(dá)到這一點(diǎn)，你可以使用下列選項(xiàng)：

• partcoef函數(shù)荚虚，來(lái)檢測(cè)是否在某一特定的m設(shè)置下薛夜，隨機(jī)數(shù)也會(huì)被聚類
• 或者直接計(jì)算

> m1 <- mestimate(yeast.s)
> m1 # 1.15

設(shè)置一個(gè)合理的聚類值c是很有挑戰(zhàn)性的，尤其是那些short time series版述，很有可能就會(huì)有overlapping clusters梯澜。我們可以設(shè)置一個(gè)最大的c值，大到最后出現(xiàn)了一個(gè)空的empty clusters（看 cselection函數(shù)）

# 不太懂repeat值代表了什么
> cselection(yeast.s,m=1.25,crange=seq(4,32,4),repeats=5,visu=TRUE)
          c:4 c:8 c:12 c:16 c:20 c:24 c:28 c:32
repeats:1   4   8   12   16   19   24   27   31
repeats:2   4   8   12   16   20   23   28   30
repeats:3   4   8   12   16   20   23   28   32
repeats:4   4   8   12   16   20   24   28   31
repeats:5   4   8   12   16   20   23   27   32

在cluster centroid之間最小距離 也可以作為簇有效指數(shù)渴析。在這里晚伙，我們可以檢測(cè)不同的c值之間的。我們可以預(yù)期D.min在達(dá)到最合適值之后俭茧，下降幅度會(huì)變低咆疗。你也可以選擇

4.2 Cluster score

Membership值也可以暗示兩個(gè)向量之間的相關(guān)性。如果兩個(gè)基因?qū)τ谝粋€(gè)特定的cluster都有高的membership score母债，那么他們通常來(lái)說(shuō)表達(dá)模式是相似的午磁。我們對(duì)于高于閾值α的基因，叫做這個(gè)cluster的α-core。

membersip score的設(shè)置通忱焯撸可以作為基因的后驗(yàn)篩選牵署。我們可以用acore函數(shù)漏隐。

tmp <- acore(yeast.s,cl,min.acore = 0.5)
# 生成的似乎是個(gè)列表喧半，里面有16個(gè)。就可以知道每個(gè)簇里面含有的基因ID了青责。

5 Cluster stability

FCM參數(shù)的變化也可以體現(xiàn)出cluster的穩(wěn)健性挺据。我們認(rèn)為那些穩(wěn)健的clusters具有某個(gè)特征，即在m的變化下脖隶，也只會(huì)展現(xiàn)出很小的變化扁耐。

cl2 <- mfuzz(yeast.s,c=16,m=1.35)
mfuzz.plot(yeast.s,cl=cl2,mfrow=c(4,4),time.labels=seq(0,160,10))

6 Global clustering structures

軟聚類有趣的一點(diǎn)就是clusters之間的overlap或者coupling。在cluster k和l之間的coupling coefficient 可以定義為：

N是整個(gè)基因表達(dá)矩陣的數(shù)目产阱。如果coupling值越低婉称，說(shuō)明兩者的表達(dá)模式距離越遠(yuǎn)。如果越高构蹬，說(shuō)明表達(dá)模式越相近王暗。

O <- overlap(cl)
Ptmp <-  overlap.plot(cl,over=O,thres=0.05)

7 Mfuzzgui - the graphical user interface for the Mfuzz pack-age

mfuzz有圖形化界面，不過(guò)我沒去用庄敛。

小結(jié)

最近期末考試復(fù)習(xí)太忙了俗壹。。藻烤。绷雏。有空再加上點(diǎn)注意事項(xiàng)。