基于cFDR的缺血性卒中与低密度脂蛋白胆固醇关联SNPs的识别

毛相杰，温书景，徐 飞，高 盼，张 蔷，杨 斌，周乾宇，王 博，岳文娟，孙长青

1)郑州大学公共卫生学院社会医学与卫生事业管理学教研室 郑州 450001 2)洛阳市石化医院内一科 河南洛阳 471012 3)河南省教育厅 郑州450018

环境、遗传因素及两者交互作用在缺血性卒中(ischemic stroke，IS)病因学中发挥着重要作用[1]，其中遗传因素可解释IS疾病风险的35%～40%[2]。低密度脂蛋白胆固醇(low-density lipoprotein cholesterol，LDL-C)是IS和冠心病等心脑血管疾病的危险因素[3]，其遗传估计值为40%～80%[4]。遗传多效性是指单个基因或遗传变异与两个及以上表型相关联的现象。前期研究[5-6]显示LDL-C与IS存在共享遗传位点。全基因组关联研究(genome-wide association study，GWAS)对复杂疾病或性状遗传位点的识别只能解释一小部分的遗传度，单核苷酸多态性(single-nucleotide polymorphism，SNP)在更大程度上影响疾病发病风险，但目前尚缺乏识别SNP的可靠方法。本研究应用条件错误发现率(conditional false discovery rate，cFDR)方法[7]整合IS及LDL-C的GWAS数据，以增加统计效能，进而提高相关遗传位点的识别力，尤其是发现遗传多效性位点，旨在获得IS和LDL-C之间未知的共享生物机制，为疾病的预防提供新的基因治疗靶点和理论支撑。

1 材料与方法

1.1材料IS的GWAS数据来源于2015年美国神经学和卒中研究所SiGN(Stroke Genetics Network)项目组织的大型荟萃分析，下载网址http://cerebrovascularportal.org/informational/downloads，主要是由欧洲血统人群及少量非洲与西班牙血统人群组成，共包含16 851个病例和32 473个对照。LDL-C的GWAS数据来源于GLGC(Global Lipids Genetics Consortium)，下载网址http://lipidgenetics.org/，包含188 577个欧洲血统和7 898个非欧洲血统参与者[8-9]。SIGN和GLGC之间有3个重叠样本，重叠人数3 836(占总人数的1.56%)。使用Linux操作系统和R 3.4.3软件进行cFDR分析，使用代码包为cfdr。

1.2数据整理及质量控制提取并整理全基因组测序的SNPs相关汇总统计结果，包括SNPs名称及关联性检验P值信息。统计分析过程剔除了主要组织相容性复合体(major histocompatibility complex,MHC)区域(chr6:25652429-33368333)和染色体8p23.1区域(chr8:7242715-12483982)的SNPs，以避免基因组复杂区域内SNPs之间的连锁不平衡影响研究结果的可靠性[10]。基因组控制校正方法应用于原始GWAS和表型meta分析，以确保每个SNP的方差估计不会因人口异质性而被夸大。

1.3两性状间的多效性基因富集的评估分层Q-Q图以SNP与其中一个表型关联检验P值的累积分布函数代替理论经验累积分布函数。将SNP与疾病关联的观察值分位数定义为p值(y轴)，SNP与疾病关联的理论值分位数定义为q值(x轴)；负自然对数转换后，y轴表示为-log10(P)，x轴表示为-log10(q)。观察分布线与理论分布线(x=y)重合，表明不同性状间没有多效性基因富集；观察分布线(x轴)向左偏离表明存在多效性基因富集，偏离越大，表明富集程度越大。以主要性状名义上的-log10(P)为x轴，以主要性状基于条件性状的关联的富集倍数为y轴创建富集倍数图，通过零线向上偏移的程度来评估富集倍数。

1.4cFDR计算方法及两性状关联SNPs的筛选条件cFDR为后验概率：如果两个性状的P值都≤观察到的P值，那么一个给定的SNP对于第1个性状来说是无效SNP，表示为：FDR(P1|P2)=π0(P2)P1/F(P1|P2)，其中P1、P2分别是两个性状的P值，F(P1|P2)是条件性累积分布函数，π0(P2)是第2个性状的P值≤P2时第1个性状无义SNP的百分比，较小的cFDR值表示SNP与其中一种或两种性状关联。

为鉴定多效性基因，需计算conjunction cFDR值。对于任意一个两性状关联的SNP，假设SNP与两种性状关联的理论值都≤观察值，则conjunction cFDR为SNP与任意一种性状都没有关联的后验概率值，表示为：FDRIS&LDL-C=Max{FDR(IS|LDL-C), FDR(LDL-C|IS)}，较小的conjunction cFDR值表明SNP与两种性状关联，即SNP具有遗传多效性。

1.5多效性SNPs的鉴定条件曼哈顿图能直观显示出基因组上连锁不平衡区域内SNPs的染色体位置信息，若cFDR<0.05[即SNPs的-log10(cFDR)>1.3]，则该SNP被认为与研究的性状关联；若SNP的conjunction cFDR<0.05[即-log10(conjunction cFDR)>1.3]，则该SNP被认为与研究的两种性状关联，具有基因多效性。

1.6基因功能富集分析GO基因功能富集分析从分子功能、细胞成分和生物过程共3个方面定性描述IS和LDL-C关联基因在功能通路中的表达情况以及关联基因功能间的联系，分析网址http://geneontology.org/。

2 结果

2.1多效性基因富集的评估结果见图1。图1A显示，以LDL-C为条件性状，IS与LDL-C之间存在较大程度的基因富集；图1B显示，以IS为条件性状，LDL-C与IS间有更大程度的基因富集；图1C显示，以LDL-C为主要性状，基因富集大约有19倍的增加；图1D显示，以IS为主要性状，基因富集增加了约17倍。综上，IS与LDL-C具有相关性。

A、B：分别以LDL-C和IS为条件性状的分层Q-Q图；C、D：分别以LDL-C和IS为主要性状的基因富集倍数图

2.2IS关联的SNPs以IS为主要性状，发现9个SNPs与IS关联，其中2个SNPs已被先前的GWAS研究报道，其余7个被鉴定为新的SNPs位点。9个SNPs被注释到12个与IS相关的基因上，其中6个已被GWAS报道与IS相关，余6个属于新鉴定的基因。见表1。

表1 IS关联的SNPs

*：新发现SNPs或基因

2.3LDL-C关联的SNPs以LDL-C为主要性状，发现245个与LDL-C相关联的SNPs。其中19个SNPs已被GWAS报道，77个与其他LDL-C相关报道位点存在高度连锁不平衡(R2>0.6)，52个与血脂其他性状相关联，剩下的97个被鉴定为新的LDL-C关联SNPs。245个SNPs位于21个染色体上，被注释到241个基因上，其中103个基因已被GWAS报道与LDL-C有关，其余138个基因经鉴定为新基因。

2.4多效性SNPs位点见图2。如图2所示，红线为截断值水平。共发现8个SNPs与IS和LDL-C共同关联，2个(rs10774625和rs2238151)已被报道与两个性状共同关联，其余6个被鉴定为新的SNPs(rs10466588、rs11066301、rs4767293和rs630014被报道与LDL-C相关，rs9557951和rs7987982未发现与任何性状有关联)。8个SNPs被注释到11个多效性基因上，ATXN2、ABO与ALDH2被报道与IS和LDL-C都相关，其余8个被鉴定为新的多效性基因。见表2。

2.5功能富集分析将IS和LDL-C共同关联的基因进行GO功能富集分析，发现关联基因大多数与脂蛋白代谢过程、胆固醇平衡、大脑发育、脂肪酸代谢负调节等通路有关，可影响心脑血管疾病的发展过程。见表3。

图2 多效性位点鉴定的条件曼哈顿图

SNPs作用位置注释基因染色体SNP性状rs10466588∗基因间LOC101929011∗,BUD13∗chr11LDL-Crs10774625基因内ATXN2chr12LDL-C/ISrs11066301∗基因内PTPN11∗chr12LDL-Crs2238151基因内ALDH2chr12LDL-C/ISrs4767293∗基因间ERP29∗, NAA25∗chr12LDL-Crs9557951∗基因间METTL21EP∗,SLC10A2∗chr13-rs7987982∗基因内COL4A1∗chr13-rs630014∗基因内ABOchr9LDL-C

*：新发现SNPs或基因

表3 基因功能富集分析

3 讨论

孟德尔随机化研究[11]表明LDL-C的升高与IS风险增加存在因果关联。TIMD4-HAVCR1变异研究[12]显示该基因与LDL-C和IS风险增加都有关联。这些研究说明IS和LDL-C存在遗传多效性。本研究基于cFDR方法验证并发现8个遗传多效性SNPs，其注释基因ATXN2、ALDH2和ABO已被报道与IS和LDL-C均相关，LOC101929011和BUD13与LDL-C相关，PTPN11、NAA25和COL4A1与IS相关；而METTL21EP、SLC10A2和ERP29被报道与老年痴呆症和炎症性肠病有关联[13-14]，且以往的研究显示SLC10A2和ERP29可能在IS和LDL-C的发病机制上起重要作用。

SNP rs4767293位于NAA25和ERP29基因间的区域，且两个基因被报道与炎症性肠病有关联[13]。流行病学研究[15]表明炎症性肠病后期会导致动脉血管内皮功能障碍和血小板聚集，造成动脉粥样硬化。ERP29位于内质网中并能在各种组织和细胞类型中表达，包括N终端和C终端域，而C终端域包含一个新的螺旋折叠可以直接结合某些膜蛋白或疏水分泌蛋白[16]。ERP29通过结合跨膜蛋白参与上皮细胞的形成，调控上皮细胞间质转化，从而影响癌症的进展[17]。有研究[18]表明ERP29可能通过调节相关酶的分泌或者易位到脂筏来参与LDL的氧化。蛋白质交互网络分析[18]表明ERP29还与钙网蛋白、生物合成前体和髓过氧化物酶有关联。

SNP rs9557951位于METTL21EP和SLC10A2基因间区域，其注释基因均与老年痴呆症有关[14]。老年痴呆与IS具有共同的危险因素，包括吸烟、饮酒、糖尿病、高胆固醇血症等。Meta分析表明IS患者老年痴呆症的患病风险增加了1.6倍[19]。在老年痴呆和IS疾病进程中低灌注和认知能力下降具有相似性，且老年痴呆与小血管卒中存在共享的遗传关联基因(ATP5H、KCTD2、ICT1基因与两性状都有关)[20-21]。SLC10A2基因编码了一种钠/胆汁酸协同转运蛋白，对维持胆固醇体内平衡至关重要，该基因在突变之前已在家族性高胆固醇血症病例中被发现；SLC10A2也可在大脑中表达[14]。SLC10A2基因是人类胆汁酸转运体(ASBT)的重要成员，ASBT能够抑制血浆甘油三酯的浓度并增加高密度脂蛋白胆固醇的浓度[22]。最近研究[23]发现胆汁酸参与葡萄糖、脂质和能量代谢的调节。

本研究还存在一些局限性，首先GWAS汇总数据的重叠样本可能会造成假阳性，但cFDR应用贝叶斯统计、条件概率原理且属于非模型分析，可以降低或减少其造成的FDR。其次，该方法只能识别多效性位点而不能解决水平多效性和垂直多效性问题。

综上所述，cFDR方法通过对两种关联性状的GWAS汇总数据整合和重分析，验证并识别了9个IS关联位点和245个LDL-C关联位点，并初步鉴定出8个IS和LDL-C共享的多效性位点，结合实验验证可进一步探索关联疾病之间共同的致病机制，有助于发现早期诊断和治疗的新基因靶点。