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

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

Info

Publication number
CN107609348B
CN107609348B CN201710757362.9A CN201710757362A CN107609348B CN 107609348 B CN107609348 B CN 107609348B CN 201710757362 A CN201710757362 A CN 201710757362A CN 107609348 B CN107609348 B CN 107609348B
Authority
CN
China
Prior art keywords
expression
gene
value
data
differential expression
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.)
Expired - Fee Related
Application number
CN201710757362.9A
Other languages
English (en)
Other versions
CN107609348A (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
Shanghai Sanyu Huaxia Gene Technology 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 Shanghai Sanyu Huaxia Gene Technology Co ltd filed Critical Shanghai Sanyu Huaxia Gene Technology 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

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Complex Calculations (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)对全部样本，使用公式
Figure BDA0001392484730000021
来计算噪声的平均波动水平，其中a1和a0是回归常数，使用广义线性回归模型来估计；μ是基因的平均表达值；再对每一个基因的样本，计算该基因的平均表达值μ和基因表达方差σ，并按照CV2＝σ2/μ2计算该基因的波动水平CV2，然后选取那些
Figure BDA0001392484730000022
的基因；
(2)对于经步骤(1)选取的基因，使用主成分分析PCA方法选取差异表达基因；
(3)选取样本中包含差异表达基因的基因表达数据，构成差异表达基因表达谱。
特别地，步骤(二)中的二值化转换方法为：
对每一个差异表达基因聚类，将此基因在所有样本中的表达值分为两类，用数值“1”代替属于高表达水平那一类的基因表达值，用数值“0”代替属于低表达水平那一类的基因表达值。
特别地，步骤(四)的具体方法为：
基于给定N，先将所有样本聚类成N类，然后执行如下循环过程：
(1)设定K的初始值为2，K为差异表达基因分类数目；
(2)将所有差异表达基因分成K类，此时整个数据划分成N×K个块，数据格式是二值化的差异表达基因表达谱；
(3)创建维度为N×K的情形分析表，表中每一个单元的数值等于对应数据块中二值化数值“1”的个数，使用超几何分布分析方法，对情形分析表的每一个单元计算p值，对于第i行，第j列的p值计算如下
Figure BDA0001392484730000031
x是表中第i行第j列的数值，J是表中第j列的和，I是第i行的和，M是表中数字的总和，
Figure BDA0001392484730000032
使用“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)来计算噪声的平均波动水平
Figure BDA0001392484730000041
Figure BDA0001392484730000042
其中a1和a0是回归常数，使用广义线性回归模型来估计；μ是基因的平均表达值。再对每一个基因的样本，计算该基因的平均表达值μ和基因表达方差σ，并按照CV2＝σ2/μ2计算该基因的波动水平CV2。然后选取那些
Figure BDA0001392484730000043
的基因。
(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值计算如下
Figure BDA0001392484730000051
x是表中第i行第j列的数值，J是表中第j列的和，I是第i行的和，M是表中数字的总和，
Figure BDA0001392484730000052
使用‘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 (4)

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

Families Citing this family (1)

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

Citations (2)

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

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8340914B2 (en) * 2004-11-08 2012-12-25 Gatewood Joe M Methods and systems for compressing and comparing genomic data

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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
Coupled two-way clustering analysis of gene microarray data;Gad Getz等;《Proc. Natl. Acad. Sci., 2000》;20001024;第97卷(第22期);第12079-12084页 *
基于生物知识的生物芯片表达谱数据分析研究;方焯;《中国博士学位论文全文数据库 基础科学辑》;20080315;第1-111页 *

Also Published As

Publication number Publication date
CN107609348A (zh) 2018-01-19

Similar Documents

Publication Publication Date Title
Wirth et al. Mining SOM expression portraits: Feature selection and integrating concepts of molecular function
Raychaudhuri et al. Basic microarray analysis: grouping and feature reduction
Patruno et al. A review of computational strategies for denoising and imputation of single-cell transcriptomic data
CN110797089B (zh) 一种基于单细胞rna测序数据识别细胞类型的方法
CN112750502B (zh) 二维分布结构判定的单细胞转录组测序数据聚类推荐方法
CN115240772B (zh) 一种基于图神经网络的解析单细胞通路活性的方法
CN110377605B (zh) 一种结构化数据的敏感属性识别与分类分级方法
CN114091603A (zh) 一种空间转录组细胞聚类、分析方法
Shu et al. Performance assessment of kernel density clustering for gene expression profile data
CN111292807B (zh) 一种单细胞转录组数据中分析双细胞的方法
CN107609348B (zh) 高通量转录组数据样本分类数目估计方法
CN110544047A (zh) 一种不良数据辨识方法
CN110060735B (zh) 一种基于k-mer组群分割的生物序列聚类方法
CN116564409A (zh) 基于机器学习的转移性乳腺癌转录组测序数据识别方法
US8180775B2 (en) Computer-implemented method for clustering data and computer-readable medium encoded with computer program to execute thereof
Abu-Jamous et al. Identification of genes consistently co-expressed in multiple microarray datasets by a genome-wide Bi-CoPaM approach
CN116976574A (zh) 一种基于两阶段混合聚类算法的建筑负荷曲线降维方法
Alexander et al. Capturing discrete latent structures: choose LDs over PCs
Liu et al. Assessing agreement of clustering methods with gene expression microarray data
CN109033743B (zh) 一种降低单细胞转录组数据中技术噪声的方法
Jeuken et al. Pathway analysis through mutual information
Mohammadi et al. Estimating missing value in microarray data using fuzzy clustering and gene ontology
CN109215741B (zh) 基于双超图正则化的肿瘤基因表达谱数据双聚类方法
Costa et al. A symbolic approach to gene expression time series analysis
Matsuda et al. Distortion-free PCA on sample space for highly variable gene detection from single-cell RNA-seq data

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
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.

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

Granted publication date: 20200623

Termination date: 20210829