CN107609348A - 高通量转录组数据样本分类数目估计方法 - Google Patents

高通量转录组数据样本分类数目估计方法 Download PDF

Info

Publication number
CN107609348A
CN107609348A CN201710757362.9A CN201710757362A CN107609348A CN 107609348 A CN107609348 A CN 107609348A CN 201710757362 A CN201710757362 A CN 201710757362A CN 107609348 A CN107609348 A CN 107609348A
Authority
CN
China
Prior art keywords
gene
expression
binaryzation
value
mtr
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710757362.9A
Other languages
English (en)
Other versions
CN107609348B (zh
Inventor
邹欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Sanyu Huaxia Gene Technology Co ltd
Original Assignee
Top Chia (shanghai) Gene Bioengineering Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Top Chia (shanghai) Gene Bioengineering Co Ltd filed Critical Top Chia (shanghai) Gene Bioengineering Co Ltd
Priority to CN201710757362.9A priority Critical patent/CN107609348B/zh
Publication of CN107609348A publication Critical patent/CN107609348A/zh
Application granted granted Critical
Publication of CN107609348B publication Critical patent/CN107609348B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明涉及一种高通量转录组数据样本分类数目估计方法。首先找出差异表达基因，然后将每一个差异表达基因的表达值进行二值化转换，再对样本和差异表达基因分别进行聚类，并使用富集分析方法建立数据的二值化表达模型。对每一个给定的聚类数目，评估其相应的二值化表达模型，从而判定算法是否收敛。相对于现有技术，本发明对样本分类间相似性进行评估，能保证所得到的每一类样本都有区别于其他样本的独特特征。本发明还有很好的噪声鲁棒性，同时可以很好的降低人工判定带来的偏差。

Description

高通量转录组数据样本分类数目估计方法
技术领域
本发明涉及一种高通量转录组数据样本分类数目估计方法，是针对高通量转录组学数据的分析方法，属于生物数据分析领域。
背景技术
高通量组学分析，在生物医学研究领域有着极广泛的应用。相较于传统的分子生物学研究，组学分析可以大大提高分析通量，即一次可以对上万个目标分子进行测量。也正因为如此，组学数据的数据结构复杂性大大增加了。因此统计分析在组学数据解读中发挥了重要的作用。聚类分析方法在转录组学数据分析中得到了广泛应用。但是目前对于如何将转录组数据进行有生物学意义的分类还没有一个共识，其中关键的问题是分类数目的确定。已有的方法可以分为两大类：
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结束。
0 1 0 1
1 1 1 0
1 1 1 0
步骤(五)分析步骤(四)得到的二值化N×K表达模型矩阵。对于每一个给定N，我们评估所得到的二值化N×K矩阵中每一列的相似性。我们定义，N的值就是使二值化N×K矩阵没有相同列的最大值，即给定当前N，若N类样本中有至少两类的表达谱特征相同，也就是二值化N×K表达模型矩阵中有至少两列的值完全相同，则算法收敛，N-1即是所求的结果；否则N＝N+1，重复步骤(四)和(五)。
一个有重复列的二值化矩阵的例子如下。此时N＝4，但是其中第一列和第三列完全相同，所以N的值设为3，同时算法收敛。
0 1 0 1
1 1 1 0
由上述过程可知，本发明提出的方法可以自动判定算法收敛，这极大的降低了人工干预所造成的误差。并且，本发明提出的方法可以确保不同的样本类别间有显著的统计差异性，并且可以有效区分生物学过程造成的基因表达波动和随机噪声波动。

Claims (5)

1.一种高通量转录组数据样本分类数目估计方法，通过循环过程自动判断收敛，对数据聚类时的分类数目作出估计，其特征在于包括如下步骤：
(一)由高通量转录组数据样本集筛选出差异表达基因，由差异表达基因的表达数据构成差异表达基因表达谱；
(二)在上述差异表达基因表达谱中，对每一个差异表达基因的表达值进行二值化转换；
(三)设定N的初始值为2，N为样本分类数目；
(四)把样本分成N类，同时把差异表达基因分成K类，使用富集分析方法建立数据的二值化N×K表达模型矩阵；
(五)分析步骤(四)得到的二值化N×K表达模型矩阵，若N类样本中有至少两类的表达谱特征相同，也就是二值化N×K表达模型矩阵中有至少两列的值完全相同，则算法收敛，此时输出N的值，N＝N-1；否则，N＝N+1，重复步骤(四)和(五)。
2.根据权利要求1所述的高通量转录组数据样本分类数目估计方法，其特征在于：所述步骤(一)的具体步骤包括：
(1)对全部样本，使用公式来计算噪声的平均波动水平，其中a1和a0是回归常数，使用广义线性回归模型来估计；μ是基因的平均表达值；再对每一个基因的样本，计算该基因的平均表达值μ和基因表达方差σ，并按照CV2＝σ2/μ2计算该基因的波动水平CV2，然后选取那些的基因；
(2)对于经步骤(1)选取的基因，使用主成分分析PCA方法选取差异表达基因；
(3)选取样本中包含差异表达基因的基因表达数据，构成差异表达基因表达谱。
3.根据权利要求1所述的高通量转录组数据样本分类数目估计方法，其特征在于：所述步骤(二)中的二值化转换方法为：对每一个差异表达基因聚类，将此基因在所有样本中的表达值分为两类，用数值“1”代替属于高表达水平那一类的基因表达值，用数值“0”代替属于低表达水平那一类的基因表达值。
4.根据权利要求1所述的高通量转录组数据样本分类数目估计方法，其特征在于：所述步骤(四)的具体方法为：
基于给定N，先将所有样本聚类成N类，然后执行如下循环过程：
(1)设定K的初始值为2，K为差异表达基因分类数目；
(2)将所有差异表达基因分成K类，此时整个数据划分成N×K个块，数据格式是二值化的差异表达基因表达谱；
(3)创建维度为N×K的情形分析表，表中每一个单元的数值等于对应数据块中二值化数值“1”的个数，使用超几何分布分析方法，对情形分析表的每一个单元计算p值，对于第i行，第j列的p值计算如下
<mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>J</mi> </mtd> </mtr> <mtr> <mtd> <mi>x</mi> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>M</mi> <mo>-</mo> <mi>J</mi> </mtd> </mtr> <mtr> <mtd> <mi>I</mi> <mo>-</mo> <mi>x</mi> </mtd> </mtr> </mtable> </mfenced> </mrow> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>M</mi> </mtd> </mtr> <mtr> <mtd> <mi>I</mi> </mtd> </mtr> </mtable> </mfenced> </mfrac> </mrow>
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)。
5.根据权利要求1所述的高通量转录组数据样本分类数目估计方法，其特征在于：对样本和差异表达基因进行分类时，均使用K均值的聚类方法。
CN201710757362.9A 2017-08-29 2017-08-29 高通量转录组数据样本分类数目估计方法 Expired - Fee Related CN107609348B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710757362.9A CN107609348B (zh) 2017-08-29 2017-08-29 高通量转录组数据样本分类数目估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710757362.9A CN107609348B (zh) 2017-08-29 2017-08-29 高通量转录组数据样本分类数目估计方法

Publications (2)

Publication Number Publication Date
CN107609348A true CN107609348A (zh) 2018-01-19
CN107609348B CN107609348B (zh) 2020-06-23

Family

ID=61056290

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710757362.9A Expired - Fee Related CN107609348B (zh) 2017-08-29 2017-08-29 高通量转录组数据样本分类数目估计方法

Country Status (1)

Country Link
CN (1) CN107609348B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111681710A (zh) * 2020-06-03 2020-09-18 中国人民解放军军事科学院军事医学研究院 基于基因表达特征的细胞分类方法、装置和电子设备

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103329138A (zh) * 2011-01-19 2013-09-25 皇家飞利浦电子股份有限公司 用于处理基因组数据的方法
US20150242458A1 (en) * 2004-11-08 2015-08-27 Mitotech, Llc Methods and systems for compressing and comparing genomic data
CN105740651A (zh) * 2016-03-07 2016-07-06 吉林大学 一种特定癌症差异表达基因调控网络的构建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150242458A1 (en) * 2004-11-08 2015-08-27 Mitotech, Llc Methods and systems for compressing and comparing genomic data
CN103329138A (zh) * 2011-01-19 2013-09-25 皇家飞利浦电子股份有限公司 用于处理基因组数据的方法
CN105740651A (zh) * 2016-03-07 2016-07-06 吉林大学 一种特定癌症差异表达基因调控网络的构建方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
GAD GETZ等: "Coupled two-way clustering analysis of gene microarray data", 《PROC. NATL. ACAD. SCI., 2000》 *
方焯: "基于生物知识的生物芯片表达谱数据分析研究", 《中国博士学位论文全文数据库 基础科学辑》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111681710A (zh) * 2020-06-03 2020-09-18 中国人民解放军军事科学院军事医学研究院 基于基因表达特征的细胞分类方法、装置和电子设备

Also Published As

Publication number Publication date
CN107609348B (zh) 2020-06-23

Similar Documents

Publication Publication Date Title
CN108595913B (zh) 鉴别mRNA和lncRNA的有监督学习方法
CN109799269B (zh) 基于动态特征重要度的电子鼻气体传感器阵列优化方法
CN109934269B (zh) 一种电磁信号的开集识别方法和装置
CN109034127B (zh) 一种频谱异常检测方法、装置和电子设备
US20090222390A1 (en) Method, program and apparatus for generating two-class classification/prediction model
CN109635010B (zh) 一种用户特征及特征因子抽取、查询方法和系统
CN114091603A (zh) 一种空间转录组细胞聚类、分析方法
CN107480441B (zh) 一种儿童脓毒性休克预后预测的建模方法及系统
CN106951728B (zh) 一种基于粒子群优化和打分准则的肿瘤关键基因识别方法
CN111209939A (zh) 一种具有智能参数优化模块的svm分类预测方法
CN116504314B (zh) 基于细胞动态分化的基因调控网络构建方法
CN107609348A (zh) 高通量转录组数据样本分类数目估计方法
CN116127398B (zh) 一种基于机理模型与多源数据融合的液压泵故障诊断方法
CN113159220A (zh) 基于随机森林的混凝土侵彻深度经验算法评价方法和装置
CN117033912A (zh) 一种设备故障预测方法、装置、可读存储介质及电子设备
CN112817954A (zh) 一种基于多种方法集成学习的缺失值插补方法
CN116776245A (zh) 一种基于机器学习的三相逆变器设备故障诊断方法
CN111916204A (zh) 一种基于自适应稀疏深度神经网络的脑疾病数据评估方法
CN108872142B (zh) 一种波长选择算法中多参数的选择优化方法
CN110610203A (zh) 基于dwt和极限学习机的电能质量扰动分类方法
CN115983534A (zh) 污水处理过程的状态评价方法以及评价系统
CN111383716B (zh) 基因对的筛选方法、装置、计算机设备和存储介质
CN113326971A (zh) 一种基于PCA和Adaboost的隧道交通事故持续时间预测方法
CN107067034A (zh) 一种快速识别红外光谱数据分类的方法及系统
CN109359694B (zh) 一种基于混合协同表示的分类器的图像分类方法和装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
TA01 Transfer of patent application right

Effective date of registration: 20190110

Address after: Room 801-830, Building 53, Yingdong Village, Chenjiazhen, Chongming District, Shanghai, 202150 (Shanghai Smart Island Data Industry Park)

Applicant after: Shanghai Sanyu Huaxia Gene Technology Co.,Ltd.

Address before: Room 202, Building 518, Shenzhuan Highway, Songjiang High-tech Park, Caohejing Development Zone, Xuhui District, Shanghai, 2003

Applicant before: SHANGZHENGDA (SHANGHAI) GENE BIOENGINEERING CO.,LTD.

TA01 Transfer of patent application right
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200623

Termination date: 20210829

CF01 Termination of patent right due to non-payment of annual fee