CN115651973B - 一种可传代细胞的高保真甲基化位点的分离与分析方法 - Google Patents

Abstract

本申请提供一种可传代细胞的高保真甲基化位点的分离与分析方法，包括步骤：(A)将可传代细胞进行传代培养多代后，取不同代数的细胞进行基因组DNA提取，作为实验样本；(B)将spike‑in加入实验样本中，进行重亚硫酸盐转化，构建简并代表性重亚硫酸盐测序文库，并进行测序；所述spike‑in为甲基化和未甲基化的lambda DNA；(C)对于测序数据，根据spike‑in计算转化效率和错误转化率最终对甲基化水平的评估产生偏差，并且对每个位点的DNA甲基化水平采用二项分布近似贝叶斯模型进行定量检测，通过比较不同代别的甲基化水平得到高保真甲基化位点。

Description

一种可传代细胞的高保真甲基化位点的分离与分析方法

技术领域

本发明涉及分子生物学和生物信息学技术领域，更具体地，涉及一种可传代细胞的高保真甲基化位点的分离与分析方法。

背景技术

哺乳动物DNA甲基化作为一种重要的表观遗传学修饰，改变遗传表现，参与哺乳生物的基因表达、胚胎发育、疾病发生等重要生物学过程。DNA甲基化通常发生在CpG位点的胞嘧啶上，在DNA甲基化转移酶的作用下，在基因组CpG二核苷酸上结合一个甲基基团。DNA甲基化是细菌识别区分内外源DNA的主要方式，是细菌识别自身的身份标记。大量研究表明DNA甲基化水平异常与肿瘤发生、恶化和细胞癌变进程有着密切的联系。

甲基化错误的积累与多种疾病的发作相关，包括癌症和各种发育、自身免疫和神经系统疾病。正常的DNA甲基化维持对于细胞组织的发育和行使功能至关重要，找到在细胞传代过程中甲基化水平稳定不变的CpG位点，即高保真的CpG甲基化位点，对研究甲基化维持的机理以及生物学意义都至关重要。并且，我们可以通过设计探针或构建序列敲除的细胞株，研究这些高度保守的甲基化位点的特性与可能发挥的重要生物学功能，对于理解生物的发育和衰老、疾病的发生和发展都具有重要的价值。

通常认为，DNA甲基化在细胞分裂过程中能够稳定的保持下去，细胞通过不同的甲基化转移酶的协同作用来实现在DNA复制过程中DNA甲基化位点的保持，从而确保其功能在子代细胞中的延续。现有对于甲基化保真度的研究大多是基于多细胞的平均甲基化，即对于多细胞组织的甲基化检测结果来说，其某个位点的甲基化水平变化其实是该组织中某些细胞个体在该CpG 位点发生了甲基化的丢失或获得，才使该位点的平均甲基化水平发生改变，产生了部分甲基化。但由于目前的研究基本是对于单一时间点的多细胞进行甲基化检测，这样就难以找到在细胞传代过程中甲基化水平保持稳定不变的 CpG位点，而后者由于其高度的保守型往往对应了更为重要的功能。因此，提供一种基于可传代细胞的分离高保真甲基化位点的方法具有重要意义。

发明内容

本发明的目的在于，提供一种可传代细胞的高保真甲基化位点的分离与分析方法，能够分离出在细胞传代过程中，甲基化水平稳定不变的CpG位点。

一种可传代细胞的高保真甲基化位点的分离与分析方法，包括步骤：

(A)将可传代细胞进行传代培养多代后，取不同代数的细胞进行基因组DNA提取，作为实验样本；

(B)将spike-in加入实验样本中，进行重亚硫酸盐转化，构建简并代表性重亚硫酸盐测序(RRBS)文库，并进行测序；所述spike-in为甲基化和未甲基化的lambda DNA(λ噬菌体中的DNA)；

(C)对于测序数据，根据spike-in计算转化效率和错误转化率最终对甲基化水平的评估产生偏差，并且对每个位点的DNA甲基化水平采用二项分布近似贝叶斯模型进行定量检测，通过比较不同代别的甲基化水平得到高保真甲基化位点。

在一些实施方式中，在步骤(A)中，所述可传代细胞为哺乳动物细胞或细胞系，来自于生物组织、含有细胞的体液；所述传代培养的细胞必须是同源且表型单一的。

进一步的，将传代培养10代-40代的细胞作为起始样本，从中选择2-4 个培养了不同代数的细胞样本，将各细胞样本同步化到G1/G0期后进行细胞收集和基因组DNA提取，作为实验样本。

进一步的，所述细胞样本的细胞数量不少于10⁵，实验样本的基因组DNA 的含量不少于4μg，每次取不少于2μg用于建库。

在一些实施方式中，在步骤(B)中，将实验样本和spike-in分别进行 DNA片段化处理；然后分别进行末端修复及加A、连接接头、USER酶消化及片段筛选纯化；再将甲基化和未甲基化的lambda DNA混合后按照0.5％～1％的比例加入到实验样本基因组DNA中进行重亚硫酸盐转化，扩增纯化后得到重亚硫酸盐测序文库，上机进行高通量测序。

进一步的，所述DNA片段化处理采用超声波打断技术，所述连接接头使用的接头为甲基化专用的接头。

在一些实施方式中，在步骤(C)中，利用lambda DNA计算得到实验样本的每个CpG位点的重亚硫酸盐平均转化效率CR和平均错误转化率ICR，再根据平均转化效率CR和平均错误转化率ICR计算实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平ML和成功率θ，根据成功率θ计算观测值Y和先验值X，将先验值X与观测值Y进行比较，得到后验分布，分类得到高保真甲基化位点。

进一步的，利用未甲基化的lambda DNA计算实验样本的每个CpG位点的重亚硫酸盐平均转化效率CR，计算公式为S1，利用甲基化的lambda DNA 计算实验样本的每个CpG位点的重亚硫酸盐平均错误转化率ICR，计算公式为S2，spike-in的重亚硫酸盐转化效率CR和重亚硫酸盐错误转化率ICR与实验样本相同，计算公式S1和S2分别为：

其中，U_c为未甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，U_mc为未甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目， M_c为甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，M_mc为甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目。

进一步的，利用平均转化效率CR和平均错误转化率ICR计算实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平ML，spike-in的重亚硫酸盐甲基化水平与实验样本相同，重亚硫酸盐转化后的甲基化水平ML的计算公式为S3：

ML＝(CR-ICR)×ml+1-CR (S3)

其中，ML为实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平，ml为实验样本的每个CpG位点的重亚硫酸盐转化前的甲基化水平。

进一步的，由于传代培养的细胞均来源于同一个细胞，实验样本在重亚硫酸盐转化前高保真位点的甲基化水平ml理论上应为100％、50％或0％，引入平均转化效率CR和平均错误转化率ICR的参数，根据公式S3得到的重亚硫酸盐转化后的甲基化水平ML即为成功率θ。

进一步的，根据成功率θ计算观测值Y，利用二项分布似然函数S4进行模拟，设每个CpG位点的测序深度为n，即以成功率θ进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为观测值Y，计算观测值Y的公式为S4，p(x＝k)是事件成功的次数等于观测值Y时的概率，公式S4为：

其中，k为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，θ为事件成功的概率即重亚硫酸盐转化后的甲基化水平ML，n为独立实验的次数即该CpG位点测序深度。

进一步的，由于每个CpG位点重亚硫酸盐转化后的甲基化水平ML为先验分布,先验概率为π(θ)，根据先验概率随机采样一个参数θi，设每个CpG 位点的测序深度为n，即以参数θi进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为先验值X，计算先验值X的公式为S5，p(x|ML)是事件成功的次数等于先验值X时的概率，公式S5为：

其中，x为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，重亚硫酸盐转化后的甲基化水平ML为事件成功的概率，n为独立实验的次数即该CpG位点测序深度。

进一步的，将观测值Y与先验值X进行比较，若先验值X与观测值Y 相等，则保留随机采样的参数θi作为后验分布的样本；若先验值X与观测值 Y不相等，则舍弃，返回上一步重新选择一个参数θi，直到有N＝1000个样本被保留。

进一步的，我们得到后验分布即所有甲基化保真位点的水平分布，用R 语言的quantile函数取95％的置信下限作为阈值，若CpG位点的甲基化水平大于该阈值则为保真位点，反之为动态位点，该保真位点即为高保真的甲基化CpG位点。

与现有技术相比，本发明的优势为：

1.将可传代细胞进行传代培养多代，传代培养多代的细胞作为起始样本，动态位点甲基化水平发生了较大的变化，在统计学上检测到动态位点甲基化水平变化的统计功效较强，从而分离出高保真甲基化位点和动态位点。

2.为了避免不同细胞表型中甲基化的异质性带来的影响，我们设计了所提供进行传代培养的细胞需要来源于同一个细胞表型。

3.不同于大多数检测CpG甲基化的研究只测量了一个时间点的甲基化水平，而没有比较不同代之间的水平，本发明选择不同的两个细胞传代时期,通过比较两个时间点从而分离出随着细胞传代严格保持甲基化信息的高保真位点。

4.建立了定量分类甲基化高保真位点的数据分析方法，对于每个位点的 DNA甲基化水平的检测为二项分布的采样过程，符合二项分布，根据lambda DNA的CpG位点转化参数得到重亚硫酸盐转化后的甲基化水平呈现一个先验分布，近似贝叶斯的计算方法植根于贝叶斯统计可以被用来估计模型参数的后验分布，因此采用二项分布近似贝叶斯的算法，从而在模型上考虑到了计算所需的参数不是明确和固定的特殊特点，使得检测结果更准确。

附图说明

结合以下附图一起阅读时，将会更加充分地描述本申请内容的上述和其他特征。可以理解，这些附图仅描绘了本申请内容的若干实施方式，因此不应认为是对本申请内容范围的限定。通过采用附图，本申请内容将会得到更加明确和详细地说明。

图1为两个实验样本(Jurkat 20代和30代)文库纯化后得到的文库片段分布图。

图2为两个实验样本(Jurkat 20代和30代)的CpG位点甲基化水平散点图。

具体实施方式

描述以下实施例以辅助对本申请的理解，实施例不是也不应当以任何方式解释为限制本申请的保护范围。

下列实施例中未注明具体条件的实验方法，按照常规实验条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数按重量计算。除非有特别说明，否则实施例所用的材料均为市售产品。

实施例1：实验细胞样本的传代培养

本申请的可传代细胞为哺乳动物细胞或细胞系，来自生物样品，例如生物组织、含有细胞的体液。该可传代细胞可以通过已知的方式获取，例如通过血液、血浆、口腔粘膜样本、鼻腔样本、胚胎组织、内皮细胞等获取，生物样品可以是任何适当的形式提供，例如可以是液氮冷冻的形式、石蜡包埋的形式、新鲜分离的形式等。

本实施例采用分离的人体细胞系，具体的，为人急性白血病T淋巴细胞株Jurkatclone E6-1，然后进行细胞传代培养，步骤为：

1.将人急性白血病T淋巴细胞株Jurkat clone E6-1稳定传代后使用5％二甲亚砜冻存液保存于液氮中，然后进行复苏。

2.将冻存的细胞株置于37℃无菌水浴锅中快速解冻，随后将体系转移至 1.5mL离心管中，800rpm离心3分钟，小心吸弃上清。

3.将步骤2中的细胞沉淀用1mL 1×PBS重悬，用移液枪将细胞轻柔吹散，使细胞充分润洗，离心后吸除PBS，以去除瓶中的细胞碎片以及活力较差的细胞。

4.加入2mL新鲜培养液(RPMI1640，含有10％FBS和1％双抗，提前至恢复室温)至6cm皿中，将步骤3离心后的细胞用1mL新鲜培养液重悬，并转移到含有新鲜培养液的培养瓶中，轻柔摇匀，镜下观察细胞密度，确保密度合适且细胞较为分散即可。将细胞悬液置于37℃5％CO2的培养箱中培养，定期观察生长状况。

5.当细胞培养至第20代时，取出约四分之三的细胞，将其同步化到G1/G0 期后进行样品收集和DNA提取，剩余部分的细胞继续培养，在培养至30代时进行同步化，随后完成样品收集和DNA提取。

6.支原体检测：本研究在细胞培养过程中需要定期进行支原体检测，一旦发生污染，要及时丢弃样品并对培养环境进行支原体的全面检测和查杀，防止污染扩大。具体检测操作如下：

(1)在传代前收取培养时间约48h的细胞培养液，冻融后，作为模板 DNA，进行支原体的PCR检测。体系如下：

(2)同时设置阳性对照(支原体基因组作为模板)和阴性

对照(ddHO作为模板)，PCR程序如下：

步骤	温度(℃)	时间	循环次数
				1	95	10min	1
2	95	30s	35
				3	55	30s	35
4	72	1min	35
				5	72	5min	1
6	4	hold	1

(3)PCR结束后进行琼脂糖凝胶电泳：1％Agrose Gel，0.5％TBE缓冲液，120V30min。

实施例2：实验细胞样本同步化到G1/G0期

前期研究指出，DNA甲基化转移酶的修饰作用需要一定的时间，在细胞分裂的S期，新合成的DNA链其甲基化修饰具有滞后性，导致S期其基因组DNA甲基化水平显著低于其他时期(G1/G0和G2/M时期)。因此，为避免部分细胞在S期引起甲基化异质性，我们将上述细胞培养至第20代和第 30代，然后进行血清饥饿24h，使细胞同步到G1/G0期，然后收集细胞。同步化的具体步骤为：

1.用精密天平称取52mg的洛伐他汀化合物，溶于1.04mL 95％乙醇中，加入813uL的1mol/L的NaOH，充分混匀溶解，即可将洛伐他汀由失活的内酯前药转化为活性的二羟基开放酸。所得溶液用1mol/L HCl中和至pH 7.2，然后用蒸馏水或其他介质调至13mL的体积，并在-20℃的多等分中储存直至使用，储存液浓度为10mM。

2.将悬浮培养的实验细胞悬液转移到15mL离心管中，800rpm离心3分钟，弃上清，并用1×PBS清洗一次，然后加入5mL新鲜的RPMI 1640培养基(含有10％FBS和1％青链霉素)，并在其中加入15μL 10mM的洛伐他汀溶液，终浓度为30μM，充分混匀，置于37℃5％CO2的培养箱中培养24h。

3.细胞收集：24h后，将实验细胞悬液离心，弃去上清液，并用1×PBS清洗两次，随后取不少于2×105个细胞用500μL 1×PBS重悬，用于流式检测。

4.细胞固定：配制5mL 80％乙醇于15mL离心管中，混合均匀后打开盖子，将其垂直置于涡旋混匀仪上，低速涡旋，将步骤3中重悬的细胞用移液枪逐滴加入到乙醇中，使其充分分散均匀。将离心管封口，置于-20℃冰箱固定12h以上(最多可以保存两周)，待进行流式检测。

实施例3：实验细胞样本的基因组DNA提取

细胞传代培养至20代和30代的实验细胞样本，分别进行同步化到G1/G0 期后，分别进行基因组DNA提取，作为实验样本，基因组DNA提取的步骤具体为：

1.将-20℃保存的细胞沉淀取出，室温解冻后，使用350μL 1×PBS缓冲液将细胞悬浮起来，移液枪轻柔吹打。

2.向EP管中加入0.8μL RNase A，中速涡旋15s后室温下静置1min。

3.向EP管中加入150μL Buffer CL和8μL蛋白酶K(试剂盒提供，注意不要将蛋白酶K与Buffer CL溶液预混；第一次使用时用试剂盒提供的Buffer PK将蛋白酶K溶解后再使用)。涡旋1min，瞬时离心后，56℃水浴孵育10 min。

4.取出后，直接向EP管中加入350μL Buffer PD，于涡旋混匀仪上中速涡旋30s，室温离心，12000×g 10min。

5.将步骤4中的上清液小心吸取至制备管中(注意不要吸到沉淀)，放入收集管内，室温离心，12000×g 1min。

6.弃滤液，制备管内加入500μL Buffer W1，室温离心，12000×g 1min。

7.弃滤液，制备管内加入700μl Buffer W2(首次使用时，需按瓶身指示，加入相应体积无水乙醇至母液中)，室温离心，12000×g 1min。

8.重复上步(7)操作，再洗涤一次。

9.弃滤液，室温离心，12000×g 1min。

10.将制备管放入洁净的1.5mL离心管中，取30μL Eluent Buffer(提前预热至65℃提高洗脱效率)加入制备管膜中央，室温下静置1min。

11.室温离心，12000×g 1min，洗脱基因组DNA。重复用30μL Eluent Buffer再洗脱一次，提高产物得率。

12.使用NanoDrop One进行定量和质检：首先，用2μL双蒸水清洗点样孔，轻轻放下探头，接触5秒后抬起，用低尘纸轻轻吸去双蒸水。点击屏幕，选择核酸→DNA，吸取1μLEluent Buffer加入基座点样孔处，作为空白检测溶液，降下检测臂，点击空白检测后等待完成。随后抬起检测臂，用低尘精密擦拭纸擦干点样孔，取1μL待测DNA样品加入点样孔处，降下检测臂等待，检测完成后，记录好数据。将检测臂抬起，低尘纸擦净点样孔和检测臂探头，轻轻放回，结束实验。基因组DNA通常吸光度比值A260/280在1.8-2.0 之间。

实施例4：Spike-in的合成与处理

λDNA是来源于噬菌体的DNA，基因组片段总长为48502bp。我们选择其14556～15893bp之间的序列作为甲基化spike in片段的扩增模板, 19681～21559bp作为非甲基化spike in片段的扩增模板。根据选定区间的序列设计PCR引物(用于扩增甲基化区间的引物合成时已进行甲基化修饰)，使用dNTP Mix来扩增未甲基化的片段，使用5-MethylcytosinedNTP Mix来扩增甲基化的片段。跑胶验证后切胶回收，然后进行打断，具体步骤如下：

1.片段扩增。体系如下：

将体系分别混合均匀后，瞬时离心，置于PCR仪中运行以下程序：

步骤	温度(℃)	时间	循环次数
				1	95	30s	1
2	95	30s	20
				3	61	30s	20
4	68	30s	20
				5	48	5min	1
6	4	hold	1

结束后，进行2％琼脂糖凝胶电泳，以120V运行50min，紫外灯下切取扩增出的条带片段，转移到1.5mL离心管中。本文实验使用lambda DNA的 14556～15893bp序列扩增甲基化片段,19681～21559bp扩增非甲基化片段。

2.使用Axygen Gel Recovery Kit试剂盒，按照试剂盒的说明书步骤，对扩增片段进行胶回收。

3.使用Covaris打断仪将第2步得到的DNA进行片段化：

(1)将第2步得到的DNA加入到Covaris专用的超声破碎管内，补齐到 50μL(浓度不超过20ng/μL)。打开电脑和仪器开关，打开舱门，向样品槽内加双蒸水至水线，关闭舱门。

(2)打开软件，选择所需要的方案：

DNA_0300_bp_050μL_screw_cap_microTubeDNA_0300_bp_050μL_scre w_cap_microTube。

(3)待主页面的状态Water Tem，Water level，Door均为绿色时，移开金属头，将样品管放进样品槽内，放回探头。点击软件“Run”。

(4)运行结束后，将片段化的spike-in转移到新管中，1.5％凝胶电泳检测超声效率，确定已完成片段化后-20℃保存。

实施例5：构建简并代表性重亚硫酸盐测序(RRBS)文库

1.实验样本(20代和30代)提取的基因组DNA，分别进行MspI酶切。在无菌的200μLPCR管中配制如下酶切体系：

涡旋混匀仪上混合均匀，瞬时离心，置于PCR仪中，37℃作用2小时。酶切结束后取2μL酶切体系进行1.5％琼脂糖凝胶电泳，检测其酶切效率。当 DNA片段呈弥散状且分布均匀，说明酶切效率较好。

2.酶切后，分别使用AxyPrep PCR Clearn kit试剂盒，根据试剂盒说明书的步骤，回收酶切产物。

3.使用Ultra^TMDNA Library Prep Kit for Illumina将两个实验样本酶切后的DNA片段和spike-in分别进行文库构建。具体步骤为：

(1)末端修复及加A。将DNA 片段用Tris-HCl补齐至50μL，在200μL PCR管中配制一下体系：

将配好的体系置于涡旋混匀仪上混匀10秒，瞬时离心，置于PCR仪中，热盖温度设置为≥75℃，反应结束后将体系取出置于冰上或-20℃保存，反应条件为：

步骤	温度(℃)	时间
			1	20	30min
2	65	30min
			3	4	hold

(2)接头连接。将末端修复完成的体系取出，加入接头和连接酶，注意此处使用的接头为甲基化专用的接头。将配好的体系置于涡旋混匀仪上混匀 10秒，瞬时离心，置于PCR仪中，20℃反应15分钟，热盖关闭。体系如下：

(3)USER酶消化。反应结束后，将PCR管取出，在超净台内加入3μL USER酶，涡旋仪上混匀后离心，置于PCR仪中，37℃反应15分钟，热盖温度设置为≥47℃。

(4)片段筛选。将4℃保存的Beckman Ampure XP Beads提前取出，恢复室温至少30分钟。将建库体系用ddH2O补足到100μL，加入60μL Ampure XP Beads(0.6×)，充分涡旋混匀，室温静置5分钟，瞬时离心后置于磁珠吸附磁力架上，富集5分钟，待溶液澄清后，将上清转移至新管。

(5)向上清液新管中加入25μL Ampure XP Beads(0.25×)，充分涡旋混匀，室温静置5分钟，瞬时离心后置于磁力架上富集5分钟，小心吸弃上清，加入新配制的80％乙醇润洗磁珠，30秒后吸弃乙醇。重复操作一次。充分吸净乙醇后，晾干磁珠至表面无水膜，加入21μL10mM Tris-HCl(pH 8.0)，涡旋混匀，使磁珠充分分散在溶剂中，室温静置5分钟，瞬时离心后，置于磁力架上富集5分钟，小心吸取20μL转移至新管。

(6)建库的甲基化和未甲基化的spike-in片段由于已按照200bp进行了 Covaris打断，因此可直接用1.8×磁珠进行纯化，无需进行片段筛选，定量后将已建库的甲基化和未甲基化的spike-in按照1:1混合。

(7)Qubit定量。取198μL Qubit Buffer于配套的EP管中，加入1μL荧光染料，混合均匀，吸取1μL实验样本文库DNA加入到混合体系中，充分混合均匀，避光孵育2min后，进行Qubit定量。选择dsDNA high sensitivity模式，体积选择1μL，点击运行，记录结果(Qubit荧光分光光度计最低可检测到0.1ng/μL的浓度)。

(8)根据测定的结果，加入1％的已建库的spike-in(甲基化和未甲基化混合后)到上一步已完成定量的实验样本DNA文库中，混合均匀。

4.将上述混合均匀的混合物，使用Zymo Research EZ Gold Methylation Kit，按照试剂盒的说明书步骤，对上述混合物进行BS转化。

5.纯化

对上述BS转化体系进行纯化，具体步骤如下：

(1)加入600μL M-Binding Buffer到吸附柱中，将150μL的转化体系加入管中，吹打混匀。12000×g离心30s，弃液。

(2)加入100μL M-wash Buffer,12000×g离心30s，弃液。

(3)加入200μL M-Desulphonation Buffer,室温(20-30℃)静置 15-20min，12000×g离心30s，弃液。

(4)再加入200μL M-wash Buffer,12000×g离心30s，弃液。

(5)将吸附柱放入新的1.5mL EP管中，加入10μL M-Elution Buffer(加热到65℃可提高洗脱效率)，静置1min，12000×g离心30s。为保证DNA体系的纯净度，在试剂盒回收洗脱后，选择5×Ampure XP Beads再次纯化一次，操作同Ampure XP Beads纯化的一般操作。最终用20μL 10mM Tris-HCl洗脱。

(6)再加入200μL M-wash Buffer,12000×g离心30s，弃液。

(7)将吸附柱放入新的1.5mL EP管中，加入10μL M-Elution Buffer(加热到65℃可提高洗脱效率)，静置1min，12000×g离心30s。

(8)为保证DNA体系的纯净度，在试剂盒回收洗脱后，选择5×Ampure XP Beads再次纯化一次，操作同Ampure XP Beads纯化的一般操作。最终用 20μL 10mM Tris-HCl洗脱。

6.扩增与片段筛选

对上述纯化后的DNA进行扩增与片段筛选，步骤如下：

(1)取200μL离心管，在冰上配制PCR扩增体系：

混匀，瞬离，PCR程序为：

(2)反应结束后，用1.8×Ampure XP Beads对上述体系进行纯化，溶解于20μLTris-HCl中。

(3)为了更好的区分条带，获得更精确的片段分布结果，我们选择进行 6％TBE-PAGE垂直凝胶电泳。将Novex预制胶取出后，撕开胶条，开孔朝外，组装到Novex配套的垂直电泳槽中，倒入1×TBE电泳液，确保装置密封性较好，无缓冲液漏出。轻轻拔掉胶上部的梳子，将已经混合好Loading Buffer 的样品小心加入到电泳孔内，选择50bp DNA Ladder作电泳指示条带。打开电泳仪开关，设置恒压模式，120V电压电泳90min。

(4)电泳结束后，将胶板取出，撬开胶板外侧支撑板，轻柔取出胶块，放到装有50mL1×TBE溶液(含1/10000GelRed染料)的染色槽中，置于摇床上低速震荡5分钟，避光，随后使用凝胶成像仪，观察胶块上片段的分布。

(5)片段筛选与回收

a.自制离心套管：取1mL离心管，酒精灯加热1mL注射器针头，快速在离心管底扎4-5个孔，将扎孔的离心管置于1.5mL离心管中配套使用。

b.在紫外切胶台下，小心切取270-350bp范围的胶块(插入片段约 150-220bp)，转移到离心套管中。12000rpm离心5min，弃去1mL离心管，在1.5mL管中加入350μL 10mM Tris-HCl，充分震荡，并置于金属震荡仪上16℃孵育4小时。结束后，15000×g离心5min，吸取300μL上清到新管中。

在管内加入以下试剂：

c.混合均匀后封口膜封口，置于-80℃超低温冰箱1小时，并预冷离心机到4℃。随后，20000×g离心1小时，弃上清，用1mL预冷的75％乙醇清洗沉淀两次，20000×g离心5min。弃上清后再次瞬时离心，充分洗净残留的乙醇，开盖晾干至沉淀半透明。加入20μL 10mMTris-HCl，静置3min，取1μL 进行Qubit定量，其余4℃保存。

(6)Jurkat 20代和30代RRBS文库片段分布，文库纯化后，使用2100 AnalysisSysterm对文库的片段分布进行检测，得到文库的片段分布图，图1，如图1所示，可见回收的文库片段分布主峰在260bp左右，且无接头二聚体，符合实验及上机要求，用于上机进行高通量测序。

实施例6：数据分析

高通量测序后，得到的数据，采用二项分布近似贝叶斯的算法，建立了定量分类甲基化高保真位点的数据分析方法，数据分析步骤如下：

1.甲基化转化参数计算

利用未甲基化的lambda DNA计算实验样本的每个CpG位点的重亚硫酸盐平均转化效率CR，计算公式为S1，利用甲基化的lambda DNA计算实验样本的每个CpG位点的重亚硫酸盐平均错误转化率ICR，计算公式为S2， spike-in的重亚硫酸盐转化效率CR和重亚硫酸盐错误转化率ICR与实验样本相同，计算公式S1和S2分别为：

其中，U_c为未甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，U_mc为未甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目， M_c为甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，M_mc为甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目。两个实验样本 20代的CR和ICR平均值分别是99.2％和1.9％，30代的CR和ICR平均值分别是99.2％和1.8％。

2.计算实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平ML 和成功率θ

利用平均转化效率CR和平均错误转化率ICR计算实验样本的每个CpG 位点的重亚硫酸盐转化后的甲基化水平ML，spike-in的重亚硫酸盐甲基化水平与实验样本相同，重亚硫酸盐转化后的甲基化水平ML的计算公式为S3：

ML＝(CR-ICR)×ml+1-CR (S3)

其中，ML为实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平，ml为实验样本的每个CpG位点的重亚硫酸盐转化前的甲基化水平。由于传代培养的细胞均来源于同一个细胞，实验样本的重亚硫酸盐转化前的甲基化水平ml为100％、50％或0％，引入平均转化效率CR和平均错误转化率ICR 的参数，根据公式S3得到的重亚硫酸盐转化后的甲基化水平ML即为成功率θ。两个实验样本20代和30代的成功率θ值分别是98.1％、49.4％、0.8％和98.2％、49.5％、0.8％。

3.计算观测值Y

根据成功率θ计算观测值Y，利用二项分布似然函数S4进行模拟，设每个CpG位点的测序深度为n，即以成功率θ进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为观测值Y，计算观测值Y的公式为S4，p(x＝k)是事件成功的次数等于观测值Y时的概率，公式S4 为：

其中，k为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，θ为事件成功的概率即重亚硫酸盐转化后的甲基化水平ML，n为独立实验的次数即该CpG位点测序深度。

4.计算先验值X

由于每个CpG位点重亚硫酸盐转化后的甲基化水平ML为先验分布,先验概率为π(θ)，根据先验概率随机采样一个参数θi，设每个CpG位点的测序深度为n，即以参数θi进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为先验值X，计算先验值X的公式为S5， p(x|ML)是事件成功的次数等于先验值X时的概率，公式S5为：

其中，x为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，重亚硫酸盐转化后的甲基化水平ML为事件成功的概率，n为独立实验的次数即该CpG位点测序深度。

5.得到后验分布

将观测值Y与先验值X进行比较，若先验值X与观测值Y相等，则保留随机采样的参数θi作为后验分布的样本；若先验值X与观测值Y不相等，则舍弃，返回上一步重新选择一个参数θi，直到有N＝1000个样本被保留。

6.分类得到高保真甲基化位点

我们得到后验分布即所有甲基化保真位点的水平分布，用R语言的 quantile函数取95％的置信下限作为阈值，若CpG位点的甲基化水平大于该阈值则为保真位点，反之为动态位点，该保真位点即为高保真的甲基化CpG 位点。我们取20代和30代两个时期均为保真位点的甲基化位点作为真正意义上的高保真甲基化位点。

得到两个实验样本(Jurkat 20代和30代)的CpG位点甲基化水平散点图，如图2所示，将细胞培养到20代和30代的甲基化水平绘制到二维散点图上，CpG位点在(0％，0％)、(50％，50％)和(100％，100％)的坐标位置说明该位点随着细胞的传代甲基化水平保持稳定，而从图2中可以看出一部分CpG位点偏离了这三个点，比如，在20代时这个CpG位点的甲基化水平为100％，但是到了30代水平降到了80％左右，或者在20代时甲基化水平为0％，但是到了30代升到了20％，说明这些位点的甲基化水平发生了变化。因此我们比较两个时间点从而分离出随着细胞传代严格保持甲基化信息的高保真位点。

根据我们的统计模型，我们对CpG位点的甲基化水平进行计算，100％、 50％、0％高保真甲基化位点以及动态位点分类统计结果展示如下，其中100％高保真位点占总位点数比例为25.29％，50％高保真位点比例为6.52％，0％保真位点比例为24.78％。动态位点数目占总位点数的比例为43.4％，具体统计数据如下表：

上述得到的高保真的CpG甲基化位点，可以通过设计探针或构建敲除的细胞株，进一步研究这些高度保守的甲基化位点的特性和与疾病相关的DNA甲基化变化，对于理解疾病的发生和发展具有重要的意义。

尽管本申请已公开了多个方面和实施方式，但是其它方面和实施方式对本领域技术人员而言将是显而易见的，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。本申请公开的多个方面和实施方式仅用于举例说明，其并非旨在限制本申请，本申请的实际保护范围以权利要求为准。

Claims (3)

1.一种可传代细胞的高保真甲基化位点的分离与分析方法，其特征在于，包括步骤：

(A)将可传代细胞进行传代培养多代后，所述可传代细胞为哺乳动物细胞或细胞系，取不同代数的细胞进行基因组DNA提取，作为实验样本，将传代培养至20代-30代的细胞作为起始样本，然后进行血清饥饿24h，将各细胞样本同步化到G1/G0期后分别进行细胞收集和基因组DNA提取，作为实验样本；

(B)将spike-in加入实验样本中，将实验样本和spike-in分别进行DNA片段化处理，然后分别进行末端修复及加A、连接接头、USER酶消化及片段筛选纯化；再将甲基化和未甲基化的lambda DNA混合后按照0.5％～1％的比例加入到实验样本基因组DNA中进行重亚硫酸盐转化，扩增纯化后得到重亚硫酸盐测序文库，上机进行高通量测序，构建简并代表性重亚硫酸盐测序(RRBS)文库，并进行测序；所述spike-in为甲基化和未甲基化的lambda DNA；

(C)对于测序数据，采用高通量测序后得到的数据，采用二项分布近似贝叶斯的算法，建立了定量分类甲基化高保真位点的数据分析方法，数据分析步骤如下：

(1)甲基化转化参数计算

利用未甲基化的lambda DNA计算实验样本的每个CpG位点的重亚硫酸盐平均转化效率CR，计算公式为S1，利用甲基化的lambda DNA计算实验样本的每个CpG位点的重亚硫酸盐平均错误转化率ICR，计算公式为S2，spike-in的重亚硫酸盐转化效率CR和重亚硫酸盐错误转化率ICR与实验样本相同，计算公式S1和S2分别为：

其中，Uc为未甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，Umc为未甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目，Mc为甲基化的lambda DNA进行重亚硫酸盐转化的CpG位点数目，Mmc为甲基化的lambda DNA未进行重亚硫酸盐转化的CpG位点数目。两个实验样本20代的CR和ICR平均值分别是99.2％和1.9％，30代的CR和ICR平均值分别是99.2％和1.8％；

(2)计算实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平ML和成功率θ，利用平均转化效率CR和平均错误转化率ICR计算实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平ML，spike-in的重亚硫酸盐甲基化水平与实验样本相同，重亚硫酸盐转化后的甲基化水平ML的计算公式为S3：

ML＝(CR-ICR)×ml+1-CR (S3)

其中，ML为实验样本的每个CpG位点的重亚硫酸盐转化后的甲基化水平，ml为实验样本的每个CpG位点的重亚硫酸盐转化前的甲基化水平。由于传代培养的细胞均来源于同一个细胞，实验样本的重亚硫酸盐转化前的甲基化水平ml为100％、50％或0％，引入平均转化效率CR和平均错误转化率ICR的参数，根据公式S3得到的重亚硫酸盐转化后的甲基化水平ML即为成功率θ；

(3)计算观测值Y

根据成功率θ计算观测值Y，利用二项分布似然函数S4进行模拟，设每个CpG位点的测序深度为n，即以成功率θ进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为观测值Y，计算观测值Y的公式为S4，p(x＝k)是事件成功的次数等于观测值Y时的概率，公式S4为：

其中，k为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，θ为事件成功的概率即重亚硫酸盐转化后的甲基化水平ML，n为独立实验的次数即该CpG位点测序深度；

(4)计算先验值X

由于每个CpG位点重亚硫酸盐转化后的甲基化水平ML为先验分布,先验概率为π(θ)，根据先验概率随机采样一个参数θi，设每个CpG位点的测序深度为n，即以参数θi进行n次的伯努利实验，计算得到每次实验成功的次数即每个CpG位点胞嘧啶的数量为先验值X，计算先验值X的公式为S5，p(x|ML)是事件成功的次数等于先验值X时的概率，公式S5为：

其中，x为事件成功的次数，即该CpG位点检测到的甲基化序列的数量，重亚硫酸盐转化后的甲基化水平ML为事件成功的概率，n为独立实验的次数即该CpG位点测序深度；

(5)得到后验分布

将观测值Y与先验值X进行比较，若先验值X与观测值Y相等，则保留随机采样的参数θi作为后验分布的样本；若先验值X与观测值Y不相等，则舍弃，返回上一步重新选择一个参数θi，直到有N＝1000个样本被保留；

(6)分类得到高保真甲基化位点

将得到后验分布即所有甲基化保真位点的水平分布，用R语言的quantile函数取95％的置信下限作为阈值，若CpG位点的甲基化水平大于该阈值则为保真位点，反之为动态位点，该保真位点即为高保真的甲基化CpG位点。

2.如权利要求1所述的可传代细胞的高保真甲基化位点的分离与分析方法，其特征在于，在步骤(A)中，来自于生物组织、含有细胞的体液；所述传代培养的细胞是同源且表型单一的。

3.如权利要求2所述的可传代细胞的高保真甲基化位点的分离与分析方法，其特征在于，从起始样本中选择2-4个培养了不同代数的细胞样本，将各细胞样本同步化到G1/G0期后进行细胞收集和基因组DNA提取，作为实验样本。