阿尔茨海默症相关基因的预测方法与流程

本发明具体涉及一种阿尔茨海默症相关基因的预测方法。

背景技术：

阿尔茨海默病(ad)在美国是痴呆症最常见的病因，其特点是各种认知功能的逐步下降。年龄是ad中最重要的风险因素，但认知功能的下降是个体特异性的，还受环境因素，个体经验和遗传因素的影响。与早期发病相关的基因通常为突变的淀粉样前体蛋白基因app，早衰老性psen1与psen2基因。与晚发性ad相关的候选基因包括a2m，abca1/2，clu，picalm，sorl1和trem2。apoe，app，tau也有与晚发性ad相关的可变剪切变体。了解不同的基因表达可以最终回答有关ad发病机理的问题，并确定可能的疾病治疗目标。

在寻找与疾病相关基因的过程中，传统的分析途径，如差异分析软件edger，deseq2，limma的voom函数等，是在基因水平上对表达数据进行分析。但是，现有的分析方法所找到的基因并不能与病症进行良好的对应，从而影响了人们对于病症的认识。

技术实现要素：

本发明的目的在于提供一种能够对阿尔茨海默症的相关基因进行预测的阿尔茨海默症相关基因的预测方法。

本发明提供的这种阿尔茨海默症相关基因的预测方法，包括如下步骤：

s1.获取ad患者和非ad患者的数据信息；

s2.计算ad患者和非ad患者的剪切异构体(isoform)表达的fpkm值，根据得到的fpkm值对剪切异构体进行筛选，并对fpkm值进行数据处理得到fpkm1值；

s3.根据步骤s2得到的fpkm1值，再次对剪切异构体进行筛选；

s4.以步骤s3得到的剪切异构体的信息为基准，对ad患者和非ad患者进行聚类，并剔除离群值；

s5.以步骤s4得到的ad患者和步骤s3得到的剪切异构体构建ad患者无尺度网络，同时以步骤s4得到的非ad患者和步骤s3得到的剪切异构体构建非ad患者无尺度网络；

s6.根据步骤s5得到的ad患者无尺度网络和非ad患者无尺度网络，分别构建各自的拓扑重叠矩阵tom-ad和tom-非ad；

s7.将步骤s6得到的拓扑重叠矩阵tom-ad和tom-非ad，各自转换为节点相异性矩阵disstom-ad和disstom-非ad，并对节点相异性矩阵disstom-ad和disstom-非ad各自进行聚类得到若干个模块；

s8.对步骤s7得到的若干个模块进行相似性评价并合并；

s9.将步骤s6得到的拓扑重叠矩阵tom-ad和tom-非ad各自转换为网络边文件；

s10.根据步骤s9得到的网络边文件生成无向网络，计算无向网络的网络参数，并根据网络参数筛选得到具有网络结构差异的剪切异构体；

s11.对步骤s8得到的若干个模块进行合并，并判断合并后模块的差异性，从而得到具有差异性的剪切异构体；

s12.将步骤s10得到具有网络结构差异的剪切异构体和步骤s11得到的具有差异性的剪切异构体取并集，得到最终预测的与阿尔茨海默症相关的剪切异构体。

步骤s1所述的ad患者和非ad患者的数据信息，具体包括ad患者的患病信息和剪切异构体在患者内的表达信息，以及非ad患者的患病信息和剪切异构体在患者内的表达信息。

步骤s2所述的根据得到的fpkm值对剪切异构体进行筛选，并对fpkm值进行数据处理得到fpkm1值，具体为采用如下步骤进行筛选和处理：

(1)采用如下规则进行剪切异构体的筛选：

对于所有ad患者，若剪切异构体i满足如下公式，则保留剪切异构体i，否则删除剪切异构体i：

式中m为ad患者的总数，num(fpkmiad＞1)为剪切异构体i在所有ad患者中的fpkm值大于1的次数；

对于所有非ad患者，若剪切异构体i满足如下公式，则保留剪切异构体i，否则删除剪切异构体i：

式中n为ad患者的总数，num(fpkmi＞1)为剪切异构体i在所有非ad患者中的fpkm值大于1的次数；

(2)对于所有ad患者和非ad患者，采用如下算式计算得到fpkm1值：

fpkm1iad＝log2(fpkmiad+1)

fpkm1i＝log2(fpkmi+1)

式中fpkmiad为ad患者的原始的第i个fpkm值，fpkm1iad为ad患者的处理后的第i个fpkm1i值；fpkmi为非ad患者的原始的第i个fpkm值，fpkm1i为非ad患者的处理后的第i个fpkm1i值。

步骤s3所述的根据fpkm1值再次对剪切异构体进行筛选，具体为采用如下步骤进行筛选：

1)对于ad患者，采用如下规则进行剪切异构体的筛选：

若剪切异构体i在所有ad患者中的fpkm1值的方差不为0，则保留剪切异构体i；若剪切异构体i在所有ad患者中的fpkm1值的方差为0，则删除剪切异构体i；

2)对于非ad患者，采用如下规则进行剪切异构体的筛选：

若剪切异构体i在所有非ad患者中的fpkm1值的方差不为0，则保留剪切异构体i；若剪切异构体i在所有非ad患者中的fpkm1值的方差为0，则删除剪切异构体i。

步骤s4所述的对ad患者和非ad患者进行聚类，具体为对ad患者以剪切异构体的信息为基准进行层次聚类，同时对非ad患者以剪切异构体的信息为基准进行层次聚类。

步骤s5所述的构建无尺度网络，具体为对剪切异构体对的相关系数进行次方的幂指数加权得到的邻接矩阵构建。根据无尺度网络原则确定最佳加权系数(最佳软阈值)，即：连接节点个数取对数log(k)和节点出现概率对数值log(p(k))之间的相关系数至少达到0.8。此处加权系数进行网格搜索的范围为12-20，间隔为2。

步骤s6所述的根据无尺度网络构建拓扑重叠矩阵，具体为根据步骤s5构建无尺度网络时的最佳加权系数(最佳软阈值)构建邻接矩阵，并将构建的邻接矩阵转换为拓扑重叠矩阵。

步骤s8所述的对步骤s7得到的若干个模块进行相似性评价并合并，具体为采用动态树切割算法对各个模块进行检测并去除离群值，然后对各个模块进行相似性评价和合并。

所述的对各个模块进行相似性评价，具体为对每个模块进行pca(principalcomponentanalysis,主成分分析)降维，然后计算降维后各个模块的右奇异向量的不相似度，并根据计算得到的右奇异向量的不相似度进行聚类。

所述的根据右奇异向量的不相似度进行聚类，具体为设定不相似度的阈值为0.25，并根据设定的不相似度的阈值进行聚类。

步骤s9所述的将拓扑重叠矩阵tom-ad和tom-非ad各自转换为网络边文件，具体为采用如下规则进行转换：

对于拓扑重叠矩阵tom-ad，保留矩阵中值大于设定值t的行和列得到中转矩阵，将中转矩阵中行名称作为网络边文件的第一列，中转矩阵中的列名称作为网络边文件的第二列，将中转矩阵中的值作为网络边文件的第三列，从而得到拓扑重叠矩阵tom-ad的网络边文件；

对于拓扑重叠矩阵tom-非ad，保留矩阵中值大于设定值t的行和列得到中转矩阵，将中转矩阵中行名称作为网络边文件的第一列，中转矩阵中的列名称作为网络边文件的第二列，将中转矩阵中的值作为网络边文件的第三列，从而得到拓扑重叠矩阵tom-非ad的网络边文件。

所述的设定值t为0.4。

步骤10所述的计算无向网络的网络参数并筛选得到具有网络结构差异的剪切异构体，具体为采用如下步骤进行计算和筛选：

针对ad患者对应的无向网络，以及非ad患者对应的无向网络，分别采用如下步骤进行计算和筛选：

ⅰ.计算每一个剪切异构体的连接度degree和聚集系数cc；

ⅱ.删除连接度degree为0或聚集系数cc为0的剪切异构体；

ⅲ.根据剩下的剪切异构体的连接度degree计算剩余的每一个剪切异构体的fc值fcdegree，同时根据剩下的剪切异构体的聚集系数cc计算剩余的每一个剪切异构体的fc值fccc；

ⅳ.保留|log2fcdegree|＞1且|log2fccc|＞1的剪切异构体，从而得到具有网络结构差异的剪切异构体。

步骤s11所述的对步骤s7得到的若干个模块进行合并，判断合并后模块的差异性从而得到具有差异性的剪切异构体，具体为采用如下步骤进行合并和判断：

ⅰ.对于步骤s8得到的ad患者对应的若干个模块，以及非ad患者对应的若干个模块，采用如下规则进行合并：

对于任意的两个模块，若该两个模块共有的转录体的数目大于该两个模块中具有较少数目的转录体的模块的转录体数目的p％，则将该两个模块合并；所述的两个模块，为两个模块均为ad患者对应的模块，或者均为非ad患者对应的模块，或者一个属于ad患者对应的模块且另一个属于非ad患者对应的模块；p为正实数；

ⅱ.针对步骤ⅰ得到的合并后的模块，采用如下公式计算各个模块的差异性指标d：

d＝1-(e(inter)/e(sub_ad)+e(inter)/e(sun_con))/2

式中，e()为计算网络的边的数目的函数，sub_ad为ad患者模块在步骤s10所生成总网络中提取到的子网络，sun_con为非ad患者模块在步骤s10所生成总网络中提取到的子网络，inter为从sub_ad和sun_con中提取出的公共网络；

ⅲ.对每个模块，若模块的差异性指标大于事先设定的阈值q，则认定该模块具有差异性，该模块内所包含的剪切异构体为具有差异性的剪切异构体。

所述的p为80。

本发明提供的这种阿尔茨海默症相关基因的预测方法，基于剪切异构体层面的差异性表达分析，同时结合模块检测以及网络构建，综合考虑硬阈值转化相似矩阵的缺陷，提出了软阈值的概念，同时为了减少噪声和假关联的影响，还提出了采用拓扑重叠矩阵tom来替代邻接矩阵构建网络，并假设构建网络的前提是基因间相互作用关系满足无尺度网络原则，从而使得本发明构建的网络更近似于生物网络；因此，本发明方法能够对阿尔茨海默症的相关基因和剪切异构体进行预测，为后续的工作提供数据支撑和理论基础。

附图说明

图1为本发明方法的方法流程图。

具体实施方式

如图1所示为本发明方法的方法流程图：本发明提供的这种阿尔茨海默症相关基因的预测方法，包括如下步骤：

s1.获取ad患者和非ad患者的数据信息；具体包括ad患者的患病信息和剪切异构体在患者内的表达信息，以及非ad患者的患病信息和剪切异构体在患者内的表达信息；

s2.计算ad患者和非ad患者的剪切异构体(isoform)表达的fpkm值，根据得到的fpkm值对剪切异构体进行筛选，并对fpkm值进行数据处理得到fpkm1值；具体为采用如下步骤进行筛选和处理：

(1)采用如下规则进行剪切异构体的筛选：

对于所有ad患者，若剪切异构体i满足如下公式，则保留剪切异构体i，否则删除剪切异构体i：

式中m为ad患者的总数，num(fpkmiad＞1)为剪切异构体i在所有ad患者中的fpkm值大于1的次数；

对于所有非ad患者，若剪切异构体i满足如下公式，则保留剪切异构体i，否则删除剪切异构体i：

式中n为ad患者的总数，num(fpkmi＞1)为剪切异构体i在所有非ad患者中的fpkm值大于1的次数；

(2)对于所有ad患者和非ad患者，采用如下算式计算得到fpkm1值：

fpkm1iad＝log2(fpkmiad+1)

fpkm1i＝log2(fpkmi+1)

式中fpkmiad为ad患者的原始的第i个fpkm值，fpkm1iad为ad患者的处理后的第i个fpkm1i值；fpkmi为非ad患者的原始的第i个fpkm值，fpkm1i为非ad患者的处理后的第i个fpkm1i值；

s3.根据步骤s2得到的fpkm1值，再次对剪切异构体进行筛选；具体为采用如下步骤进行筛选：

1)对于ad患者，采用如下规则进行剪切异构体的筛选：

若剪切异构体i在所有ad患者中的fpkm1值的方差不为0，则保留剪切异构体i；若剪切异构体i在所有ad患者中的fpkm1值的方差为0，则删除剪切异构体i；

2)对于非ad患者，采用如下规则进行剪切异构体的筛选：

若剪切异构体i在所有非ad患者中的fpkm1值的方差不为0，则保留剪切异构体i；若剪切异构体i在所有非ad患者中的fpkm1值的方差为0，则删除剪切异构体i；

s4.以步骤s3得到的剪切异构体的信息为基准，对ad患者和非ad患者进行聚类，并剔除离群值；具体为对ad患者以剪切异构体的信息为基准进行层次聚类，同时对非ad患者以剪切异构体的信息为基准进行层次聚类；

s5.以步骤s4得到的ad患者和步骤s3得到的剪切异构体构建ad患者无尺度网络，同时以步骤s4得到的非ad患者和步骤s3得到的剪切异构体构建非ad患者无尺度网络；其中，无尺度网络依据对剪切异构体对的相关系数进行次方的幂指数加权得到的邻接矩阵构建；根据无尺度网络原则确定最佳加权系数(最佳软阈值)，即：连接节点个数取对数log(k)和节点出现概率对数值log(p(k))之间的相关系数至少达到0.8；此处加权系数进行网格搜索的范围为12-20，间隔为2；

s6.根据步骤s5得到的ad患者无尺度网络和非ad患者无尺度网络，分别根据各自的无尺度网络构建时的最佳加权系数(最佳软阈值)构建各自的邻接矩阵，并将各自构建的邻接矩阵转换为各自的拓扑重叠矩阵；

s7.将步骤s6得到的拓扑重叠矩阵tom-ad和tom-非ad，各自转换为节点相异性矩阵disstom-ad和disstom-非ad，并对节点相异性矩阵disstom-ad和disstom-非ad各自进行聚类得到若干个模块；

s8.对步骤s7得到的若干个模块，采用动态树切割算法对各个模块进行检测并去除离群值，然后对各个模块进行相似性评价和合并；

同时，对各个模块进行相似性评价，具体为对每个模块进行pca(主成分分析)降维，然后计算降维后各个模块的右奇异向量的不相似度，并根据计算得到的右奇异向量的不相似度，设定不相似度的阈值为0.25，并根据设定的不相似度的阈值进行聚类；

s9.将步骤s6得到的拓扑重叠矩阵tom-ad和tom-非ad各自转换为网络边文件；具体为采用如下规则进行转换：

对于拓扑重叠矩阵tom-ad，保留矩阵中值大于设定值t(优选为0.4)的行和列得到中转矩阵，将中转矩阵中行名称作为网络边文件的第一列，中转矩阵中的列名称作为网络边文件的第二列，将中转矩阵中的值作为网络边文件的第三列，从而得到拓扑重叠矩阵tom-ad的网络边文件；

对于拓扑重叠矩阵tom-非ad，保留矩阵中值大于设定值t(优选为0.4)的行和列得到中转矩阵，将中转矩阵中行名称作为网络边文件的第一列，中转矩阵中的列名称作为网络边文件的第二列，将中转矩阵中的值作为网络边文件的第三列，从而得到拓扑重叠矩阵tom-非ad的网络边文件；

s10.根据步骤s9得到的网络边文件生成无向网络，计算无向网络的网络参数，并根据网络参数筛选得到具有网络结构差异的剪切异构体；具体为采用如下步骤进行计算和筛选：

针对ad患者对应的无向网络，以及非ad患者对应的无向网络，分别采用如下步骤进行计算和筛选：

ⅰ.计算每一个剪切异构体的连接度degree和聚集系数cc；

ⅱ.删除连接度degree为0或聚集系数cc为0的剪切异构体；

ⅲ.根据剩下的剪切异构体的连接度degree计算剩余的每一个剪切异构体的fc值fcdegree，同时根据剩下的剪切异构体的聚集系数cc计算剩余的每一个剪切异构体的fc值fccc；

ⅳ.保留|log2fcdegree|＞1且|log2fccc|＞1的剪切异构体，从而得到具有网络结构差异的剪切异构体；

s11.对步骤s8得到的若干个模块进行合并，并判断合并后模块的差异性，从而得到具有差异性的剪切异构体；具体为采用如下步骤进行合并和判断：

ⅰ.对于步骤s8得到的ad患者对应的若干个模块，以及非ad患者对应的若干个模块，采用如下规则进行合并：

对于任意的两个模块，若该两个模块共有的转录体的数目大于该两个模块中具有较少数目的转录体的模块的转录体数目的p％(优选为80％)，则将该两个模块合并；所述的两个模块，为两个模块均为ad患者对应的模块，或者均为非ad患者对应的模块，或者一个属于ad患者对应的模块且另一个属于非ad患者对应的模块；

ⅱ.针对步骤ⅰ得到的合并后的模块，采用如下公式计算各个模块的差异性指标d：

d＝1-(e(inter)/e(sub_ad)+e(inter)/e(sun_con))/2

式中，e()为计算网络的边的数目的函数，sub_ad为ad患者模块在步骤s10所生成总网络中提取到的子网络，sun_con为非ad患者模块在步骤s10所生成总网络中提取到的子网络，inter为从sub_ad和sun_con中提取出的公共网络；

ⅲ.对每个模块，若模块的差异性指标大于事先设定的阈值q，则认定该模块具有差异性，该模块内所包含的剪切异构体为具有差异性的剪切异构体；

s12.将步骤s10得到具有网络结构差异的剪切异构体和步骤s11得到的具有差异性的剪切异构体取并集，得到最终预测的与阿尔茨海默症相关的剪切异构体。

以下，结合具体的实施例，对本发明进行进一步说明：

首先获取基础数据；

然后计算fpkm指标，并做fpkm1＝log2(fpkm+1)处理，仅保留fpkm1>1的样本数大于总样本数80％的转录体行。依据样本分组信息，将ad患者和非ad患者对应的数据形成ad_fpkm和control_fpkm文件；

读取上述的ad_fpkm和control_fpkm文件并转置，检查缺失条目的数据以及零方差的转录体，并筛选和更新数据；

进行层次聚类，采用默认参数，作图检测是否有离群值，若有则去除离群点并更新数据框。设定待选加权系数(软阈值)范围为12～20，间隔为2，阈值范围可根据情况适当扩宽，得到使转录体间关系符合无尺度网络原则的最佳软阈值；

选用合适软阈值，构建邻接矩阵。为减小噪声和假关联的影响，将邻接关系转换为拓扑重叠矩阵(tom)，并转换成节点相异性矩阵disstom。对disstom进行聚类，并使用动态树切割模式识别检测并去除离群值。为了量化整个模块的共同表达相似性，对每个模块作pca降维。计算降维后各模块的右奇异向量的不相似度并进行聚类，设定不相似度阈值为0.25，进行模块合并，从而得到ad患者的模块为43个，control(非ad患者对应的模块)为11个；

以tom作为输入，保留边的权重阈值设为0.4，得到第一二列为源终节点，第三列为边权重的网络边文件；

网络边文件作为输入，生成无向网络。分别计算网络属性—连接度degree和聚集系数cc。合并ad患者和control(非ad患者)网络计算出来的网络属性，去除属性为na和0的转录体行并求取fc(foldchange)值，筛选出两个属性均满足|fc|＞1的转录体行，得到差异表达亚型，共488个；

对模块采用如下合并原则：两个模块共有转录体数大于较小模块转录体数的80％。根据原则进行模块合并(合并后共50个模块)，合并后对每个模块提取出其在ad患者和control(非ad患者)两个网络中的子网络。计算模块的差异性指标并筛选，得到差异表达亚型，共549个

同时，根据ad患者和非ad患者的数据信息，采用deseq2算法(算法参数为：fdr<0.05且|logfc|>1阈值)进行剪切异构体差异表达分析，从而得到显著差异表达的剪切异构体，作为对比结果1；

根据ad患者和非ad患者的数据信息，按照tmm(trimmedmeanofm-values)方法进行标准化，然后进行离散度检测，将离散度检测结果进行似然比检验(似然比检验参数：fdr<0.05且|logfc|>1阈值)，得到显著差异表达的剪切异构体，作为对比结果2；

通过网络fc得到的结果，共488个转录体，模块差异性得到的结果，共549个转录体，将其合并后共1033个转录体；

在isoform层面deseq2和edger的分析结果分别为1548个和1844个转录体，再合并共3178个转录体；

将转录体id转为相应的基因id，共1793个基因。

本发明结果与gene-centric的差异基因表达分析方法相比，共得到1043个新型基因；

本发明的效果验证：

go富集验证新型基因。

根据1043个新型基因，目前所得的是基因的ensemble，go富集所需为基因的geneid，因此需要从ncbi中下载ensemble与geneid的映射文件，将ensemble转换成geneid。使用y叔clusterprofiler包中的enrichgo，输入为geneid，参数orgdb为人类数据库org.hs.eg.db；ont指定goterm的三大分支(bp：生物过程，mf：分子功能，cc：细胞组成)，这里指定ont为bp；pvaluecutoff和qvaluecutoff分别为筛选goterm的pvalue和qvalue阈值，我们将其设为0.01，执行函数即可获得go富集结果。

ppi网络验证新型基因。

结合genecards，omim，文献共得到85个ad相关基因，从hprd(humanproteinreferencedatabase)筛选得到人类(taxid＝9606)中与阿尔茨海默病相关的全部ppi网络。将85个ad相关基因比对到ppi网络中，得到ppi_net，共有731条边。将新型差异表达基因比对到ppi_net，共有43个基因匹配成功，其中4个在85个ad相关基因中，39个在ppi_net中与85个ad相关基因有边相邻(即存在直接或间接关联)。