本发明涉及一种筛选急性髓系白血病dna甲基化预后标志物的方法，属于生物医药领域。
背景技术：
白血病是我国最常见恶性肿瘤之一,在儿童及40岁以下年轻癌症患者中，白血病死亡率高居首位。2016年我国学者发表重要文章《cancerstatisticsinchina,2015》，揭示白血病的发病率排在第12位，死亡率则上升至第9位，每年新发白血病病例约5万人，且发病率和死亡率均呈现逐年上升的趋势[1]。在白血病中，以急性髓系白血病(acutemyeloidleukemia,aml)为主，约占急性白血病发病率的70％[2]。aml发病机制复杂、临床异质性大、复发死亡率高[3-7]，治疗费用昂贵，给家庭和社会带来沉重的经济负担。攻关aml的精准诊疗是改善疗效的必然选择，成为当下研究的热点和难点。
aml目前主要依据骨髓形态学(morphology，m)、免疫学(immunology，i)、细胞遗传学(cytogenetics，c)和分子生物学(molecular，m)(micm)特征进行诊断分型和预后分层[8,9]。基于目前的预后分层体系，仍有近50％的患者为正常核型且缺乏预后相关基因突变[10-12]，无法实现精准预后指导下的精准治疗。
研究表明，dna甲基化、组蛋白乙酰化、染色质重塑等表观遗传学修饰与aml的发生发展密切相关[13]。dna甲基化异常是aml的一个显著特征[13-15]，在指导aml的分型及预后分层方面具有重要作用[16-18]。然而，由于目前dna甲基化检测技术的不足，dna甲基化检测未能在临床上常规应用，如全基因组重亚硫酸盐测序(whole-genomebisulfitesequencing，wgbs)价格昂贵、测序深度不足、数据分析困难；甲基化450k芯片(illuminahumanmethylation450beadchiparray，illumina450karray)并非直接测定碱基、检测位点有限且数据分析系统过于封闭；agilentsureselect甲基化测序分子多样性低、起始样本量大、仅单链捕获；简化表观重亚硫酸盐测序(reducedrepresentationbisulfitesequencing，rrbs)受限于酶切位点、区域固定及存在数据丢失等[19,20]。因此，亟需一种全新的、适合于临床的dna甲基化检测方法在aml中推广应用，通过准确测定aml的dna甲基化状态，完善预后分层，推动精准诊疗。
2015年发表在naturecommunication杂志上的“全基因组cpg位点甲基化捕获测序技术(genome-cpgs-scalemethylc-capturesequencing，mcc-seq)”，是基于二代测序平台的全新的dna甲基化捕获检测技术[20]。此技术既弥补了450k芯片cpg位点覆盖不足的缺陷，又解决了wgbs覆盖面过大、测序深度不足、位点精准度不够、耗时费力成本高的问题，并具有独特探针设计、能够双链捕获、起始样本量低、准确定量检测等优点，是一种准确、经济、高效的dna甲基化检测方法，将促进dna甲基化研究的推广和应用。迄今为止，在aml中，尚无开展mcc-seq测序技术的相关研究报道，更无应用此技术进行aml预后研究的报道。
综上，鉴于aml异质性大、目前预后分层体系尚有不足、dna甲基化与aml密切相关、既往甲基化研究方法存在缺陷、mcc-seq检测甲基化具有多重优势，我们有充分依据认为，mcc-seq可用于检测aml的dna甲基化，精确描绘aml的甲基化特征，通过dna甲基化差异分析，寻找与预后不良相关的甲基化调控基因，完善aml预后分层体系，指导aml的精准诊疗，从而提高疗效，改善患者预后，为国家和社会减轻经济负担，更为推动精准医学的发展做出贡献。
全基因组cpg位点甲基化捕获测序技术(genome-widemethylc-capturesequencing，mcc-seq)是基于罗氏nimblegenseqcapepi序列富集系统的一种全新的dna甲基化检测技术，覆盖了550万个甲基化位点，还加入了cpg岛以外的cpg位点及疾病相关区域的位点。相比于illumina450k芯片，seqcapepi包括了dna芯片上的所有位点，并多出12倍以上的cpg位点，并且公共位点的甲基化水平分析结果一致性高。此外，还可以发现甲基化修饰改变的新位点、获得等位基因特异性的甲基化模式[21,22]。mcc-seq检测方法可以有效的避免样品损失，样品起始量仅1μg，且关键的测序指标有30％的改善，如均一性，这反映了文库复杂度的改善。另外，其针对一个区段同一个位置设计多种探针来保证不同情况的甲基化都能被有效捕获，并针对转化后的正反双链，均设计了探针，通过正负链的测序结果作为参考，以不发生甲基化的链的测序结果作为参考，来判断胞嘧啶(c)转变为胸腺嘧啶(t)是基因组发生的突变还是由于未甲基化的胞嘧啶经重亚硫酸盐处理的结果，从而提高甲基化分析的正确性。
2015年，naturecommunication杂志首先报道了mcc-seq用于人类脂肪组织表观基因组关联分析的研究(epigenome-wideassociationstudy，ewas)，证实其是一种准确、经济、高效的dna甲基化检测方法，尤其适用于疾病易感性/表型相关功能性差异甲基化基因的研究[20]。
技术实现要素：
本发明所要解决的技术问题是：鉴于aml发病率及死亡率高、异质性大、目前预后分层体系尚有不足、dna甲基化与aml密切相关、既往甲基化研究方法存在缺陷、mcc-seq检测甲基化具有多重优势，本发明应用mcc-seq技术，对aml骨髓标本进行dna甲基化测序，定量检测aml全基因组cpg位点甲基化状态，并通过dna甲基化差异分析，寻找与预后不良相关的甲基化调控基因，指导aml的精准诊疗，改善患者预后，为国家和社会减轻经济负担，更为推动精准医学的发展做出贡献。
本发明提供的技术方案是：一种筛选急性髓系白血病dna甲基化预后标志物的方法，包括如下步骤：
(1)采集样本和收集临床资料；
(2)制备样本；
(3)提取基因组dna；
(4)全基因组cpg位点甲基化捕获测序，步骤如下：
a)首先，构建文库：进行dna的片段化，末端修复，加“a”，接头连接，bisulfite处理，pcr扩增及产物纯化，最后得到纯化产物；
b)目标区域捕获测序：杂交前文库与芯片温浴杂交，目的dna片段洗脱，目的dna片段的pcr扩增，测序；
(5)甲基化测序数据分析：
首先，进行数据过滤及比对：将步骤(4)中的测序数据进行过滤，去掉低质量数据，得到可用数据，数据检测合格后，将可用数据与参考基因组进行比对，得到比对结果，在确认比对质量合格后，使用唯一比对数据得到c碱基甲基化信息，进行信息分析处理，得到标准信息分析结果和个性化分析结果；
然后，捕获区域的深度和覆盖度分析：分析不同reads测序深度下的捕获区域的覆盖度，并分析c碱基有效测序深度的累积分布，甲基化c碱基在基因组上的分布包含三种形式即cg，chg和chh，其中h代表a或t或c碱基；
进行甲基化位点识别及计算甲基化水平，计算各类型mc的位点数目及其在全部mc的位点中所占的比例，并计算甲基化水平；
分析差异甲基化区域：在多个不同分组间、不同基因区域寻找符合条件的差异甲基化区域dmrs，将找到的dmrs注释到基因体，比较dmrs在基因体不同区域分布情况；
筛选差异甲基化基因：对差异甲基化基因进行go富集性分析和kegg富集性分析，确定差异甲基化基因所参与的主要生化代谢途径及信号转导途径，最后筛选出一组12个与遗传学预后相关、受dna甲基化调控的功能性差异甲基化基因，即：bard1，bcl9l，clec11a，defb1，foxd2，igf1，il18，itih1，lsp1，p2rx6，rnase3和tubgcp2。
所述的方法，进一步包括步骤(6)，即：临床意义验证，根据筛选到的12个dmgs，计算每一例患者每一个dmgs的dna甲基化水平，同一患者12个dmgs甲基化水平的平均值记为该患者的综合甲基化值m-value，用于代表患者的甲基化水平。
所述的方法，步骤(5)中，每一个甲基化c碱基的甲基化水平均按如下公式进行计算：c位点的甲基化水平＝支持甲基化的reads的数目/(支持甲基化的reads的数目+支持非甲基化的reads的数目)×100％；
不同基因区域的平均甲基化水平计算公式如下：
某区域平均甲基化水平＝该区域支持甲基化的reads的数目/该区域总的reads的数目×100％。
步骤(5)中，dmr符合以下三个条件：至少包含5个cpgs；其中至少有3个差异甲基化位点，区域差异甲基化水平≥20％。
步骤b)中，扩增目的dna片段的qpcr引物1，qpcr引物2，
引物1序列：上游：gttaggtagggaagaagggagtagt
下游：cccaaaaatcaaataatcaaaaaaa
引物2序列：上游：gtggttaattaatttttgagttttgt
下游：tattaccctataaccaccatcacc步骤(7)中，数据过滤的处理步骤如下：
(1)去除接头污染的reads，即：去掉reads中接头污染的碱基数大于5bp；
(2)去除低质量的reads；
(3)去除含n比例大于5％的reads。
所述的方法，步骤(5)中，可用数据与参考基因组进行比对，在信息分析过程中，测序的结果和参考基因组都进行了c-to-t(forward)和g-to-a(reverse)的转换，将转换后的reads和转换后的基因组序列进行比对，唯一比对reads将用于甲基化信息的分析，用于比对的基因组数据库位universityofcaliforniasantacruzhg19。
本发明还提供所述的方法筛选得到的急性髓系白血病dna甲基化预后标志物，即dna甲基化调控的功能性差异甲基化基因：bard1，bcl9l，clec11a，defb1，foxd2，igf1，il18，itih1，lsp1，p2rx6，rnase3和tubgcp2。
同时，本发明还提供一种一种用于筛选急性髓系白血病dna甲基化预后标志物的试剂盒，该试剂盒包括：dna末端修复体系：该体系包含有10xkapaendrepairbuffer7μl，kapaendrepairenzyme5μl，ddh2o8μl，总体积70ul；加“a”体系：该体系包含有10xkapaa-tailingbuffer5μl，apaa-tailingenzyme3μl，总体积50ul；连接体系：该体系包含有5xkapaligationbuffer10μl，kapat4dnaligase5μl，methylationadapter(10μm)5μl，总体积50ul；pcr扩增体系：该体系包含有kapahifihotstartreadymix25μl，kapalibrary扩增引物2μl，去离子水13μl，总体积50ul，其中，扩增引物序列：上游：agtggttaattaattttcgagtttc，下游：tattattaccctataaccaccatcg；目的片段捕获后的pcr扩增体系：该体系包含有qpcr引物1和qpcr引物2各1μl，hifidna多聚酶混合物12.5μl，总量25μl，其中，qpcr引物1序列：上游：gttaggtagggaagaagggagtagt，下游：cccaaaaatcaaataatcaaaaaaa，qpcr引物2序列：上游：gtggttaattaatttttgagttttgt，下游：tattaccctataaccaccatcacc。
本发明具有以下有益效果：
本发明主要是从临床的实际角度出发，以解决临床问题为导向，应用全基因组cpg位点甲基化捕获测序技术检测急性髓系白血病dna甲基化并进行预后分析。本发明方法具有良好的转换率(>99.5％)、q30(95.03％)、比对率(92.90％)及经技术重复样本验证具有很好的稳定性。本发明通过整合遗传学和表观遗传学信息，筛选到了具有预后指导意义的受dna甲基化调控的基因标志物，有助于完善aml预后分层体系，指导aml的精准诊疗，从而提高疗效，改善患者预后。本发明应用mcc-seq技术检测aml骨髓样本数据准确、可靠，具有良好的可行性和实用性。aml患者dna甲基化水平是一个评价甲基化状态的稳定、可靠的指标，aml患者启动子区具有最重要的功能性差异甲基化区域，适用于筛选预后相关的dna甲基化调控基因。因此，本发明筛选得到的12个dmgs的平均甲基化水平m-value是评价患者dna甲基化水平的稳定、可靠的指标，12个dmgs可作为aml预后分层的新的生物标志物，高m-value是诱导缓解、无病生存和总生存的预后不良因素。
附图说明
图1全基因组cpg位点甲基化水平与各临床因素的关系。
图2基因启动子区具有最重要的功能性dna甲基化特征，其中：图a，初治aml和正常对照骨髓样本不同基因区域的dna甲基化水平的比较；图b，8对初治-缓解配对aml骨髓样本不同基因功能区域的dna甲基化水平的比较；图c，初治aml和正常对照骨髓样本比较不同基因功能区域的差异甲基化区域；图d，初治aml骨髓样本不同遗传学预后分组间比较的不同基因功能区域的差异甲基化基因。
具体实施方式
下面通过具体实施方式的详细描述来进一步阐明本发明，但并不是对本发明的限制，仅仅作示例说明。
实施例中如无特殊说明，所采用的实验方法，均为常规方法；所用的实验材料、试剂等，均可通过商业途径得到。
一、实验材料
1、mcc-seq测序病例
本发明选取aml患者，按照纳入标准：(1)明确诊断为非m3型急性髓系白血病；(2)年龄大于18岁；(3)至少完成一疗程化疗且评价疗效者；(4)micm诊断资料、临床治疗资料完整；(5)有突变靶向测序结果(以便准确进行分子生物学预后分层)；(6)至少保存有初治时的骨髓标本。最终选取21例初治aml患者(编号为c01，c02，c03……c21)，共检测44例次骨髓样本，详细信息见表1。
表1dna甲基化测序患者及样本信息总结表
“7+3”方案,以阿糖胞苷和蒽环类抗生素为基础的标准诱导方案，具体参见nccn指南；地西他滨联合方案，地西他滨20mg/m2d1-5,阿糖胞苷10mg/m2q12hd1-5,阿克拉霉素20mgd1,3,5,粒细胞刺激因子300μg/d从化疗前一天用至中性粒细胞恢复。
2、tcga数据库验证病例
为了更好地验证从自测病例中筛选到的差异甲基化基因的临床意义，本发明人从tcga数据库调取了illumina450k芯片检测甲基化的aml患者的数据(https://tcga-data.nci.nih.gov/tcga/)[23]。根据临床信息及dna甲基化检测数据的完整性进行筛选，最终共有169例aml患者纳入本研究。患者信息见表2。
表2tcga数据库患者信息总结表
二、实验方法
1、样本采集及临床资料收集
所有病例(21例初治aml、2例复发aml、5例健康对照)签署知情同意书后，通过骨髓穿刺采集3-5ml骨髓样本。收集患者的临床诊治信息，包括一般信息、骨髓micm信息、aml诊断、治疗及疗效、生存情况等。
2、技术重复样本制备
我们随机选取了2例复发的aml患者(c22和c23)，分别采集新鲜骨髓标本并提取单个核细胞，然后将同一患者的单个核细胞均分为2份，从而得到4份技术重复样本(s22-重复1和s22-重复2，s23-重复1和s23-重复2)，用于提取dna和mcc-seq测序，以评估mcc-seq检测骨髓样本的稳定性。
3、基因组dna提取
依据genomidnaextractionkit(promega)的说明书进行骨髓dna提取。
送检样本前应用紫外分光光度计测定吸光度值(a)以确定其含量及纯度，并用琼脂糖凝胶电泳检测dna的完整性。本实验对dna样本的要求：dna总量≥3μg，浓度≥50ng/μl，纯度od260/280≥1.8。此外，在dna完整性方面，要求经琼脂糖凝胶电泳检测无拖尾、无弥散带。
4、全基因组cpg位点甲基化捕获测序
4.1文库构建
一、dna的片段化
预冷超声打断仪(4℃)，将5.8μl稀释过的硫化转换对照试剂加入至1μgdna中，与足量的eb混合至50μl，在超声作用下打断为200bpdna片段。具体操作方法见《bioruptor标准操作流程》。然后用琼脂糖电泳对dna片段大小进行检测是否合格。
二、末端修复
在ep管中配制修复体系，具体成分及用量见表3。混合后于20℃温浴约半小时，然后用磁珠对dna回收纯化，并加入42μleb将其溶解。
三、加“a”
在ep管中配制加“a”体系，具体成分及用量见表4。配好后混合，于37℃温浴约半小时，然后用磁珠对dna回收纯化，并加入20μleb溶解。
四、接头连接
用加完“a”的dna产物进一步配制连接体系，见表5。配好后混合，于
16℃温浴约半小时，然后用磁珠对dna回收纯化，并加入30μleb溶解。
表3修复体系
表4加“a”体系
表5连接体系
五、硫化处理
使用zymo公司ezdnamethylation-goldtmkit对连接产物进行亚硫酸盐处理，使得未甲基化的c转变为t。具体步骤参照ezdnamethylation-goldtmkit步骤进行操作，最终溶解于20ulddh2o中。
六、pcr扩增及产物纯化
将加完接头后的dna产物用pcr进行扩增并进一步纯化。
(1)pcr扩增体系见表6
表6pcr扩增体系
引物序列：上游：agtggttaattaattttcgagtttc
下游：tattattaccctataaccaccatcg
(2)pcr反应
设定pcr程序：
(3)pcr产物的回收纯化
用磁珠对pcr扩增后的dna产物进行回收纯化，溶于ddh2o中待下一步使用。
4.2目标区域捕获测序
一、杂交前文库与芯片温浴杂交
(1)将等量的待杂交文库进行混合，制成样本文库(1μg)。
(2)将每个样本文库对应的index等摩尔混合，制成index文库(1000pmol/l)。
(3)将10μl硫化捕获增强剂与步骤(1)中的混合物加入ep管中，并向其中加入10μl公共引物(1000pmol/l)和10μl步骤(2)中的混合物(1000pmol/l)。
(4)用密闭膜封闭步骤(3)中带有杂交样本的ep管，并在膜上打几个孔，置于真空干燥浓缩仪上(温度50℃)抽干，直至完全干燥。
(5)向ep管的干粉中加入以下两种成分：2×杂交缓冲液和杂交组分a(表7)，震荡后离心10秒，然后放于95℃温浴仪10分钟使其变性。
表7杂交样本的组分及含量
(6)到时间后快速离心，将步骤(5)中混合物加至含有nimblegen芯片的pcr管中。震荡3秒后离心，置于47℃加热模块上68-72h，设置热盖至57℃。
二、目的dna片段洗脱
(1)将编号为试管1、试管2、试管3、试管4的洗涤液原液配成1×缓冲液(表7)。
(2)将dynabeadsm280提前置于室温半小时，实验前震荡1分钟，将其充分混匀，分装50μl于离心管中，并放于磁力架上静置5分钟去掉上清，然后将两倍磁珠体积的1×bead洗涤缓冲液加入，震荡10秒后放于架上，将上清弃去。用50μl1×bead洗涤缓冲液重悬并转入新的已标记好的ep管中，再次放于磁力架上去上清。
表8磁珠洗涤液成分及含量
(3)向步骤(2)的ep管中加入杂交的dna文库，将加热模块设定47℃加热，并间隔震荡混匀45分钟。然后加入47℃50μl体积的1×洗涤缓冲液i，并混匀。将混匀产物转入新的已标记好的ep管中，放至磁力架去上清。再加入100μl47℃的1×stringent洗涤缓冲液，将其混匀后放于47℃加热模块上温浴5min。
(4)将1×洗涤缓冲液i100μl加入到步骤(3)的ep管中，混合后放至磁力架去上清。然后加入1×洗涤缓冲液ii100μl，混合后放至磁架去上清。再加入1×洗涤缓冲液iii100μl，混合后放至磁力架去上清。最后向ep管中加入50μlddh2o，洗脱beads并捕获dna使之溶解。将beads-dna混合物保存于-20℃。
三、目的dna片段的pcr扩增
(1)制备捕获后的pcr混合物(表9)
表9捕获后的pcr混合物
引物1序列：上游：gttaggtagggaagaagggagtagt
下游：cccaaaaatcaaataatcaaaaaaa
引物2序列：上游：gtggttaattaatttttgagttttgt
下游：tattaccctataaccaccatcacc
(2)pcr反应
设定pcr程序：
(3)pcr产物的回收纯化
用磁珠将扩增后的dna产物进行回收纯化，并溶于ddh2o中待测序。
四、hiseq2500测序
杂交捕获后的样本用hiseq2500pe100进行测序。
5、甲基化测序数据分析
数据下机后，首先进行数据过滤，去掉低质量数据，得到可用数据。数据检测合格后，将可用数据与参考基因组进行比对，得到比对结果。在确认比对质量合格后，使用唯一比对数据得到c碱基甲基化信息，进行信息分析处理，得到标准信息分析结果和个性化分析结果。
5.1数据过滤及比对
一、原始数据
illumina高通量测序结果最初以原始图像数据文件存在，经casava软件进行碱基识别后转化为原始测序序列，其结果以fastq(简称为fq)文件格式存储。fastq文件包含每条测序序列的名称、碱基序列以及其对应的测序质量信息。在fastq格式文件中，每个碱基对应一个碱基质量字符，每个碱基质量字符对应的ascii码值减去33(sanger质量值体系)，即为该碱基的测序质量得分。不同测序质量得分代表不同的碱基测序错误率，如测序质量得分值为20和30分别表示碱基测序错误率为1％和0.1％。
二、数据过滤
对测序得到的原始序列进行过滤，得到高质量的cleanreads，再进行后续分析，后续分析都基于cleanreads。
数据处理步骤如下：
(1)去除接头污染的reads(reads中接头污染的碱基数大于5bp。对于pe，若一端受到接头污染，则去掉两端的reads)；
(2)去除低质量的reads(reads中质量值q≤19的碱基占总碱基的15％以上，对于pe，若一端为低质量reads，则会去掉两端reads)；
(3)去除含n比例大于5％的reads(对于pe，若一端含n比例大于5％，则会去掉两端reads)。
数据过滤统计指标如下：
readslength(bp)：序列长度；
rawreads：原始下机的序列数；
rawbases：原始下机序列的碱基数；
cleanreads：过滤后得到的高质量的序列数；
cleanreadsrate(％)：过滤后得到的高质量序列占原始下机序列的比例。这个值越大，说明测序质量或者文库质量越好；
cleanbases：过滤后的高质量序列的碱基数；
low-qualityreads：由于低质量碱基过多，被去掉的序列数；
low-qualityreadsrate(％)：由于低质量碱基过多，被去掉的序列占原始下机序列的比例；nsreads：由于含n过高，被去掉的序列数；
nsreadsrate(％)：由于含n过高，被去掉的序列占原始下机序列的比例；
adapterpollutedreads：由于某些序列受到接头的污染，被去掉的序列数；
adapterpollutedreadsrate(％)：由于某些序列受到接头的污染，被去掉的序列占原始下机序列的比例；
originalq30basesrate(％)：rawreads中测序质量值大于30(错误率小于0.1％)的碱基占总碱基的比例；
cleanq30basesrate(％)：cleanreads中测序质量值大于30(错误率小于0.1％)的碱基占总碱基的比例。
三、测序数据质量值分布
测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。每个碱基测序错误率是通过测序phred数值(phredscore，qphred)通过公式转化得到，而phred数值是在碱基识别(basecalling)过程中通过一种预测碱基判别发生错误概率模型计算得到的，对应关系如下表10所示：
表10illuminacasava碱基识别与phred分值之间的简明对应关系
四、测序数据碱基分布
以过滤后序列的碱基位置作为横坐标，以每个位置的atcg碱基(其中n表示未知的碱基)的比例作为纵坐标，得到碱基分布图。
五、比对分析
在信息分析过程中，测序的结果和参考基因组都进行了c-to-t(forward)和g-to-a(reverse)的转换，将转换后的reads和转换后的基因组序列进行比对，唯一比对reads将用于甲基化信息的分析。用于比对的基因组数据库位universityofcaliforniasantacruzhg19。使用软件及参数[24](版本：v0.9.0)：bismark-p5-n1–directional-odirectory-1fq1-2fq2。
5.2甲基化位点识别及计算甲基化水平
一、甲基化位点识别
不同分布类型的甲基化c位点在不同物种基因组中出现比例不同[25]，因此，各类型mc(mcg、mchg和mchh)的位点数目及其在全部mc的位点中所占的比例(例：mchg所占比例＝mchg数目/mc的总数)，在一定程度上反映了特定物种的全基因组甲基化图谱的特征。
二、计算甲基化水平
每一个甲基化c碱基的甲基化水平均按如下公式进行计算：c位点的甲基化水平＝支持甲基化的reads的数目/(支持甲基化的reads的数目+支持非甲基化的reads的数目)×100％。捕获区域平均甲基化水平反应了基因组甲基化图谱的总体特征。这里只考虑覆盖深度不小于5的甲基化位点。
不同基因区域的平均甲基化水平计算公式如下：某区域平均甲基化水平＝该区域支持甲基化的reads的数目/该区域总的reads的数目×100％。
5.3差异甲基化区域分析
本研究在多个不同分组间、不同基因区域寻找差异甲基化区域(dmrs)，dmr需要符合以下三个条件：至少包含5个cpgs；其中至少有3个差异甲基化位点(差异甲基化位点即差异甲基化水平≥20％的位点)；区域差异甲基化水平≥20％。使用软件及参数(v0.9.2)：r包methylkit[26]和edmr[27]。
将找到的dmrs注释到genebody，比较dmrs在genebody不同区域分布情况，更好的理解甲基化修饰的改变对基因调控的影响。
5.4差异甲基化基因筛选及临床意义验证
本发明的目的是寻找aml患者中高甲基化预后不良的基因，因此以细胞遗传学和分子遗传学预后分组为基础，在不同预后分组间进行差异甲基化比较，寻找在危险度高的组中高甲基化差异区域，并根据差异区域注释差异基因，进一步进行基因功能分析，筛选确定目标基因，然后在研究病例和tcga数据库病例中进行临床意义验证，包括与临床因素之间的关系、缓解率及生存意义分析等。
6、统计学方法
连续变量采用‘学生t检验、方差分析、mann-whitneyu秩和检验；分类变量采用fisher’s精确性检验或卡方检验。生存分析采用kaplan-meier曲线。总生存率(overallsurvival，os)定义：从初诊到患者死亡或随访终止；无病生存率(diseasefreesurvival，dfs)定义为从完全缓解到复发、死亡、或者随访终止。双侧检验p<0.05为显著性统计差异。上述统计分析在spss19.0(ibmcorp.，armonk，ny，usa)、graphpadprism6(graphpadsoftwareinc.，sandiego,california，usa)等软件上进行。
三、实验结果
1、测试数据覆盖度结果
采用所有上机测序的44例次样本进行测序数据质控分析。我们捕获的目标覆盖全基因组240,513个区域、包含大于5mb的cpg位点，共有80m的碱基对。测序读长为125bp。覆盖度不低于1×,5×,10×,20×的cpg位点平均百分比分别为92.73％,80.32％,67.20％,和43.81％。所有44例次dna样本覆盖度超过30×的数据量有430gb。所有样本的转化率均超过99.5％。
2、测试数据过滤结果
过滤剔除质量差的数据后，44例次样本平均的总cleanreads数目为58,147,036(范围：40,663,694-79,302,058)，测序质量值cleanq30baserate平均为95.03％(范围：92.65％-98.05％)。
3、测试数据比对结果
与参考基因组比对，cleanreads平均比对到基因组的目标捕获cpgs为72.32％(范围：38.54％-85.31％)，总体比对率为92.90％(范围：80.52％-96.40％)。
4、技术重复样本结果
为了验证mcc-seq检测骨髓样本dna甲基化的稳定性，我们对两例aml患者的样本复本进行了独立、平行的检测。结果显示，不同测序深度下同一患者的两个复本dna甲基化测序结果具有高度的一致性，并且随测序深度的增加，相关系数逐渐增大，相关性更强(s22-重复1和s22-重复2比较：在测序深度的界值分别为1×，5×，10×，20×时，r＝0.959，0.971，0.978，0.985；s23-重复1和s23-重复2比较：在测序深度的界值分别为1×，5×，10×，20×时，r＝0.954，0.967，0.974，0.982)。
5、aml患者dna甲基化水平影响因素
我们计算了21例初治aml样本的全基因组范围内捕获的cpgs的甲基化水平(dnamethylationindicator，dmi)，并分析了各临床特征与dmi的关系，以确定是否有特定的因素影响患者的dmi。结果显示，dmi不受骨髓样本原始细胞比例、患者年龄、性别、fab分型、细胞遗传学危险度分组、分子遗传学危险度分组以及体细胞突变的数目的影响。同时发现，年龄≥50岁的患者与年龄<50岁的患者相比，dmi显著增高(49.39％±3.43％vs.46.75％±2.31％，p＝0.048)(图1)。
6、不同基因区域dna甲基化水平
我们分析了不同基因区域，如所有cpg位点捕获区域(allcpgs)，cpg岛(cpgislands)，启动子区(promoters)，外显子区(exons)，第一外显子(exon1)，内含子(introns)，增强子(enhancers)，5’非编码区(5’untranslatedregion,5’utr)等的dna甲基化水平。通过比较初治aml患者和正常对照骨髓这些区域的甲基化水平发现，仅启动子区和增强子区的dna甲基化水平(dmi)显著高于正常对照(p＝0.025和p＝0.021)(图2a)；进一步通过8对初治-缓解配对的样本分析发现，在诱导缓解期，仅启动子区的dna甲基化水平显著低于初治(p＝0.018)，而增强子区dna甲基化水平并无明显变化(p＝0.145)(图2b)；另外3对初治-复发配对样本的结果提示，复发时启动子区的dna甲基化水平与初治时无显著差异(p＝0.305)。上述结果显示，aml患者启动子区dna甲基化水平高于正常对照，且缓解期显著下降，复发时再度升高，说明启动子区dna甲基化是aml发病的一个显著特征，且与临床治疗反应密切相关。
为了寻找具有功能的差异甲基化基因，我们首先在初治aml样本和正常对照样本间进行了相关基因功能区域(如启动子、第一外显子、增强子和5’非编码区，即promoter，exon1，enhancer，5’utr)的差异甲基化区域(dmr)比较。结果发现，启动子区具有最多的差异甲基化区域(60.9％，669/1099，p<0.001))，且以高甲基化差异区域为主(5.0％，502/669，p<0.001)(图2c)；从整合遗传学预后和表观遗传学预后的角度出发，我们进一步比较了初治aml不同细胞遗传学预后和分子遗传学预后分组间的差异甲基化区域(dmr)，并根据不同基因区域的dmr注释相应的差异甲基化基因(differentiallymethylatedgenes，dmgs)，结果显示，注释到启动子区的dmgs最多(p<0.001)(图2d)。上述结果提示，启动子区具有最重要的功能性差异甲基化区域，并且与已知的遗传学预后分组具有联系。基于启动子区的差异甲基化进行甲基化预后相关基因的筛选具有实用性。
7、差异甲基化基因筛选结果
为了筛选启动子区高甲基化而导致患者预后不良的基因，确定一组整合了遗传学和表观遗传学的预后基因标志物，以便研发用于判断aml预后的试剂盒，我们进行了筛选和验证。首先，按照nccn指南预后危险度分层标准，将进行了mcc-seq测序的21例初治aml分成不同的细胞遗传学和分子遗传学预后分组，并比较不同预后分组间启动子区的甲基化差异，比较的组别方向是：高危vs.中危，中危vs.低危，高危vs.低危。细胞遗传学组间比较共找到100个高差异甲基化基因(hyper-dmgs)，分子遗传学组间比较共找到44个高差异甲基化基因。然后进一步采用“取交集”的方法，即选择共同出现在细胞遗传学组间比较的100个hyper-dmgs中和分子遗传学组间比较的44个hyper-dmgs中的hyper-dmgs，从而得到了在细胞遗传学组间比较和分子遗传学组间比较中预后均不好的18个hyper-dmgs，我们进一步通过检索基因数据库(ncbigenedatabase)、检索相关基因的文献报道、对基因进行go功能分析和kegg通路分析，明确筛选到的18个基因的特征，并最终确定其中的12个基因用于后续临床意义验证。18个基因中，剔除3个假基因(gucy1b2，hnrnpa1p33和tuba3fp)和2个无分子功能的基因(mir3150b和mir4638)，及1个在tcga数据库中无dna甲基化检测结果的基因(plec)。最终有6个基因被剔除，确定了一组12个与遗传学预后相关、受dna甲基化调控的功能性差异甲基化基因(bard1，bcl9l，clec11a，defb1，foxd2，igf1，il18，itih1，lsp1，p2rx6，rnase3和tubgcp2)，并在研究病例和tcga病例中进行临床意义验证，证明了其具有判断aml预后的意义，可用于后续预后分层试剂盒的研发。
8、差异甲基化基因临床意义验证结果
根据筛选到的12个dmgs，我们计算了每一例患者每一个dmgs的dna甲基化水平。同一患者12个dmgs甲基化水平的平均值记为该患者的综合甲基化值(m-value)，用于代表患者的甲基化水平。然后，我们在21例mcc-seq的研究病例中和169例tcga病例中验证了m-value的临床意义。
结果发现，m-value不受患者年龄、性别和样本原始细胞比例的影响，这与捕获范围内的cpgs的甲基化水平不受这些因素影响的结果类似(图1)。
我们评价了m-value与临床短期疗效，即诱导缓解率的关系。结果显示，完全缓解患者的平均m-value值低于未缓解者(37.42％±15.79％vs.49.69％±12.51％，p＝0.09)。为了更进一步评价m-value与cr的关系，我们以21例患者m-value值的中位数为界，将患者分为低m-value组，即m-value值≤中位数(n＝11)；和高m-value组，即m-value值>中位数(n＝10)。在低m-value组，有高达90.9％(10/11)的患者达到了cr，缓解率显著高于高m-value组(40.0％,4/10；p＝0.024)。并且，6名预后中危(intermediate-riskaml，ir-aml)且低m-value的患者，有5例达到了cr，缓解率为83.3％，而8例ir-aml且高m-value的患者，仅3例达到了cr，缓解率仅为37.5％。上述结果提示，无论对于所有非m3aml还是ir-aml，低m-value的患者更容易获得cr。
评价m-value与远期疗效，即生存的关系：我们将研究病例及tcga数据库病例均按其各自的m-value的中位数分为低m-value组和高m-value组。首先，在21例研究病例中，总的中位总生存期(os)、中位无病生存期(dfs)、1年累积os、1年累积dfs分别为23.8个月、未达到、78.9％和69.1％。在低m-value组和高m-value组，中位os分别为未达到和14.93个月(p＝0.062)；中位dfs分别为未达到和10.97个月(p＝0.039)；1年累积os分别为88.9％和68.6％(p＝0.145)；1年累积dfs分别为90.9％和30.0％(p<0.001)。生存分析显示，高m-value是dfs的危险因素(hr：6.83，95％ci：1.07-40.28；p＝0.039)；对于os，高m-value有预后更差的趋势。在tcga病例中，高m-value组患者的os显著差于低m-value组(medianos：15.1monthsvs.16.4months；hr：1.491，95％ci：1.043-2.151，p＝0.038)。对于dfs，高m-value有不良预后的趋势。并且，高m-value患者的2年累积os和dfs均显著低于低m-value组(os：35.9％vs.45.9％，p＝0.001；dfs：30.2％vs.44.2％，p<0.001)。上述结果表明，高m-value是预后不良因素。
参考文献
1.chen,w.,etal.,cancerstatisticsinchina,2015.cacancerjclin,2016.
2.hasserjian,r.p.,acutemyeloidleukemia:advancesindiagnosisandclassification.intjlabhematol,2013.35(3):p.358-66.
3.walter,r.b.,etal.,significanceoffabsubclassificationof"acutemyeloidleukemia,nos"inthe2008whoclassification:analysisof5848newlydiagnosedpatients.blood,2013.121(13):p.2424-31.
4.prada-arismendy,j.,j.c.arroyave,ands.rothlisberger,molecularbiomarkersinacutemyeloidleukemia.bloodrev,2017.31(1):p.63-76.
5.papaemmanuil,e.,etal.,genomicclassificationandprognosisinacutemyeloidleukemia.nengljmed,2016.374(23):p.2209-21.
6.taylor,j.,w.xiao,ando.abdel-wahab,diagnosisandclassificationofhematologicmalignanciesonthebasisofgenetics.blood,2017.130(4):p.410-423.
7.dohner,h.,etal.,diagnosisandmanagementofacutemyeloidleukemiainadults:recommendationsfromaninternationalexpertpanel,onbehalfoftheeuropeanleukemianet.blood,2010.115(3):p.453-74.
8.bennett,j.m.,etal.,proposalsfortheclassificationoftheacuteleukaemias.french-american-british(fab)co-operativegroup.brjhaematol,1976.33(4):p.451-8.
9.daniela.arber,etal.,the2016revisiontotheworldhealthorganizationclassificationofmyeloidneoplasmsandacuteleukemia.blood,2016.127(20):p.2391-405.
10.bullinger,l.,etal.,identificationofacquiredcopynumberalterationsanduniparentaldisomiesincytogeneticallynormalacutemyeloidleukemiausinghigh-resolutionsingle-nucleotidepolymorphismanalysis.leukemia,2010.24(2):p.438-49.
11.mardis,e.r.,etal.,recurringmutationsfoundbysequencinganacutemyeloidleukemiagenome.nengljmed,2009.361(11):p.1058-66.
12.shen,y.,etal.,genemutationpatternsandtheirprognosticimpactinacohortof1185patientswithacutemyeloidleukemia.blood,2011.118(20):p.5593-603.
13.gutierrez,s.e.andf.a.romero-oliva,epigeneticchanges:acommonthemeinacutemyelogenousleukemogenesis.jhematoloncol,2013.6:p.57.
14.figueroa,m.e.,etal.,leukemicidh1andidh2mutationsresultinahypermethylationphenotype,disrupttet2function,andimpairhematopoieticdifferentiation.cancercell,2010.18(6):p.553-67.
15.schoofs,t.,w.e.berdel,andc.muller-tidow,originsofaberrantdnamethylationinacutemyeloidleukemia.leukemia,2014.28(1):p.1-14.
16.figueroa,m.e.,etal.,dnamethylationsignaturesidentifybiologicallydistinctsubtypesinacutemyeloidleukemia.cancercell,2010.17(1):p.13-27.
17.deneberg,s.,etal.,gene-specificandglobalmethylationpatternspredictoutcomeinpatientswithacutemyeloidleukemia.leukemia,2010.24(5):p.932-41.
18.deneberg,s.,etal.,prognosticdnamethylationpatternsincytogeneticallynormalacutemyeloidleukemiaarepredefinedbystemcellchromatinmarks.blood,2011.118(20):p.5573-82.
19.sun,z.,etal.,baseresolutionmethylomeprofiling:considerationsinplatformselection,datapreprocessingandanalysis.epigenomics,2015.7(5):p.813-28.
20.allum,f.,etal.,characterizationoffunctionalmethylomesbynext-generationcapturesequencingidentifiesnoveldisease-associatedvariants.natcommun,2015.6:p.7211.
21.nimblegen.www.nimblegen.com.
22.walker,d.l.,etal.,dnamethylationprofiling:comparisonofgenome-widesequencingmethodsandtheinfiniumhumanmethylation450beadchip.epigenomics,2015.7(8):p.1287-302.
23.cancergenomeatlasresearch,n.,genomicandepigenomiclandscapesofadultdenovoacutemyeloidleukemia.nengljmed,2013.368(22):p.2059-74.
24.krueger,f.ands.r.andrews,bismark:aflexiblealignerandmethylationcallerforbisulfite-seqapplications.bioinformatics,2011.27(11):p.1571-2.
25.lister,r.,etal.,humandnamethylomesatbaseresolutionshowwidespreadepigenomicdifferences.nature,2009.462(7271):p.315-22.
26.akalin,a.,etal.,methylkit:acomprehensiverpackagefortheanalysisofgenome-widednamethylationprofiles.genomebiol,2012.13(10):p.r87.
27.li,s.,etal.,anoptimizedalgorithmfordetectingandannotatingregionaldifferentialmethylation.bmcbioinformatics,2013.14suppl5:p.s10.