cfDNA与甲基化
写在前面的
cfDNA相关研究大部分会经常提到两个关键词：Methylation和Bisulfite。前者是甲基化的意思，后者是检测甲基化常用的测序方法中的一种重要材料亚硫酸盐。表观遗传学领域研究发现哺乳动物基因组中常常出现甲基化修饰，一般体现在胞嘧啶C的5号位碳原子在甲基转移酶的作用下结合一个甲基基团（CH3）形成5mC。甲基化修饰被证明会影响基因的表达，从而调控细胞命运进而影响疾病的发生。疾病产生早期，无论是细胞的自主裂解还是病变细胞被免疫系统识别后的被动裂解都能释放出cfDNA，释放的cfDNA在一定程度上记录了疾病产生的表观遗传信息。因此如果能及时识别出cfDNA上携带的表观遗传信息则对于疾病早筛具有重要意义。因此cfDNA相关研究中甲基化表观遗传信息挖掘是一个重要的组成部分。
基于重亚硫酸盐的甲基化测序
除了传统的基于PCR式的甲基化检查技术，常规的甲基化检测技术包含基于Bisulfite转化的检测技术、基于限制性内切酶的检测技术和通过免疫共沉淀的检测技术。
这些技术中基于限制性内切酶的检测技术也称为Methylation sensitive restriction endonuclease, MRSE，这种技术通常会使用对甲基化位点敏感的酶进行切割序列。目前大约有320种酶对甲基化位点敏感。这类酶的工作原理是如果酶切位点的碱基发生甲基化，则酶无法切割该位点，通过DNA印迹法分析即可完成对甲基化位点的识别。这种技术的优势在于无需了解靶标DNA序列的信息，可以直接提供甲基化位点的评价。但是缺点在于对于核酸质量有较高要求，并且受限于酶切序列的特异性，无法做到全基因组水平的甲基化位点识别。
基于免疫共沉淀的甲基化检测技术通常使用5’-甲基胞嘧啶抗体特异性富集基因组上发生甲基化的DNA片段，然后利用高通量测序技术对全基因组CpG岛高度甲基化的区域进行分析，但是受制于分辨率问题，该方法对单碱基位点的甲基化识别能力较弱。
基于重亚硫酸盐（Bisulfite based）测序的技术巧妙在，重亚硫酸盐可以对未发生甲基化的胞嘧啶发生脱氨基反应进而转化为尿嘧啶（下图所示），已经发生甲基化的胞嘧啶则无法脱氨基生成尿嘧啶。经过反应后，理论上剩余的胞嘧啶大部分来源于发生甲基化的胞嘧啶。通过高通量测序后，则能识别出甲基化。但是重亚硫酸盐测序有两种缺点：1.胞嘧啶5号为发生羟基化修饰的胞嘧啶同样可以在重亚硫酸盐环境下保留，尽管甲基化修饰是较常见的，但是现有的文献发现5-甲基修饰和5-羟基修饰具有不同的生物学意义。2.重亚硫酸盐处理后的样本中大部分DNA容易被降解，测序后序列的复杂度显著降低。尽管重亚硫酸盐测序存在一些弊端，但是这种方法的简便性和成本以及相关生物信息分析方法的完善驱动了重亚硫酸盐测序在过去的几十年中是甲基化测序行业的金标准技术。
重亚硫酸盐测序策略
RRBS
RRBS是Reduced Representation Bisulfite Sequencing的缩写，RRBS技术通常使用对甲基化碱基不敏感的MspI酶，该酶的识别切割靶标为5’-C^CGG-3’。根据酶切位点的特性，可以推出这种酶通常在CpG岛富集区间具有良好的效果。因为人类基因组上CpG岛存在两种分布方式，第一种是主要的分布方式，CpG岛集中分布于启动子区域或在基因内部充当替代启动子。第二种分布方式为离散地分布在基因间，这种分布方式较少。因此MspI酶主要靶向于CpG岛富集区。当使用MspI酶切后，特异性地富集启动子及其CpG岛区域就能靶向性地获得CpG岛的甲基化图谱。根据RRBS的策略，可以推出这种方式成本低，但是缺点在于仅能识别CpG岛富集区，对于基因组其他区域的甲基化识别能力受限。
WGBS
WGBS是Whole Genome Bisulfite Sequencing的缩写，顾名思义WGBS技术是面向全基因组水平的高通量甲基化检测技术，这种技术不仅可以识别CpG岛富集区的甲基化水平，还具备将CpG岛离散区的少量甲基化位点检出的能力。2008年后随着MethylC-seq测序技术的革新，WGBS逐渐成为了甲基化测序领域的金标准技术。
重亚硫酸盐测序原理
在理解重亚硫酸盐测序原理之前，有一个基本的概念需要阐述。特别是针对我这种博士期间以病毒为主要研究对象的外行而言，这个概念的理解对于我更好地理解该领域的下游知识非常重要。DNA分子由两条链构成，以前仅仅只是简单地称呼其双链分别为正链和负链。但是还有一种命名方式。为了纪念Watson和Crick发现DNA双螺旋结构，通常情况下将转录为mRNA的那条链，也就是常说的反义链（antisense strain）或者模板链（template strain）也就是我以前常听说负链称为Watson链，而和mRNA序列一致的那条链，也就是正义链（sense strain）或编码链（coding strain），也就是我以前常听说的正链称为Crick链。如下图所示，我们假设有一条双链DNA，其中Watson链和Crick链各有一个胞嘧啶发生甲基化，记为Cm，首先使用高温变性DNA，再利用Bisulfite处理变性后的单链DNA，在Bisulfite的作用下未发生甲基化的胞嘧啶发生脱氨基反应变成尿嘧啶（U），两个甲基化的胞嘧啶仍然是胞嘧啶。此时，Bisulfite环境下发生序列改变的两条链名称发生了变化，Watson变成了Bisulfite Watson（BSW缩写），Crick变成了Bisulfite Crick（BSC缩写），由于构建测序文库前需要经过PCR扩增核酸量，因此在PCR后根据碱基互补原则，BSW的互补链中原本应该添加为鸟嘌呤（G）的位点全部在尿嘧啶（U）的作用下变成了腺嘌呤（A），此时BSW的互补链Bisulfite Watson Reverse Complementation，得到了一个缩写名BSWC，而BSC的互补链Bisulfite Crick Reverse Complementation（BSCR）中原本添加为鸟嘌呤（G）的位点也变成了腺嘌呤（A），所以在Bisulfite的作用下，BSW和BSC上的C全部变成了U，对应于测序数据中C>T转换，而BSCR和BSWR中的G全部变成了A，对应测序数据中G>A转换。因此重亚硫酸盐测序直接引入了C>T及G>A的转换，影响了测序数据中4碱基的分布。
重亚硫酸盐测序数据分析策略
由于Bisulfite的改变，BS-seq数据中reads的碱基被人为改变，那么传统的序列比对算法则无法处理BS-seq数据。因此BS-seq的序列比对相比于传统NGS数据分析的算法增加了复杂度。为了解决这个问题生物信息学家提出了2种适合于BS-seq序列比对的算法：
三碱基比对算法
三碱基比对算法在序列比对前会人为进行碱基转换，具体操作为分别对reads和参考基因组各进行两次C>T和G>A的变换。转换后再使用传统的序列比对算法进行比对。比对前处理操作如下图所示
为了便于表述，将U统一称为T，与NGS测序数据中的碱基一致。对于BSW链，未经Bisulfite处理时，原始链和基因组一致，经过处理未甲基化的C变为T，甲基化的C保留，所以此刻只需将甲基化的C转变为T即可与经过C>T转换的参考基因组（Watson链）一致。反之，当一条reads经过C>T转换时与经过C>T转换的参考基因组有最佳对齐结果则说明该reads为BSW。
对于BSWR链，由于BSWR链是BSW链的互补链，BSW链被处理前未甲基化的C被转换为T，这些T在BSWR链中全部变为A，相当于参考基因组的互补链（Crick链）对应位置的G变为了A，相当于参考基因组的链对应位置的C变为了T。对于BSW甲基化的C，在BSWR中仍然互补配对为G，相当于这些G位点仍然和参考基因组的互补链的G对应，也就是仍然与参考基因组链的C互补。因此BSWR链经过G>A处理后理论上与参考基因组经过C>T处理后序列互补。反之，当一条reads经过G>A转换，再互补与参考基因组经过C>T处理后具有最好的对齐，则说明该reads为BSWR。
对于BSC链，由于BSC链是BSW链的互补链，因此BSC链在被处理前与参考基因组的互补链一致，BSC处理前未甲基化的C在处理后变为T，相当于参考基因组互补链对应位置的C变为T，相当于参考基因组链对应位置的G变为A，对于BSC链甲基化的C仍然被保留，这些C仍然于参考基因组链对应位点互补。因此BSC链经过C>T处理后理论上与参考基因组经过G>A处理后序列互补。反之，当一条reads经过C>T转换，再互补与参考基因组经过G>A处理后具有最好的对齐，则说明该reads为BSC。
对于BSCR链，由于BSCR链来自于BSC链的互补链，BSC链又是BSW的互补链，所以BSCR链方向与BSW链一致。BSC链经过甲基化处理后未甲基化的碱基C（对应BSW处理前的G）全部变成T，因此BSCR相应于BSW处理前的这些位点全部变成了A，BSC链甲基化的碱基C保留，并在BSCR变为G，和BSW处理前的G一致，所以BSCR链将自身的G全部变成A后可以与参考基因组经过G>A替换的序列一致。反之，当一条reads经过G>A转换时与经过G>A转换的参考基因组有最佳对齐结果则说明该reads为BSCR。
当完成碱基转换后找到每条reads最佳对齐位置后再回溯原本的碱基识别潜在的甲基化位点。
通配符比对算法
与三碱基比对算法不同，通配符比对算法直接默认参考基因组的C可以与reads的C或T配对，而参考基因组的其他碱基只能与reads中相对应碱基配对。这是因为根据Bisulfite测序原理，当reads中存在的T可能是未甲基化的C，也可能是经过PCR从祖先序列中的T配对而来的。而reads中存在的C可能是甲基化的C，也可能是经过PCR从祖先序列中的G配对而来的。参考基因组中无论是Watson链还是Crick链只有C碱基都对应了两种可能，尽管还存在G>A的变化，但是G>A并不是主动发生的，而是C>T变化后通过互补配对原则伴随发生的。因此通配符比对算法相对于三碱基比对更直观，它运行C可以与C或T配对，这样reads和参考基因组的比对就变成了reads和Watson链或Crick链的在允许一定容错（C/T）下的最佳配对选择。比对情况则是将任意一条reads直接与Watson或Crick链比对，或者反向互补后与Watson或Crick链比对。选择最佳配对结果即可确定属于什么链。下图大概描述了通配符配对算法下不同reads与参考基因组最佳对齐结果。完成序列比对后即可根据匹配的C碱基进行潜在甲基化位点的筛选。
通过上面算法原理的简单分析，三碱基算法是通过人为引入替换将重亚硫酸盐测序数据比对转换为普通NGS比对，这种策略化未知为已知，但是缺点在于分析速度较慢。而通配符比对的方法显著提升了速度，但是也有可能因为测序的错误率提高假阳性甲基化位点的检出概率。
参考
Jiang P, Sun K, Lun FM, Guo AM, Wang H, Chan KC, Chiu RW, Lo YM, Sun H. Methy-Pipe: an integrated bioinformatics pipeline for whole genome bisulfite sequencing data analysis. PLoS One. 2014 Jun 19;9(6):e100360. doi: 10.1371/journal.pone.0100360. PMID: 24945300; PMCID: PMC4063866.
Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics. 2009 Jul 27;10:232. doi: 10.1186/1471-2105-10-232. PMID: 19635165; PMCID: PMC2724425.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011 Jun 1;27(11):1571-2. doi: 10.1093/bioinformatics/btr167. Epub 2011 Apr 14. PMID: 21493656; PMCID: PMC3102221.
写在后面的
甲基化测序在cfDNA及其疾病早筛过程中具有非常重要的价值。相关研究非常多仍然需要不断学习。