高通量转录组数据样本分类数目估计方法与流程

本发明涉及一种高通量转录组数据样本分类数目估计方法，是针对高通量转录组学数据的分析方法，属于生物数据分析领域。
背景技术：
：高通量组学分析，在生物医学研究领域有着极广泛的应用。相较于传统的分子生物学研究，组学分析可以大大提高分析通量，即一次可以对上万个目标分子进行测量。也正因为如此，组学数据的数据结构复杂性大大增加了。因此统计分析在组学数据解读中发挥了重要的作用。聚类分析方法在转录组学数据分析中得到了广泛应用。但是目前对于如何将转录组数据进行有生物学意义的分类还没有一个共识，其中关键的问题是分类数目的确定。已有的方法可以分为两大类：1、人工检查法。这一类方法将高维数据投影到低维空间中(通常是2维或3维)，然后通过人工审查判定聚类数目。这一类方法受研究者主观因素影响很大。并且由于受到噪声的干扰，不同类别之间的界限有可能被混淆，从而造成误判。2、统计度量法。通过计算某些与分类数目相关的统计量来预测分类数目。例如，通过计算类内的稀疏度和熵来检测分类数目是否合适，或者将分类数目设置为描述数据特征所使用的特征值的数目，等等。这一类方法为分类提供了客观的参考依据，但是，目前很多的统计度量法都还需要人工来判断统计量随着分类数目增加而出现的拐点，从而估计恰当的分类数目，这就增加了主观误差。更重要的是，现有的统计量都无法保证所得到的分类之间是显著统计差异的，因此可能会造成数据解读困难。技术实现要素：针对现有技术存在的问题，本发明提出一种新的基于统计学理论的高通量转录组数据样本分类数目估计方法。本发明方法使用递归循环计算，每次循环将聚类数目n加1，算法收敛时即可估计n的值。针对目前已有方法不能保证分类结果的统计显著性的问题，本发明提出新的方法对样本分类间相似性进行评估，新方法能保证所得到的每一类样本都有区别于其他样本的独特特征。本发明有很好的噪声鲁棒性，同时可以很好的降低人工判定带来的偏差。本发明的技术方案如下：一种高通量转录组数据样本分类数目估计方法，通过循环过程自动判断收敛，对数据聚类时的分类数目作出估计，其特征在于包括如下步骤：(一)由高通量转录组数据样本集筛选出差异表达基因，由差异表达基因的表达数据构成差异表达基因表达谱；(二)在上述差异表达基因表达谱中，对每一个差异表达基因的表达值进行二值化转换；(三)设定n的初始值为2，n为样本分类数目；(四)把样本分成n类，同时把差异表达基因分成k类，使用富集分析方法建立数据的二值化n×k表达模型矩阵；(五)分析步骤(四)得到的二值化n×k表达模型矩阵，若n类样本中有至少两类的表达谱特征相同，也就是二值化n×k表达模型矩阵中有至少两列的值完全相同，则算法收敛，此时输出n的值，n＝n-1；否则，n＝n+1，重复步骤(四)和(五)。特别地，步骤(一)包括如下步骤：(1)对全部样本，使用公式来计算噪声的平均波动水平，其中a1和a0是回归常数，使用广义线性回归模型来估计；μ是基因的平均表达值；再对每一个基因的样本，计算该基因的平均表达值μ和基因表达方差σ，并按照cv2＝σ2/μ2计算该基因的波动水平cv2，然后选取那些的基因；(2)对于经步骤(1)选取的基因，使用主成分分析pca方法选取差异表达基因；(3)选取样本中包含差异表达基因的基因表达数据，构成差异表达基因表达谱。特别地，步骤(二)中的二值化转换方法为：对每一个差异表达基因聚类，将此基因在所有样本中的表达值分为两类，用数值“1”代替属于高表达水平那一类的基因表达值，用数值“0”代替属于低表达水平那一类的基因表达值。特别地，步骤(四)的具体方法为：基于给定n，先将所有样本聚类成n类，然后执行如下循环过程：(1)设定k的初始值为2，k为差异表达基因分类数目；(2)将所有差异表达基因分成k类，此时整个数据划分成n×k个块，数据格式是二值化的差异表达基因表达谱；(3)创建维度为n×k的情形分析表，表中每一个单元的数值等于对应数据块中二值化数值“1”的个数，使用超几何分布分析方法，对情形分析表的每一个单元计算p值，对于第i行，第j列的p值计算如下x是表中第i行第j列的数值，j是表中第j列的和，i是第i行的和，m是表中数字的总和，使用“1”代表p<0.05的情形分析表单元，“0”代表其他单元，如此，整个数据就被一个n×k的二值化的情形分析表矩阵所描述，称为二值化n×k表达模型矩阵；(4)检测二值化n×k表达模型矩阵是否存在相同的行，如果存在，k＝k-1，循环结束；否则k＝k+1，重复步骤(2)至(4)。特别地，对样本和差异表达基因进行分类时，均使用k均值的聚类方法。附图说明图1是本发明方法流程图；图2是表达谱数据特征二值化描述的过程示意图。具体实施方式下面结合图1，对本发明方法流程进行具体描述。实验时，采用转录组rna-seq数据集，该数据集包含y个样本，每个样本表示一个生物学个体，每个样本的维度是x，每一维度表示一个基因。各步骤的具体实施方法以示例的方式说明如下：步骤(一)的目的是从高通量转录组数据样本中找出差异表达基因。因为数据中通常含有几万个基因的表达数据，而差异表达基因只占其一小部分。我们只关心差异表达基因，而其他基因因为没有任何信息含量，要从剩下的分析中滤除。差异基因的筛选方法已有很多，如阈值法、t检验法等，都可以用于实现本发明的步骤(一)。本实施例采取的具体方法是：(1)对全部样本，使用下列公式(1)来计算噪声的平均波动水平其中a1和a0是回归常数，使用广义线性回归模型来估计；μ是基因的平均表达值。再对每一个基因的样本，计算该基因的平均表达值μ和基因表达方差σ，并按照cv2＝σ2/μ2计算该基因的波动水平cv2。然后选取那些的基因。(2)对于经步骤(1)选取的基因，使用主成分分析pca方法选取差异表达基因，我们可以选取前200个在前4个pca基上系数最大的基因。(3)选取样本中包含差异表达基因的基因表达数据，构成差异表达基因表达谱。图2表达谱数据特征二值化描述的过程示意图。其中，(a)表示原始表达谱数据，每一行代表一个样本，每一列代表一个基因。上方的色带表示样本的分类，左方的色带表示基因的分类。(b)表示二值化的基因表达谱，对应的是步骤(二)。(c)表示二值化n×k表达模型矩阵，对应是步骤(四)。步骤(二)的目的是将基因表达值进行二值化处理，由实数变为“1”或“0”。具体方法为：对每一个差异表达基因使用k-means聚类，将此基因在所有样本中的表达值分为两类。用数值“1”代替属于高表达水平那一类的基因表达值，用数值“0”代替属于低表达水平那一类的基因表达值。如此，数据中连续性的基因表达值被二值化了。比如，一共有5个样本，基因a在这5个样本中的表达值为[0.1,0.12,0.09,1.34,2.53],那么二值化后的基因表达值为[0,0,0,1,1]。步骤(三)是设定循环初值。步骤(四)的目的是得到一个二值化n×k表达模型矩阵，其中n是样本分类数目，k是差异表达基因分类数目。具体方法为：基于给定n，我们先将所有样本使用k-means方法聚类成n类，然后执行如下循环过程：(1)设定k的初始值为2，k为差异表达基因分类数目；(2)使用k-means聚类，将所有的差异表达基因聚类成k个元基因组(metagene)，此时整个数据划分成n×k个块，数据格式是二值化的差异表达基因表达谱；同时，使用k-means方法把所有差异表达基因聚成k类；(3)创建维度为n×k的情形分析表，表中每一个单元的数值等于对应数据块中二值化数值‘1’的个数，使用超几何分布分析方法，对情形分析表的每一个单元计算p值，对于第i行，第j列的p值计算如下x是表中第i行第j列的数值，j是表中第j列的和，i是第i行的和，m是表中数字的总和，使用‘1’代表p<0.05的情形分析表单元，‘0’代表其他单元，如此，整个数据就被一个n×k的二值化的情形分析表矩阵所描述；(4)检测二值化n×k矩阵是否存在相同的行，如果存在，k＝k-1，循环结束；否则k＝k+1，重复步骤(2)至(4)。一个有重复行的二值化矩阵的例子如下，此时k＝3，但是其中第二行和第三行完全相同，所以k的值设为2，同时步骤3结束。010111101110步骤(五)分析步骤(四)得到的二值化n×k表达模型矩阵。对于每一个给定n，我们评估所得到的二值化n×k矩阵中每一列的相似性。我们定义，n的值就是使二值化n×k矩阵没有相同列的最大值，即给定当前n，若n类样本中有至少两类的表达谱特征相同，也就是二值化n×k表达模型矩阵中有至少两列的值完全相同，则算法收敛，n-1即是所求的结果；否则n＝n+1，重复步骤(四)和(五)。一个有重复列的二值化矩阵的例子如下。此时n＝4，但是其中第一列和第三列完全相同，所以n的值设为3，同时算法收敛。01011110由上述过程可知，本发明提出的方法可以自动判定算法收敛，这极大的降低了人工干预所造成的误差。并且，本发明提出的方法可以确保不同的样本类别间有显著的统计差异性，并且可以有效区分生物学过程造成的基因表达波动和随机噪声波动。当前第1页12