聚类分析之k-means
聚类分析，是我们在各种组学分析中经常需要使用到的一种方法。它包括两类聚类方法，分别为层次聚类和非层次聚类。层次聚类包括合并法、分解法和树状图，非层次聚类则包括划分聚类和谱聚类。k-means聚类属于非层次聚类法，已经在多篇文献中被使用，是一种最常见的聚类分析方法。正好，本公众号最近在写一篇转录组的文章，需要使用到这一方法。因而，本文将结合文献解读和代码解读相结合的方式来介绍其中的一种聚类方法k-means。
百度百科对k-means的定义为：k均值聚类算法（k-means clustering algorithm）是一种迭代求解的聚类分析算法，其步骤是，预将数据分为K组，则随机选取K个对象作为初始的聚类中心，然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心。聚类中心以及分配给它们的对象就代表一个聚类。听起来好像很复杂，但其实原理很简单。抽象来看，就是预先设置几个点，再计算其余点与这几个点的距离，从而对这些点进行分析。随后，计算已经分类的这几个点的中心，再重新对这些点进行再分类，直到满足某一条件从而终止。终止的条件可以是所有点不会被分配给不同的聚类或者聚类中心不发生变化。
实现这一方法的途径有很多，本公众号则希望介绍两种比较具有代表性的实现方法：一是直接使用mev软件进行k-means分类；二是通过编写python脚本实现k-means分类。
一、使用Mev软件进行k-means分类
MeV（MultiExperiment Viewer），通用微阵列分析工具，是一款TIGR推出的微阵列分析软件包之一。该软件可以运用各种算法对格式化好的微阵列数据进行聚类、统计、显示、分析，这其中也包括了k-means分类。发表在Journal of Plant GrowthRegulation期刊的“Regulatory Factors of Leaf Senescence areAffected in Arabidopsis Plants Overexpressing the Histone MethyltransferaseSUVH2”、Journal of Experimental Botany期刊的“Systems analysis of cis-regulatory motifs in C4 photosynthesis genesusing maize and rice leaf transcriptomic data during a process of de-etiolation”、Journal of Plant Physiology期刊的“Transcriptomicanalysis reveals the contribution of auxin on the differentially developedcaryopses on primary and secondary branches in rice”和发表在Frontierin plant scenice期刊的“Genome-Wide Identification of lncRNAs Involved in FertilityTransition in the Photo-Thermosensitive Genic Male Sterile Rice Line”等一系列的文章都有使用到mev软件进行k-means分类。这也表明使用Mev进行k-means分析是可靠的。那么我们该如何使用这一软件呢？。接下来，本文将就使用Mev软件进行详细的介绍。
1、打开mev软件，选择方框1中‘File->Load Data’选项。
2、弹出图2的对话框，点击方框2中的Browse选项。
3、弹出图3的文件选择对话框，在文件名即方框3中填入文件的绝对路径。
4、接着图4中会加载相应的数据，其中方框4显示的文件路径，方框5显示的是文件内容。
5、接着图5中会加载相应的数据，其中方框6就是以热图的形式显示你所选择的数据。
6、在图6的方框7中选择“Analysis->Clustering->FOM”选项。
7、在图7的方框8中点击“ok”选项。
8、在图8的方框9中则可以选择弹出FOM图，其中方框10中展示的FOM图，在这里我们选择下降趋势变缓的最优簇数，即9。
9、在图9的方框11中则可以选择“Cluster->k-Means/Medians Cluster”选项，弹出相应对话框。
10、在图10的方框12中的“Number of clusters”中填入上一步确定的cluster数，而在“Maximum iterations”中填入迭代次数，这里选择的是默认值。
11、在图11的方框13中的展示的是各簇的表达情况，而在方框13中的方框14展示的则是所有簇的表达情况。
最终，通过上述步骤即可完成相应的k-means分类。
二、编写python脚本实现k-means分类
Scikit-learn 是开源的 Python 库,通过统一的界面实现机器学习、预处理、交叉验证及可视化算法。这一方法也被发表在Plant Journal期刊的“Transcriptome and epigenome analyses of vernalization in Arabidopsisthaliana”和发表在planta期刊的“TranscriptionFactor Networks in Leaves of Cichorium endivia: New Insights into theRelationship between Photosynthesis and Leaf Development”所使用。因为，本推文也将介绍如何使用python编写脚本，从而实验k-means分类。
对于k-means分类而言，其中最难的一点就是如何确定最优簇数。上文提高的Mev软件使用的是FOM进行最优簇数的确定，而对于自己编写脚本而言，如何得到最优簇数呢。这里主要介绍两种鉴定最优簇数的方法，一种是SSE(Sum of Square Error)即误差平方和法，一种轮廓法。轮廓分析可用于研究生成的聚类之间的分离距离。对于轮廓法而言，轮廓系数更加偏向于接近于1，则表明该簇数最优。而对于SSE计算方法而言，通常来说SSE越小，精确度越高。
接下来，我们将使用sikit-learn包并结合轮廓系数和SSE两者对数据进行k-means分类。
1、加载相应的包
2、定义相应的函数其中，float_f函数用于将数值转换为浮点数，line_strip函数用于将文件每行的换行符去除，write_line_end函数用于添加相应的换行符，savefile2函数用于保存数据，loaddatas16函数用于加载数据，其中加载完数据，会输出三个对象，分别为head,charac,target，其中head是首行，charac则是用于分类的元数据，target则对应相应的基因。
所有的分析都离不开数据，因此，准备数据就是非常重要的一件事了。对于本推文而言，数据中首行（T_1至T_8）为用于分类的指标，首列（gene1_gene3）为分类的数据，具体的格式详见下表：
3、轮廓系数和SSE计算
本步骤得到的轮廓数值如下表，其中首列为假定的分类数，第二类则为假定分类数下计算得到的轮廓系数，第三列则为假定分类数下计算得到的SSE。比如本例中就设定了2-18簇，那么就会依次计算2-18簇的轮廓系数和SSE。
4、根据轮廓系数和SSE数值对数据进行k-means分类。
这一步输出的结果文件为cluster分类的结果，其中，第一列为分类的数据，而第二列则为分类的簇数。部分数据如下表：
如果需要对数据进行可视化，可以参照官方指南对进行https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html。但是这一步非必须，所以本文未进行展示。但是对于本文的结果而言，如果展示分为9簇的轮廓图则为下图：
最后，综上所示，即可完成k-means的分类。
三 惯例小结
本公众号开发的相关软件，Multi-omics Hammer软件和Multi-omics Visual软件欢迎大家使用。文末是本公众号在其他平台的账户，也欢迎大家关注并多提意见。
简书：WJ的生信小院
博客园：生信小院
最后，也欢迎各位大佬能够在本平台上：1传播和讲解自己发表的论文；2：发表对某一科研领域的看法；3：想要达成的合作或者相应的招聘信息；4：展示自己以寻找博后工作或者博士就读的机会;5:博导提供博后工作或者博士攻读机会，都可以后台给笔者留言。希望本平台在进行生信知识分享的同时，能够成为生信分析者的交流平台，能够实现相应的利益互补和双赢（不一定能实现，但是梦想总得是有的吧）。
