CN110910955A - 一种易感基因罕见变异位点纵向分析模型的建立方法 - Google Patents

一种易感基因罕见变异位点纵向分析模型的建立方法 Download PDF

Info

Publication number
CN110910955A
CN110910955A CN201911002493.1A CN201911002493A CN110910955A CN 110910955 A CN110910955 A CN 110910955A CN 201911002493 A CN201911002493 A CN 201911002493A CN 110910955 A CN110910955 A CN 110910955A
Authority
CN
China
Prior art keywords
gene
rare
mutation
function
regression
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
CN201911002493.1A
Other languages
English (en)
Other versions
CN110910955B (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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN201911002493.1A priority Critical patent/CN110910955B/zh
Publication of CN110910955A publication Critical patent/CN110910955A/zh
Application granted granted Critical
Publication of CN110910955B publication Critical patent/CN110910955B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/50—Mutagenesis
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Bioethics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明生物信息数据处理技术领域，提出一种易感基因罕见变异位点纵向分析模型的建立方法，包括以下步骤：获取待分析的病人样本的全基因组序列变异数据；观察统计病人样本中基因上的基因数量突变观察值，对基因进行截断负二项回归，并构建广义线性回归函数；采用最大似然估计函数计算截断负二项回归的系数以及基因罕见变异等位基因数估算值的期望；计算基因的突变观察值与回归估计基线突变数的标准化偏移残差；将所述标准化偏移残差转换为统计显著性程度；根据预设的阈值剔除基因中的显著基因，然后重复上述步骤对截断负二项回归系数进行重新拟合，至病人样本中的所有显著基因被剔除，得到易感基因罕见变异负荷的纵向分析模型。

Description

一种易感基因罕见变异位点纵向分析模型的建立方法
技术领域
本发明涉及生物信息数据处理技术领域，更具体地，涉及一种疾病易感基因罕见变异位点鉴定纵向分析模型的建立方法。
背景技术
罕见基因变异被认为是复杂遗传疾病的主要成因之一，包括二型糖尿病、心力衰竭、骨质疏松症等遗传疾病。例如，基因TREM2和APP中的罕见变异被报道与老年痴呆相关等。但目前用于发掘罕见易感位点的统计方法还不够强大，因此迫切需要更强大有效的新方法发现更多人类复杂疾病中的罕见变异。
现有的复杂疾病罕见变异分析方法主要采用病例对照策略，它继承了基于芯片的全基因组关联分析方法的思想，利用全基因组罕见变异分析。考虑在单个罕见变异位点统计检验低功效问题，研究人员也提出了基于多个变异的关联分析策略，即区域关联性检测。该策略通常联合同一基因或区域内的多个变异位点一起考虑等位基因与疾病相关性，比较病人与正常人突变数的差异评估相关性，或者根据统计模型中方差成分，鉴定易感基因。基因水平多位点关联检测通常比单个位点关联检测更有效。但总体而言这些分析方法的统计效能偏低，难以有效地检测出疾病的罕见易感突变。此外，由于病例、对照分析策略既要测序大量病人样本，又要测序等量正常人样本，成本较高。
发明内容
本发明为克服上述现有技术所述的基因水平检测分析统计效能低的缺陷，提供一种易感基因罕见变异位点纵向分析模型的建立方法。
为解决上述技术问题，本发明的技术方案如下：
一种易感基因罕见变异位点的纵向分析模型建立方法，包括以下步骤：
S1：获取待分析的病人样本的全基因组测序或外显子组序列变异数据；
S2：观察统计病人样本中基因i上全部罕见变异等位基因数量的突变观察值yi，对所述病人样本中基因组的所有基因进行截断负二项回归，构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归模型；
S3：根据所述广义线性回归函数采用最大似然估计函数计算截断负二项回归的系数，以及基因i罕见变异等位基因数估算值的期望；
S4：根据所述基因i罕见变异等位基因数估算值的期望，计算所述病人样本中基因i的突变观察值与回归估计基线突变数的标准化偏移残差；
S5：将所述标准化偏移残差转换为统计显著性程度；
S6：通过预设的阈值剔除所述基因i中的显著基因，然后跳转执行S2步骤，重新拟合得到截断负二项回归系数，至病人样本中的所有显著基因被剔除，得到易感基因罕见变异负荷的纵向分析模型。
本技术方案中，采用截断负二项分布模型拟合罕见变异分布的基准线，用于评估易感基因的相对罕见变异负荷，具体地，对病人样本中的基因组上所有基因做截断负二项回归，然后基于最大似然估计函数计算回归系数，并计算每个基因实际突变数与回归估计基线突变数的偏移残差值，将该偏移残差值在基因组水平纵向地标准化，然后将偏移残差值转换为统计显著性程度p-值，其中，p-值越小表示显著性程度越高，也表示该基因相对于估计的基线有更多的突变量，同时提示该基因与疾病的相关性越强。根据计算得到的统计显著性p-值，通过预设的阈值剔除相应显著基因，然后，用病人样本中余下的基因重复上述步骤，重新拟合出截断负二项回归的系数并计算统计显著性p-值，重复执行上述步骤，直至病人样本中没有可以被剔除的显著基因。由最后迭代所得到的截断负二项回归系数所建立的易感基因罕见变异位点的纵向分析模型最接近零假设，能够减少有高突变数的易感基因对非易感突变基线值估算的影响。
优选地，S1步骤中的病人样本中包括基因i，且基因i中含有mi个罕见变异位点，一个罕见变异位点j中含有ni,j个变异等位基因，其中，i、mi、ni,j分别取正整数。
优选地，S2步骤中，其具体步骤如下：
S21：观察统计基因i上每个变异位点的在病人样本中的加权变异等位基因数ci,j，其计算公式如下：
ci,j＝ni,j*wi,j
其中，ni,j表示罕见变异位点j在病人样本中的实际变异等位基因数量；wi,j表示变异位点j上的功能评分，且wi,j为的正整数，如果加权值缺失则取1；
S22：根据所述加权变异等位基因数ci,j，计算基因i上全部罕见变异等位基因数量的突变观察值yi，其计算公式如下：
Figure BDA0002241764820000031
其中，所述基因i的突变观察值yi服从期望为μi，分布参数为θ的负二项分布；
S23：计算负二项分布的概率质量函数，其计算公式如下：
Figure BDA0002241764820000032
Figure BDA0002241764820000033
其中，Γ(·)表示伽马函数；e是自然对数符号；β为待拟合的回归系数；γ表示分布参数系数；
S24：截取所述基因i中罕见变异等位基因数为0至t的基因，让变异等位基因数服从截断负二项分布模型，其中，在t点截断的概率质量函数如下：
Figure BDA0002241764820000034
其中，g(y|μi,θ,t)表示基因i上的变异等位基因数目为y时的概率，y＝t+1,t+2,...，且t＝0,1,2,...；
S25:构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归函数，所述广义线性回归的连接函数如下式所示：
Figure BDA0002241764820000035
其中，x1,i是基因i的编码区长度，x2,i是基因i的频率评分，x3,i是基因i编码区长度与频率评分的乘积，x4,i是基因i的错义突变保守评分，x5,i是基因i的失去功能突变保守评分，x5,i是基因i的鸟嘌呤和胞嘧啶所占的比率；参数β0,...,β6表示分别的回归系数；EXP(yi)表示基因i的突变观察值yi的期望。
优选地，S3步骤中，所述基于截断负二项分布的最大似然估计函数的计算公式如下：
Figure BDA0002241764820000041
所述最大似然估计函数取对数为
Figure BDA00022417648200000411
li＝lnp(yi|μi,θ)-ln(1-F(X≤t|μi,θ)
其中，
Figure BDA0002241764820000042
优选地，S3步骤中，其具体步骤如下：
S31：根据所述最大似然估计函数分别对期望μi和分布参数θ进行求导，估计截断负二项回归系数参数
Figure BDA0002241764820000043
和参数
Figure BDA0002241764820000044
其中，对所述最大似然估计函数对μi求导的计算公式如下：
Figure BDA0002241764820000045
令
Figure BDA0002241764820000046
得到待拟合的回归系数β的求导计算公式：
Figure BDA0002241764820000047
对所述最大似然估计函数对θ求导的计算公式如下：
Figure BDA0002241764820000048
其中，ψ(·)为Digamma函数，且
Figure BDA0002241764820000049
令θ＝eγ，得到对分布参数系数γ的求导计算公式：
Figure BDA00022417648200000410
根据上述求导计算公式，当导数等于0时，即得到使似然估计函数的对数最大的截断负二项回归系数参数
Figure BDA0002241764820000051
和参数
Figure BDA0002241764820000052
的值；
S32：根据所述参数
Figure BDA0002241764820000053
和参数
Figure BDA0002241764820000054
计算基因i中回归估计基线突变数
Figure BDA0002241764820000055
的期望
Figure BDA0002241764820000056
其计算公式如下：
Figure BDA0002241764820000057
优选地，S4步骤中，计算所述病人样本中基因i的基因突变观察数yi与回归估计基线突变数
Figure BDA0002241764820000058
的标准化偏移残差
Figure BDA0002241764820000059
其计算公式如下：
Figure BDA00022417648200000510
Figure BDA00022417648200000511
Figure BDA00022417648200000512
其中，ei表示基因i的偏移残差，
Figure BDA00022417648200000513
和
Figure BDA00022417648200000514
分别为估算的基因i中罕见变异等位基因数的均值和标准偏差；ri表示原始残差，sign(·)为标注符号函数，ll(·)为分布的自然对数似然函数，
Figure BDA00022417648200000515
表示观察均值；
Figure BDA00022417648200000516
表示估算的基因i中罕见变异等位基因数，均值
Figure BDA00022417648200000517
通过求解所述估算的罕见变异等位基因数方程可得到。
优选地，S5步骤中，采用标准正态分布将所述标准化偏移残差
Figure BDA00022417648200000518
转换为统计显著性程度p-值，其计算公式如下：
Figure BDA00022417648200000519
其中，Φ(·)表示标准正太分布的积累分布函数。
优选地，S6步骤中，所述预设的阈值采用伪发现率FDR，且所述伪发现率FDR≥0.2。
与现有技术相比，本发明技术方案的有益效果是：采用截断负二项分布模型准确拟合罕见变异分布的基准线，并通过基准线有效检测复杂疾病的易感基因，能够有效减少高突变数的易感基因对非易感突变基线值估算的影响，从而有效提高基因水平检测分析统计效能。
附图说明
图1为本实施例的易感基因罕见变异位点的纵向分析模型建立方法流程图。
图2为本实施例中病人样本量为75、100、150、200时的p-值对比图。
图3为本实施例中与对照方法的统计效能比较图。
具体实施方式
附图仅用于示例性说明，不能理解为对本专利的限制；
为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；
对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
下面结合附图和实施例对本发明的技术方案做进一步的说明。
如图1所示，为本实施例的易感基因罕见变异位点的纵向分析模型建立方法的流程图。
本实施例提出一种易感基因罕见变异位点的纵向分析模型建立方法，包括以下步骤：
S1：获取待分析的病人样本的全基因组数据。
本实施例中，病人样本中包括基因i，且基因i中含有mi个罕见变异位点，一个罕见变异位点j中含有ni,j个变异等位基因。
本实施例中，对每个基因i设置有功能评分wi,j＝1,2,3,...。
S2：观察统计病人样本中基因i上全部罕见变异等位基因数量的突变观察值yi，对所述病人样本中基因组的所有基因进行截断负二项回归，并构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归函数。其具体步骤如下：
S21：观察统计基因i上每个变异位点的在病人样本中的加权变异等位基因数ci,j，其计算公式如下：
ci,j＝ni,j*wi,j
其中，ni,j表示罕见变异位点j在病人样本中的实际变异等位基因数量；wi,j表示变异位点j上的功能评分，且wi,j为的正整数，如果加权值缺失则取1；
S22：根据所述加权变异等位基因数ci,j，计算基因i上全部罕见变异等位基因数量的突变观察值yi，其计算公式如下：
Figure BDA0002241764820000061
其中，所述基因i的突变观察值yi服从期望为μi，分布参数为θ的负二项分布；
S23：计算负二项分布的概率质量函数，其计算公式如下：
Figure BDA0002241764820000071
Figure BDA0002241764820000072
其中，Γ(·)表示伽马函数；e是自然对数符号；β为待拟合的回归系数；γ表示分布参数系数；
S24：截取所述基因i中罕见变异等位基因数为0至t的基因，让变异等位基因数服从截断负二项分布模型，其中，在t点截断的概率质量函数如下：
Figure BDA0002241764820000073
其中，g(y|μi,θ,t)表示基因i上的变异等位基因数目为y时的概率，y＝t+1,t+2,...，且t＝0,1,2,...；
S25:构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归函数，所述广义线性回归的连接函数如下式所示：
Figure BDA0002241764820000074
其中，x1,i是基因i的编码区长度，x2,i是基因i的频率评分，x3,i是基因i编码区长度与频率评分的乘积，x4,i是基因i的错义突变保守评分，x5,i是基因i的失去功能突变保守评分，x5,i是基因i的鸟嘌呤和胞嘧啶所占的比率；参数β0,...,β6表示分别的回归系数；EXP(yi)表示基因i的突变观察值yi的期望。
本实施例中，由于大部分基因上不存在或存在极少的罕见突变，低突变的基因数的比例将会膨胀，这样往往会扭曲假定的负二项分布，因此在本实施例中截取所述全基因组数据中罕见变异等位基因数为0至t的基因，使变异等位基因数服从截断负二项分布模型。
S3：根据所述广义线性回归函数采用最大似然估计函数计算截断负二项回归系数，以及基因i罕见变异等位基因数估算值的期望。
本步骤中，基于截断负二项分布的最大似然估计函数的计算公式如下：
Figure BDA0002241764820000081
将最大似然估计函数取对数，得到以下公式：
Figure BDA0002241764820000082
li＝lnp(yi|μi,θ)-ln(1-F(X≤t|μi,θ)
其中，
Figure BDA0002241764820000083
本步骤中，包括以下具体步骤：
S31：根据所述最大似然估计函数分别对期望μi和分布参数θ进行求导，计算截断负二项回归系数参数
Figure BDA0002241764820000084
和参数
Figure BDA0002241764820000085
其中，对所述最大似然估计函数对μi求导的计算公式如下：
Figure BDA0002241764820000086
令
Figure BDA0002241764820000087
得到待拟合的回归系数β的求导计算公式：
Figure BDA0002241764820000088
对所述最大似然估计函数对θ求导的计算公式如下：
Figure BDA0002241764820000089
其中，ψ(·)为Digamma函数，且
Figure BDA00022417648200000810
令θ＝eγ，得到对γ的求导计算公式：
Figure BDA00022417648200000811
Figure BDA0002241764820000091
根据上述求导计算公式，当导数等于0时，即得到使似然估计函数的对数最大的截断负二项回归系数参数
Figure BDA0002241764820000092
和参数
Figure BDA0002241764820000093
的值；
S32：根据所述参数
Figure BDA0002241764820000094
和参数
Figure BDA0002241764820000095
计算基因i中回归估计基线突变数
Figure BDA0002241764820000096
的期望
Figure BDA0002241764820000097
其计算公式如下：
Figure BDA0002241764820000098
S4：根据所述基因i罕见变异等位基因数估算值的期望，计算所述病人样本中基因i的突变观察值与回归估计基线突变数的标准化偏移残差。
本步骤中，计算所述病人样本中基因i的基因突变观察数yi与回归估计基线突变数
Figure BDA0002241764820000099
的标准化偏移残差
Figure BDA00022417648200000910
的计算公式如下：
Figure BDA00022417648200000911
Figure BDA00022417648200000912
Figure BDA00022417648200000913
其中，ei表示基因i的偏移残差，
Figure BDA00022417648200000914
和
Figure BDA00022417648200000915
分别为估算的基因i中罕见变异等位基因数的均值和标准偏差；ri表示原始残差，sign(·)为标注符号函数，ll(·)为分布的自然对数似然函数，
Figure BDA00022417648200000916
表示观察均值；
Figure BDA00022417648200000917
表示估算的基因i中罕见变异等位基因数，均值
Figure BDA00022417648200000918
通过求解所述估算的罕见变异等位基因数方程可得到。
S5：将所述标准化偏移残差转换为统计显著性程度。
本步骤中，采用标准正态分布将所述标准化偏移残差
Figure BDA00022417648200000921
转换为统计显著性程度p-值，其计算公式如下：
Figure BDA00022417648200000919
其中，Φ(·)表示标准正太分布的积累分布函数。
本实施例中，统计显著性程度p-值表示基因i在病人群中相对基线
Figure BDA00022417648200000920
突变数的多寡，p-值越小则表示相对突变数越多，同时也表示更有可能是复杂疾病易感基因。本实施例中的p-值为由偏移残差值ei标准化后计算得到，而偏移残差值ei同时是基因i的罕见变异等位基因数的观察值yi与估算值
Figure BDA0002241764820000101
及其均值
Figure BDA0002241764820000102
计算得到。
S6：通过预设的阈值剔除所述基因i中的显著基因，然后跳转执行S2步骤，重新拟合得到截断负二项回归系数，至病人样本中的所有显著基因被剔除，得到易感基因罕见变异负荷的纵向分析模型。
在本实施例中，采用伪发现率FDR剔除相应的显著基因，且FDR设置为0.2。
本实施例中提出的易感基因罕见变异位点纵向分析模型可以减少有高突变数的易感基因对非易感突变基线值估算的影响，因此用稳定后的模型计算全部基因的偏移残差和p值，能够实现通过纵向比较基因之间的校正后的相对突变数的差异来鉴定复杂疾病易感基因，突破了现有的分析方法受限于对照样本的要求。
本实施例中，采用随机样本的方法检验本实施例提出的方法与使用原有的位点变异数没有加权的方法鉴定易感基因的1型错误，并采用现有的CMC、SKAT、Price、KBAC等四种基于基因关联分析的罕见变异样本对照常用方法与本实施例提出的一种易感基因罕见变异位点的纵向分析模型建立方法进行对比。如图2所示，分别为当病人样本量为75、100、150、200时的p值对比图。其中，图2(a)为表示病人样本量为75时的p值对比图，图2(b)表示病人样本量为100时的p值对比图，图2(c)表示病人样本量为150时的p值对比图，图2(d)表示病人样本量为200时的p值对比图；RUNER表示本实施例提出的易感基因罕见变异位点纵向分析模型，UW-RUNER表示本实施例提出易感基因罕见变异位点纵向分析模型使用原有的位点变异数且没有加权的方法。
由图可知，RUNER的p值最接近均匀分布，即使在样本量为75时的分析效果也相对较为理想，而其他四种对照方法在小样本时出现了很明显的p值膨胀，表示现有的统计关联检测对中小样本量的罕见变异效果不佳，且耗费时间长，存在基因水平检测分析统计效能低的缺陷。
如图3所示，为本实施例与对照方法的统计效能比较图，其中，图3(a)为探测基因TIE1的效能对比图，图3(b)为探测基因TCF4的效能对比图，图3(c)为表示模拟生成100个数据集，每个数据集做一次检测，每次产生的假阳性基因数求均值。在对照实验中，在健康样本中随机地插入易感等位基因产生模拟病人样本。在模拟中，假设多个罕见错义突变发生在基因TCF4和TIE1中，其中，TCF4是精神分裂症的易感基因，TIE1对产生血管病很重要。在该对照实验中，Bonferroni校正控制多重比较谬误为0.05。由图可知，当样本量为75时，RUNER的效能能够达到58％和89％，且当样本100时，RUNER的效能达到75％和93％。相比于其他四种现有的分析方法，该四种方法在小样本时效能远不如RUNER和不加权的UW-RUNER。因此，本实施例提出的易感基因罕见变异位点纵向分析模型显然能够有效克服基因水平检测分析统计效能低的缺陷。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (8)

1.一种易感基因罕见变异位点的纵向分析模型建立方法，其特征在于，包括以下步骤：
S1：获取待分析的病人样本的全基因组测序或外显子组序列变异数据；
S2：观察统计病人样本中基因i上全部罕见变异等位基因数量的突变观察值yi，对所述病人样本中基因组的所有基因进行截断负二项回归，并构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归函数；
S3：根据所述广义线性回归函数采用最大似然估计函数计算截断负二项回归系数，以及基因i罕见变异等位基因数估算值的期望；
S4：根据所述基因i罕见变异等位基因数估算值的期望，计算所述病人样本中基因i的突变观察值与回归估计基线突变数的标准化偏移残差；
S5：将所述标准化偏移残差转换为统计显著性程度；
S6：通过预设的阈值剔除所述基因i中的显著基因，然后跳转执行S2步骤，重新拟合得到截断负二项回归系数，至病人样本中的所有显著基因被剔除，得到易感基因罕见变异负荷的纵向分析模型。
2.根据权利要求1所述的纵向分析模型建立方法，其特征在于：所述S1步骤中的病人样本中包括基因i，且基因i中含有mi个罕见变异位点，一个罕见变异位点j中含有ni,j个变异等位基因，其中i为正整数。
3.根据权利要求2所述的纵向分析模型建立方法，其特征在于：所述S2步骤中，其具体步骤如下：
S21：观察统计基因i上每个变异位点的在病人样本中的加权变异等位基因数ci,j，其计算公式如下：
ci,j＝ni,j*wi,j
其中，ni,j表示罕见变异位点j在病人样本中的实际变异等位基因数量；wi,j表示变异位点j上的功能评分，且wi,j为的正整数，如果加权值缺失则取1；
S22：根据所述加权变异等位基因数ci,j，计算基因i上全部罕见变异等位基因数量的突变观察值yi，其计算公式如下：
Figure FDA0002241764810000011
其中，所述基因i的突变观察值yi服从期望为μi，分布参数为θ的负二项分布；
S23：计算负二项分布的概率质量函数，其计算公式如下：
Figure FDA0002241764810000021
Figure FDA0002241764810000022
θ＝eγ
其中，Γ(·)表示伽马函数；e是自然对数符号；β为待拟合的回归系数；γ表示分布参数系数；
S24：截取所述基因i中罕见变异等位基因数为0至t的基因，让变异等位基因数服从截断负二项分布模型，其中，在t点截断的概率质量函数如下：
Figure FDA0002241764810000023
其中，g(y|μi,θ,t)表示基因i上的变异等位基因数目为y时的概率，y＝t+1,t+2,...，且t＝0,1,2,...；
S25:构建用于预测基因i罕见突变等位基因数的基于截断负二项分布的广义线性回归函数，所述广义线性回归的连接函数如下式所示：
Figure FDA0002241764810000024
其中，x1,i是基因i的编码区长度，x2,i是基因i的频率评分，x3,i是基因i编码区长度与频率评分的乘积，x4,i是基因i的错义突变保守评分，x5,i是基因i的失去功能突变保守评分，x5,i是基因i的鸟嘌呤和胞嘧啶所占的比率；参数β0,...,β6表示分别的回归系数；EXP(yi)表示基因i的突变观察值yi的期望。
4.根据权利要求3所述的纵向分析模型建立方法，其特征在于：所述S3步骤中，所述基于截断负二项分布的最大似然估计函数的计算公式如下：
Figure FDA0002241764810000031
所述最大似然估计函数取对数为
ln L(θ,t)＝∑ili
li＝lnp(yi|μi,θ)-ln(1-F(X≤t|μi,θ)
其中，
Figure FDA0002241764810000032
5.根据权利要求4所述的纵向分析模型建立方法，其特征在于：所述S3步骤中，其具体步骤如下：
S31：根据所述最大似然估计函数分别对期望μi和分布参数θ进行求导，估计截断负二项回归系数参数
Figure FDA0002241764810000033
和参数
Figure FDA0002241764810000034
其中，对所述最大似然估计函数对μi求导的计算公式如下：
Figure FDA0002241764810000035
令
Figure FDA0002241764810000036
得到待拟合的回归系数β的求导计算公式：
Figure FDA0002241764810000037
对所述最大似然估计函数对θ求导的计算公式如下：
Figure FDA0002241764810000038
其中，ψ(·)为Digamma函数，且
Figure FDA0002241764810000039
令θ＝eγ，得到对γ的求导计算公式：
Figure FDA00022417648100000310
根据上述求导计算公式，当导数等于0时，即得到使似然估计函数的对数最大的截断负二项回归系数参数
Figure FDA0002241764810000041
和参数
Figure FDA0002241764810000042
的值；
S32：根据所述参数
Figure FDA0002241764810000043
和参数
Figure FDA0002241764810000044
计算基因i中回归估计基线突变数
Figure FDA0002241764810000045
的期望
Figure FDA0002241764810000046
其计算公式如下：
Figure FDA0002241764810000047
6.根据权利要求5所述的纵向分析模型建立方法，其特征在于：所述S4步骤中，计算所述病人样本中基因i的基因突变观察数yi与回归估计基线突变数
Figure FDA0002241764810000048
的标准化偏移残差
Figure FDA0002241764810000049
的计算公式如下：
Figure FDA00022417648100000410
Figure FDA00022417648100000411
Figure FDA00022417648100000412
其中，ei表示基因i的偏移残差，
Figure FDA00022417648100000413
和
Figure FDA00022417648100000414
分别为估算的基因i中罕见变异等位基因数的均值和标准偏差；ri表示原始残差，sign(·)为标注符号函数，ll(·)为分布的自然对数似然函数，
Figure FDA00022417648100000415
表示观察均值；
Figure FDA00022417648100000416
表示估算的基因i中罕见变异等位基因数，均值
Figure FDA00022417648100000417
通过求解所述估算的罕见变异等位基因数方程可得到。
7.根据权利要求6所述的纵向分析模型建立方法，其特征在于：所述S5步骤中，采用标准正态分布将所述标准化偏移残差
Figure FDA00022417648100000418
转换为统计显著性程度p-值，其计算公式如下：
Figure FDA00022417648100000419
其中，Φ(·)表示标准正太分布的积累分布函数。
8.根据权利要求1～7任一项所述的纵向分析模型建立方法，其特征在于：所述S6步骤中，所述预设的阈值采用伪发现率FDR，所述伪发现率FDR≥0.2。
CN201911002493.1A 2019-10-21 2019-10-21 一种易感基因罕见变异位点纵向分析模型的建立方法 Active CN110910955B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911002493.1A CN110910955B (zh) 2019-10-21 2019-10-21 一种易感基因罕见变异位点纵向分析模型的建立方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911002493.1A CN110910955B (zh) 2019-10-21 2019-10-21 一种易感基因罕见变异位点纵向分析模型的建立方法

Publications (2)

Publication Number Publication Date
CN110910955A true CN110910955A (zh) 2020-03-24
CN110910955B CN110910955B (zh) 2024-03-01

Family

ID=69816199

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911002493.1A Active CN110910955B (zh) 2019-10-21 2019-10-21 一种易感基因罕见变异位点纵向分析模型的建立方法

Country Status (1)

Country Link
CN (1) CN110910955B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112599189A (zh) * 2020-12-29 2021-04-02 北京优迅医学检验实验室有限公司 一种全基因组测序的数据质量评估方法及其应用

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170249434A1 (en) * 2016-02-26 2017-08-31 Daniela Brunner Multi-format, multi-domain and multi-algorithm metalearner system and method for monitoring human health, and deriving health status and trajectory
CN107391962A (zh) * 2017-09-05 2017-11-24 武汉古奥基因科技有限公司 基于多组学分析基因或位点对疾病调控关系的方法
CN109346127A (zh) * 2018-08-09 2019-02-15 中山大学 一种用于检测潜在癌症驱动基因的统计分析方法
US20190228838A1 (en) * 2016-09-26 2019-07-25 Mcmaster University Tuning of Associations For Predictive Gene Scoring

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170249434A1 (en) * 2016-02-26 2017-08-31 Daniela Brunner Multi-format, multi-domain and multi-algorithm metalearner system and method for monitoring human health, and deriving health status and trajectory
US20190228838A1 (en) * 2016-09-26 2019-07-25 Mcmaster University Tuning of Associations For Predictive Gene Scoring
CN107391962A (zh) * 2017-09-05 2017-11-24 武汉古奥基因科技有限公司 基于多组学分析基因或位点对疾病调控关系的方法
CN109346127A (zh) * 2018-08-09 2019-02-15 中山大学 一种用于检测潜在癌症驱动基因的统计分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
MIAO-XIN LI 等: "GATES: A Rapid and Powerful Gene-Based Association Test Using Extended Simes Procedure", 《美国人类遗传学杂志》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112599189A (zh) * 2020-12-29 2021-04-02 北京优迅医学检验实验室有限公司 一种全基因组测序的数据质量评估方法及其应用

Also Published As

Publication number Publication date
CN110910955B (zh) 2024-03-01

Similar Documents

Publication Publication Date Title
Pouyet et al. Background selection and biased gene conversion affect more than 95% of the human genome and bias demographic inferences
US9639657B2 (en) Methods for allele calling and ploidy calling
Capra et al. A model-based analysis of GC-biased gene conversion in the human and chimpanzee genomes
US9916416B2 (en) System and method for genotyping using informed error profiles
US20170206311A1 (en) Method of characterizing sequences from genetic material samples
CN109887546B (zh) 基于二代测序的单基因或多基因拷贝数检测系统及方法
Weber et al. Species delimitation in the presence of strong incomplete lineage sorting and hybridization: Lessons from Ophioderma (Ophiuroidea: Echinodermata)
US20190338349A1 (en) Methods and systems for high fidelity sequencing
CN114724626A (zh) 母体血浆的无创性产前分子染色体核型分析
WO2020063052A1 (zh) 胎儿游离dna浓度获取方法、获取装置、存储介质及电子装置
Yao et al. Limitations of principal components in quantitative genetic association models for human studies
CN110910955A (zh) 一种易感基因罕见变异位点纵向分析模型的建立方法
Sethuraman Estimating genetic relatedness in admixed populations
Zhou et al. Eigenvalue significance testing for genetic association
Biswas et al. Biological averaging in RNA-seq
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
WO2016112539A1 (zh) 确定胎儿核酸含量的方法和装置
Bérard et al. Unsupervised classification for tiling arrays: ChIP-chip and transcriptome
Warde-Farley et al. Mixture model for sub-phenotyping in GWAS
Gymrek et al. A framework to interpret short tandem repeat variations in humans
Mao et al. abSNP: RNA-seq SNP calling in repetitive regions via abundance estimation
Talenti et al. The evolution and convergence of mutation spectra across mammals
Hirt et al. A bayesian framework for genome-wide inference of DNA methylation levels
Bercovich Szulmajster et al. Measuring linkage disequilibrium and improvement of pruning and clumping in structured populations
KR102441856B1 (ko) 중요도 샘플링을 활용한 다중변이 연관연구 방법

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
GR01 Patent grant
GR01 Patent grant