最近在进行ATAC和RNA-Seq的联合分析，由于处理的材料是时间相关的，所以time-course也是一个可以分析的点。在一位刘潜老哥的帮助下，我找到了一篇靠谱熊 转录组时间序列数据处理 的文章，里面提到了mfuzz这个包。里面的软聚类的思想非常符合我的预期，然后就决定拿这个包进行我time-course的分析。为了更好的分析，我决定先翻译下这个包，了解下这个包的大致思想。
1 Overview
这部分是我偷懒的随便写的。。。。。。
感觉time-course的方法一般就是聚类。常见的聚类分为三种，分别是层次聚类（Hierarchical Clustering）、硬聚类（hard clustering）、软聚类（soft clustering）。层次聚类好像一般就是像热图那种。看了pheatmap的文档，感觉pheatmap就是层次聚类，当然你可以设置k_means，变成硬聚类。硬聚类常见的就是k-menas。软聚类就是我们这会要用到的这个包的核心思路。
2 Installation requirements
见Bioconductor的安装方法。
3 Data pre-processing
数据集是来源于酵母细胞循环表达数据。6178个基因，横跨160分钟的17个时间点。用的是芯片数据。
> 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
第一步，去除那些有超过25%数据缺失的基因。注意这些数据缺失值应该是设为NA。
yeast.r <- filter.NA(yeast, thres=0.25)
49 genes excluded.
这里就如上面的数据一样，一行即一个基因有16个时间点的数据，如果16个时间点里面有25%，即4个时间点都是NA，则剔除这个基因。
> nrow(yeast)
Features
3000
> nrow(yeast.r)
Features
2951
Fuzzy c-means就像其他聚类算法一样，其并不允许有缺失值的存在。所以我们会对剩下那些缺失值（16个数据点里面就缺了1个2个那种）进行填充。用的是对应基因的平均表达值。
对于RNA-Seq来说，你可以加上一些pseudocount，比如0.01。
yeast.f <- fill.NA(yeast.r,mode="mean")
当然，你也可以用（weighted） k-nearest neighbour method。（mode='knn'/'wknn'）。这些方法相比较而言比上面这种简单的方法要好，但需要耗费更多的算力。
3.2 Filtering
许多已经出版的聚类分析包含过滤的步骤，从而来去除那些表达相对比较低的，或者表达不怎么变化的。通常来说，比较受欢迎的就是样本的标准差作为阈值。
tmp <- filter.std(yeast.f,min.std=0)
然而在基因低表达到高表达的过程中，变化是非常平缓的。所以给定阈值筛选并不一定是可靠的，可能是非常武断。因为现在并没有很多有说服力的筛选手段，所以我们还是避免对基因数据做提前的筛选。这可以避免损失一些有生物学重大意义的基因。
比如1,2,4,10,12,13,15。看起来变化很大，但方差可能并不如你想象中的那么大。
Standardisation
由于聚类是在欧几里德空间中进行的，因此基因的表达值被标准化为平均值为零，标准差为1。该步骤确保了在欧几里得空间中具有相似表达模式的基因是相互接近的。
yeast.s <- standardise(yeast.f)
重要的是，Mfuzz认为输入的表达数据是完全经过前期数据标准化的。standardise 并不能代替标准化步骤。注意差异：标准化是为了让不同的样品间可以比较，而Mufzz中standardisation则是让转录本或者基因间可以比较。
4 Soft clustering of gene expression data
聚类可以用来解释基因表达的调控机制。众所周知的，基因的表达并不是开和关的，而是一个逐渐变化的过程。一个聚类算法应该展现出一个基因有多么的符合dominant cluster pattern。软聚类应该是一个非常好的方法，因为其可以利用membership 衡量一个基因 i跟cluster j的关系。
其实就是说基因A跟每个cluster都有关系，无非是membership score的值不一样而已。
软聚类的mfuzz函数基于的是e1071包的fuzzy c-means算法。对于软过滤而言，聚类中心点来源于所有聚类成员的权重值。在图中的membership值可以用mfuzz.plot来展现。你也可以用mfuzz.plot2来看，其会有更多的选项。
值得注意的是，clustering只会基于表达矩阵，不会使用phenoData的任何信息。还有，在mfuzz中重复会被当作是独立的信息，所以他们应该提前被算好平均值，或者放进不同的ExpressionSet对象里面。
> 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代表的应该是你选择的16个中心点的表达模式
## 感觉可以用来画图
> 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代表的是各个聚类的基因数目
> head(cl$size,2)
[1] 175 244
# cluster代表的是基因所属的membership score最高的那个簇
> head(cl$cluster,5)
YDR132C YMR012W YLR214W YLR116W YDR203W
4 11 16 13 16
# membership代表的是每个基因对应16个簇的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可以设置颜色，min.mem可以设置membership的阈值
4.1 Setting of parameters for FCM clustering
对于fuzzy c-means来说，模糊值m和聚类数c必须提前设置好。对于m，我们应该选择一个可以防止随机数据聚类的值。值得注意的是，fuzzy 聚类可以遵守这样的准则，随机数据并不能被聚类。这相比于硬聚类（例如k-means）来说，是一个明显的优点。因为其即使在随机数据中，也可以检测到cluster。为了达到这一点，你可以使用下列选项：
- partcoef函数，来检测是否在某一特定的m设置下，随机数也会被聚类
- 或者直接计算
> m1 <- mestimate(yeast.s)
> m1 # 1.15
设置一个合理的聚类值c是很有挑战性的，尤其是那些short time series，很有可能就会有overlapping clusters。我们可以设置一个最大的c值，大到最后出现了一个空的empty clusters（看 cselection函数）
# 不太懂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之间最小距离 也可以作为簇有效指数。在这里，我们可以检测不同的c值之间的
。我们可以预期D.min在达到最合适值之后，下降幅度会变低。你也可以选择
4.2 Cluster score
Membership值也可以暗示两个向量之间的相关性。如果两个基因对于一个特定的cluster都有高的membership score，那么他们通常来说表达模式是相似的。我们对于高于阈值α的基因，叫做这个cluster的α-core。
membersip score的设置通常可以作为基因的后验筛选。我们可以用acore函数。
tmp <- acore(yeast.s,cl,min.acore = 0.5)
# 生成的似乎是个列表，里面有16个。就可以知道每个簇里面含有的基因ID了。
5 Cluster stability
FCM参数的变化也可以体现出cluster的稳健性。我们认为那些稳健的clusters具有某个特征，即在m的变化下，也只会展现出很小的变化。
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
软聚类有趣的一点就是clusters之间的overlap或者coupling。在cluster k和l之间的coupling coefficient 可以定义为：
N是整个基因表达矩阵的数目。如果coupling值越低，说明两者的表达模式距离越远。如果越高，说明表达模式越相近。
O <- overlap(cl)
Ptmp <- overlap.plot(cl,over=O,thres=0.05)
7 Mfuzzgui - the graphical user interface for the Mfuzz pack-age
mfuzz有图形化界面，不过我没去用。
小结
最近期末考试复习太忙了。。。。有空再加上点注意事项。
