CN113130004A - 一种识别阿尔茨海默症相关生物标志物的相关性分析方法 - Google Patents

Abstract

一种识别阿尔茨海默症相关生物标志物的相关性分析方法，基于联合连通性的稀疏非负矩阵分解，将结构磁共振成像SMRI、单核苷酸多态性位点SNP和基因表达数据同时投影到一个公共特征空间中，还添加了大脑每个区域的连通性信息和遗传数据作为先验知识，以识别与阿尔茨海默病患者、早期认知障碍患者相关的感兴趣区域ROI，风险SNP位点和风险基因。本发明具有更强的相关性分析能力，以及更好的抗噪性能和生物学解释性。

Description

一种识别阿尔茨海默症相关生物标志物的相关性分析方法

技术领域

本发明涉及基于阿尔兹海默症研究的影像遗传学，尤其涉及一种识别阿尔茨海默症相关生物标志物的相关性分析方法。

背景技术

影像遗传学已被广泛用于神经退行性疾病，它可以探索基因对大脑结构和功能的影响，并可以使用大脑成像来评估基因对个体的影响。最近影像遗传学在研究阿尔茨海默症的发病机理和挖掘阿尔茨海默症相关生物标志物方面取得了重大进展。影像遗传学可以揭示微观遗传学和宏观成像之间的联系，检测疾病的生物标记，但是缺乏有效和准确的相关性分析模型。

发明内容

本发明的目的在于提供一种识别阿尔茨海默症相关生物标志物的相关性分析方法，具有更强的相关性分析能力，以及更好的抗噪性能和生物学解释性。

为了达到上述目的，本发明提供一种识别阿尔茨海默症相关生物标志物的相关性分析方法，基于联合连通性的稀疏非负矩阵分解，将结构磁共振成像SMRI、单核苷酸多态性位点SNP和基因表达数据同时投影到一个公共特征空间中，还添加了大脑每个区域的连通性信息和遗传数据作为先验知识，以识别与阿尔茨海默病患者、早期认知障碍患者相关的感兴趣区域ROI，风险SNP位点和风险基因。

所述相关性分析方法的和目标函数如下：

令

和

(H_I)_ij≥0的拉格朗日乘数，

L关于W和H_I的偏导数是：

基于Karush-Kuhn-Tucher条件，Ψ_ijW_ij为0，且

为0，得到W_ij和(H_I)_ij：

其中，X₁代表输入的sMRI影像数据，X₂代表SNP数据，X₃代表基因表达数据；W是三种数据的基矩阵，H_i(i＝1，2，3...)分别是三种数据的系数矩阵，A_i(i＝1，2，3...)代表不同数据之间的邻接矩阵，B_i(i＝1，2，3...)代表每种数据的拉普拉斯矩阵；λ_i、β_i、γ_i(i＝1，2，3…)是控制各约束项稀疏性的参数；i，j分别代表矩阵的第i行和第j列；k代表共表达模块数。

本发明将结构磁共振成像(sMRI)，单核苷酸多态性位点(SNP)和基因表达数据同时投影到一个公共特征空间中，其中，在同一投影方向上具有较大系数的异类变量形成一个公共模块，此外还添加了大脑每个区域的连通性信息和遗传数据作为先验知识，以识别与阿尔茨海默病(AD)患者、早期认知障碍患者(MCI)相关的感兴趣区域(ROI)，风险SNP位点和风险基因，具有更强的相关性分析能力，以及更好的抗噪性能和生物学解释性。

附图说明

图1是本发明提供的一种识别阿尔茨海默症相关生物标志物的相关性分析方法的应用流程图。图2是前七个丰富的生物过程以及这些过程涉及的模块数量。

图3A是三个模块之间的ROI重叠。

图3B是三个模块之间的SNP重叠。

图3C是基因在三个模块之间的重叠。

图4是对模块1中选择的基因进行GO富集分析的结果，水平轴是途径中的基因数量，垂直轴是途径列表，-P值的范围为0到10。

注图中坐标：I-kappaB激酶/NF-kappaB信号转导的正调控9(Positiveregulation of I-kappaB kinase/NF-kappaB signaling9)；细胞因子的产生(Cytokineproduction)；蛋白质定位到染色体，端粒区域(Protein localization to chromosome,telomeric region)；睾丸支持细胞发育(Sertoli cell development)；通过剪接体剪接mRNA(mRNA splicing，via spliceosome)；mRNA加工(mRNA processing)；端粒维持的积极调控(Positive regulation of telomere maintenance)；核质(Nucleoplasm)；核(Nucleus)；细胞质(Cytoplasm)；核基质(Nuclear matrix)；核异染色质(Nuclearheterochromatin)；膜(Membrane)；细胞内膜结合细胞器(Intracellular membrane-bounded organelle)；核糖体(Ribosome)；核染色体，端粒区域(Nuclear chromosome,telomeric region)；聚(A)RNA结合(Poly(A)RNA binding)；蛋白质结合(Proteinbinding)；甲基化组蛋白结合(Methylated histone binding)；泛素蛋白转移酶活性(Ubiquitin-protein transferase activity)；端粒DNA结合(Telomeric DNA binding)；锌离子结合(Zinc ion binding)；组蛋白脱乙酰基酶结合(Histone deacetylasebinding)；RNA结合(RNA binding)。

图5是用模块1中选择的基因构建的PPI网络。

图6是三个重要基因的ROC曲线。

图7A和B分别是与模块1中所选标记关联的SNP-ROI和Gene-ROI的热图。

图8比较了三种算法在不同噪声水平下的重建误差和目标函数值，A和B是在小型数据集上的结果，C和D是在大型数据集上的结果。

具体实施方式

以下根据图1～图8，具体说明本发明的较佳实施例。

在本专利中，提出了一种基于联合连通性的稀疏非负矩阵分解(JCB-SNMF)的新方法，该方法将结构磁共振成像(SMRI)，单核苷酸多态性位点(SNP)和基因表达数据同时投影到一个公共特征空间中，其中，在同一投影方向上具有较大系数的异类变量形成一个公共模块，此外还添加了大脑每个区域的连通性信息和遗传数据作为先验知识，以识别与阿尔茨海默病(AD)患者、早期认知障碍患者(MCI)相关的感兴趣区域(ROI)，风险SNP位点和风险基因。

NMF是功能强大的降维计算框架，可以集成不同的组学数据。NMF可以将非负矩阵分解为两个矩阵，分别是基本矩阵W和系数矩阵H。许多专家学者在其基础上进行了创新。联合非负矩阵算法(JNMF)具有较好的稳定性，但是由于其没有加入先验知识，生物可解释性和抗噪性能较差。联合稀疏网络正则化约束的非负矩阵分解(JSNMNMF)将不同类型数据的邻接矩阵加入到算法中，提高了算法的生物学可解释性。多约束联合非负矩阵分解(MCJNMF)在JNMF的基础上加入了对系数矩阵H的正交约束，使得算法能够得到冗余较小的特征，但目前只适用于两类数据。为了将结构MRI，SNP和基因表达数据整合到ADNI数据库(http://ADNI.loni.usc.edu/)中，本专利基于目前一些算法抗噪性能差、生物学可解释性差的问题，我们提出了一种改进的联合负矩阵分解方法，称为基于联合连通性的稀疏非负矩阵分解(JCB-SNMF)。考虑到MRI，SNPs和基因的特征之间的相关性以及遗传和成像数据中的连通性信息，我们使用邻接矩阵和GraphNet正则化器作为网络正则化约束来提高算法的准确性和抗噪性。GraphNet正则化器是弹性网正则化的修改版本，可以有效地集成诸如连接性之类的生理约束。在本专利中，我们使用所提出算法对SMRI，SNPs和基因表达数据进行了详尽分析，通过对所提出算法选择的重要共表达模块的生物学分析，包括基因本体(GO)富集分析和对所选基因的PPI网络构建，我们发现所选基因参与的生物过程、所选出的ROI与神经退行性疾病显著相关。我们对相互作用最强的基因进行了接收者操作特征(ROC)分析，发现了诊断AD和MCI的潜在基因，包括剪接因子3b亚基1(SF3B1)，核糖体蛋白S20(RPS20)，RNA结合基序蛋白14(RBM14)。此外，重要模块中SNP的大多数基因也与AD密切相关。我们绘制了选定的SNP-ROI对和基因-ROI对的热图，并发现了显著(p<0.05)的关系对，这可能是AD和MCI的生物标记。

1、联合非负矩阵分解(JNMF)

NMF是传统的降维方法，其通用模型如下：

其中，X∈R^n×p代表原始特征矩阵，可以通过NMF将其分解为W和H，W∈R^n×k被称为基础矩阵，H∈R^k×p称为系数矩阵，n是样本数，p是样本特征。

NMF可以很好地减少单个数据的维数，但是它不能同时在多模式数据上执行，因此提出了JNMF来解决该问题，其模型如下：

其中，X_i∈R^n×pi(i＝1，2，..)表示不同模式的原始矩阵，W_i∈R^n×k表示通过分解获得的基本矩阵，H_I∈R^k×pi表示通过分解获得的每个原始矩阵的系数矩阵。

2、联合稀疏网络正则化非负矩阵分解(JSNMNMF)

考虑到由JNMF获得的系数矩阵具有很强的独立性，并且在实际应用中n通常小于p，提出了JSNMNMF算法。JSNMNMF在目标函数中增加了先验知识，以改善结果的生物学相关性，并且还可以通过减少较大的搜索空间来提高模块的效率。在本专利中，为了改善成像与遗传学之间的弱联系，假设A₁是MRI-SNP相互作用邻接矩阵，A₂是MRI-Gene相互作用邻接矩阵，而A₃是SNP-Gene相互作用邻接矩阵。另外，为了稀疏数据以发现数据的关键特征，JSNMNMF使用的方法来控制W和H的稀疏性。因此，其目标功能如下：

其中，参数λ₁，λ₂和λ₃是约束的权重，γ₁用于限制W，γ₂用于限制H。

3、联合连通性的稀疏非负矩阵分解(JCB-SNMF)

为了鼓励规范向量相关元素的相似性，本专利在JSNMNMF算法中引入了基于连通性的惩罚项。具体来说，如果第i个节点和第j个节点(即大脑区域或SNP位点)之间的连通性很高，它将迫使范数向量的相应元素相似。因此，将大脑的连通性信息和加权的SNP相关网络相结合，以捕获遗传网络结构，作为要添加到算法中的先验矩阵，旨在提高提取特征的生物学意义。

其中，

和

分别代表X₁、X₂和X₃的拉普拉斯矩阵，把他们作为新的约束项加入到本专利所提出算法中，并改为如下形式：

P(H₁)＝Tr(H₁ ^TB₁H₁)

P(H₂)＝Tr(H₂TB₂H₂)

P(H₃)＝Tr(H₃TB₃H₃)

其中，B₁，B₂和B₃分别代表X₁、X₂和X₃的拉普拉斯矩阵。

现在，所提出算法的目标函数如下：

令

和

(H_I)_ij≥0的拉格朗日乘数，

L关于W和H_I的偏导数是：

基于Karush-Kuhn-Tucher(KKT)条件，Ψ_ijW_ij为0且

为0，可以得到W_ij和(H_I)_ij：

图1是本发明提供的一种识别阿尔茨海默症相关生物标志物的相关性分析方法的应用流程图。具体来说，将三种数据(结构磁共振成像的ROI、单核苷酸多态性、基因表达数据)放入到本专利所提出的算法中，得到若干个共表达模块，筛选其中最显著模块(p＜0.05)后，对模块中的ROI、SNPs和基因做生物学分析，包括KEGG富集分析、PPI网络构建、对PPI网络中相互作用最多的基因进行ROC曲线分析。最后找出模块中显著的ROI-SNP和ROI-gene对。

在本发明的一个实施例中，使用阿尔兹海默神经影像计划(The Alzheime’sDisease Neuroimaging Initiative，ADNI)数据库中的样本验证本发明提供的一种识别阿尔茨海默症相关生物标志物的相关性分析方法的有效性。在这个数据库中，有180名具有成像和基因分型数据的非西班牙裔白种人参与者，包括21名健康对照(HC)，147名轻度认知障碍(MCI)和12名AD患者。

在实验中，从ADNI1下载了原始MRI图像，使用DiffusionKit软件实现了原始MRI图像的头部运动校正。在蒙特利尔神经病学研究所(MNI)的标准空间进行了配准。接下来，通过使用SPM软件包中的MATLAB软件CAT工具包来实现MRI图像的分割。具体而言，基于体素的形态计量学(VBM)提供了特定组织隔室的局部数量或体积的体素估计，通过缩放由于空间配准而引起的体积变化以调整分割，可以计算感兴趣区域(ROI)中的灰质组织的体积为特征。筛选后，最终保留了140个ROI。所有SNP均由人类610 Quad BeadChip进行了基因分型。遗传分析工具PLINK使用以下排除标准筛选基因型数据：稀有SNP(次等位基因频率(MAF)<0.05)，违反Hardy-Weinberg平衡(HWE p<10-6)，每个受试者和每个SNP标记的呼叫率(<90％)，性别检查和同级对识别。这导致最终的数据集跨越5947个SNP位点。最后，使用limma软件包筛选具有显着差异表达的基因。当p值高于0.01时，那些基因被去除，获得了1477个基因。

经过上述预处理后，在真实数据集对本发明进行了验证。ROI，SNP，GENE的原始矩阵和重构矩阵的皮尔森(Pearson)相关系数分别为0.9962、0.6019、0.9921，并获得了30个共表达模块。

为了验证本发明的相关性分析能力，对所有模块中的基因表达数据进行了KEGG富集分析。提取了七个最丰富的生物过程关键字，并计算了每个模块中涉及的关键字数量。从图2(纵坐标代表模块数量)可以看出，有19个算法选出的模块涉及所有七个关键字词(包括剪接体、非酒精性脂肪肝疾病、阿尔兹海默症、甲型流感、信号通路、利什曼病、弓形虫病)。Toll样受体可以通过激活NF-κB和其他导致促炎分子合成的转录因子激活促炎免疫反应，并在与神经发炎相关的疾病(例如AD，帕金森氏病等)中发挥重要作用。表明弓形虫病感染小鼠的免疫细胞通过细胞因子网络促进神经炎症，并增强AD小鼠的认知障碍。甲型流感也与AD密切相关。β-淀粉样蛋白(βA)的积累可引起AD恶化，但也可抑制A型流感病毒。因此，本发明可以有效地选择与疾病相关的模块，并且不同的模块具有一定程度的表示性。

选出了p<0.05的5个模块进行下一步分析，如下表一所示。其中，模块6和模块10包含的SNP超过SNP总数的60％，因此不予考虑。为其他模块中的三个数据特征绘制维恩图，并比较不同模块中三个数据特征的逃逸率。如图3A～图3C所示，模块1的逃逸率是三种数据中最小的，选择模块1进行进一步分析。

表一模块选择

表二列出了本发明从SNPs基因座中识别出的ROI和基因。如表二所示，脑白质、角状回、梭状回、枕下回、颞下回、舌状回、额中回、中颞回、前丘脑、后中央回、前中央回、上颌回、上顶叶、上颞回由本发明提供的方法来选取，已有研究证实了他们是AD的风险ROI。此外，从第一个模块的SNPs数据中识别出总共47个风险基因。其中，DHX57和SPON1与小脑年龄加快有关，分别与AD有显着关系。EOMES和RGS6也被已有文献确定为AD风险基因。PPP2R2C的表达异常可能与AD的发作有关。淀粉样蛋白β肽前体淀粉样蛋白前体蛋白(APP)在AD中起着核心作用，而CASK是APP细胞内结构域(AICD)的相互作用者。在tauopathy小鼠模型中，PCSK2作为其海马中的差异基因被上调。CNTNAP5来自接触蛋白相关蛋白(Caspr)家族，与多种神经退行性疾病有关。此外，SLC9A9也与神经精神疾病有关。在模块1中发现了几个与老年性和神经性疾病相关的SNP位点。脆弱是衰老的复杂表型。研究发现rs1324192上的风险等位基因单倍型的患病率明显高于非弱势老年人。内含子SNP rs3802890与女性自闭症特征有关。rs981975也已被发现是与抗精神分裂症剂量变化相关的最高SNP。此外，已有的机器学习方法验证了rs12185438与帕金森综合症极为相关。

表二模块1中的选定ROI和基因

此外，从基因表达数据集中选择了258个基因。如图4所示，对这些基因进行了GO富集分析。核骨架的功能障碍已被证实是阿尔茨海默氏病相关的神经变性的原因。poly(A)结合蛋白MSUT2的活性决定了哺乳动物脑中对tau蛋白的易感性。大多数其他丰富的生物过程也与AD密切相关。我们还为模块1中选择的基因构建了蛋白质-蛋白质相互作用网络(PPI)。如图5所示，我们保留了177个具有相互作用关系的基因。我们选择了十个相互作用最强的基因，即剪接因子3b亚基1(SF3B1)，溴结构域PHD手指转录因子(BPTF)，收费像受体4(TLR4)，雷帕霉素激酶的机制靶标(MTOR)，LSM3同源物U6小核RNA和mRNA降解相关(LSM3)，丝氨酸和精氨酸重复矩阵1(SRRM1)，聚(rC)结合蛋白1(PCBP1)，异质核糖核蛋白A3(HNRNPA3)，核糖体蛋白S20(RPS20)和RNA结合基序蛋白14(RBM14)。其中，已有研究证实在初始阶段，TLR4的激活对淀粉样蛋白β(Aβ)具有有效的清除作用。但是，长期激活状态将导致Aβ沉积在大脑中。可以通过调节MTOR信号通路来预防和治疗AD。LSM3是AD的主要致病基因，也是模块网络中与MCI和AD密切相关的核心基因。HNRNPA3产品的生物学功能与炎症和神经退行性疾病有关。

使用IBM SPSS Statistics 22通过接收器工作特征(ROC)曲线分析了这些重要基因。在图6和表三中显示了曲线下面积(AUC)值大于0.5和P值小于0.05。本发明提供的方法不仅验证了模块1中的几个AD和MCI风险基因，还发现了可能与AD和MCI相关的三个潜在基因。

表三四个重要基因的ROC曲线信息

Gene ID	AUC	P value	95％CI
				SF3B1	0.656	0.026	0.525-0.788
RPS20	0.641	0.045	0.504-0.777
				RBM14	0.691	0.007	0.577-0.805

根据在模块1中选择的SNP，基因和大脑ROI，图7显示了大脑ROI-SNP/基因对的成对相关热图。不出所料，大多数ROI-SNP/基因对都很强。为了发现明显更强的SNP/gene-ROI对，我们分别在表4中显示了前10对，p<0.01。从表四中可以看出rs11918049和KLHL8可能与多个大脑区域的灰质体积变化密切相关。

表四在模块1中，前十对SNP-ROI和Gene-ROI分别为p<0.01。

在本发明中，为了证实算法具有更好的抗噪性能，另外随机生成两个数据集，分别在两个数据集上进行了较大和较小的样本量和特征实验。用n代表样本数量，p₁代表SMRI特征，q₁代表SNP特征，q₂代表基因特征。在大集合中，分别设置n＝300，p₁＝2000，q₁＝2500，q₂＝1000。在小集合中，分别设置n＝100，p₁＝650，q₁＝350，q₂＝600。基矩阵W和系数矩阵H₁，H₂和H₃中的元素都是随机整数。W，H₁，H₂和H₃由等式α[n]＝{β_i|β_i＝α+lη_i(i＝1,2,…,n)}生成。α是由均匀分布U(1,10)中的随机整数组成的矩阵。然后，使用η_i代表高斯噪声，使用l代表噪声水平。从图8可以看出，在两个数据集中增加噪声的过程中，所提算法的重构误差和目标函数值均小于其他两种算法。

本发明提出了JCB-SNMF方法，该方法考虑了大脑和基因各自的连通性信息，将拉普拉斯(Laplacian)矩阵添加为JSNMNMF算法的先验知识，从而提高了该算法的抗噪性能和生物学解释性。本发明使用ADNI数据集中的SMRI，SNP和基因表达数据，仿真结果表明，本发明的抗噪声性能优于JSNMNMF和JNMF。真实数据的实验结果表明，本发明可以识别和预测与AD和MCI密切相关的风险ROI，风险SNP和风险基因。不仅如此，还发现了一些重要的SNP/Gene-ROI对，值得进行更深入的研究。

需要说明的是，在本发明的实施例中，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“顺时针”、“逆时针”、“轴向”、“径向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述实施例，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (2)

1.一种识别阿尔茨海默症相关生物标志物的相关性分析方法，其特征在于，基于联合连通性的稀疏非负矩阵分解，将结构磁共振成像SMRI、单核苷酸多态性位点SNP和基因表达数据同时投影到一个公共特征空间中，还添加了大脑每个区域的连通性信息和遗传数据作为先验知识，以识别与阿尔茨海默病患者、早期认知障碍患者相关的感兴趣区域ROI，风险SNP位点和风险基因。

2.如权利要求1所述的识别阿尔茨海默症相关生物标志物的相关性分析方法，其特征在于，所述相关性分析方法的和目标函数如下：

令

和

为W_ij≥0、(H_I)_ij≥0的拉格朗日乘数，

L关于W和H_I的偏导数是：

基于Karush-Kuhn-Tucher条件，Ψ_ijW_ij为0，且

(H_I)_ij为0，得到W_ij和(H_I)_ij：

其中，X₁代表输入的sMRI影像数据，X₂代表SNP数据，X₃代表基因表达数据；W是三种数据的基矩阵，H_i(i＝1,2,3…)分别是三种数据的系数矩阵，A_i(i＝1,2,3…)代表不同数据之间的邻接矩阵，B_i(i＝1,2,3…)代表每种数据的拉普拉斯矩阵；λ_i、β_i、γ_i(i＝1,2,3…)是控制各约束项稀疏性的参数；i,j分别代表矩阵的第i行和第j列；k代表共表达模块数。

Images