技术解读｜RRBS测序中因酶切人为引入碱基问题
前言
我们都知道，研究哺乳动物DNA甲基化最常用的、性价比最高的技术就是RRBS技术了。RRBS技术是基于酶切富集CG位点的，在建库过程中往往会因补平酶切的粘性末端而人为引入碱基。这时候酶切处引入的人工C碱基如果不处理的话，将会影响这些位点甲基化水平的计算。今天小编将会从酶切引入碱基的过程入手，带大家一步步分析其中的缘由，最后给出解决人工引入碱基问题的办法。
在常规RRBS建库过程中，在补平酶切接口时会引入人工的CG序列，这两个序列在read1中CT转化会变成3'TG,在read2中GA转化会变成5'CA，故需要去除，以免影响这个位置C碱基甲基化水平的计算。
在双酶切RRBS中，在补平酶切接口时会引入人工的CG和CWG序列，read1中CT转化会变成3'TG和3'TWG,在read2中GA转化会变成5'CA和5'CWA，故需要去除，以免影响这个位置C碱基甲基化水平的计算。
为什么是如此去除？具体分析见下文：
单酶切RRBS
单酶切RRBS技术采用的是MspI 单酶切建库，MspI的酶切位点如下：
我们回顾一下单酶切RRBS建库过程：
我们分析一下整个建库过程中酶切位点周围的碱基变化：
原序列：
5'-C|CGG=========C|CGG-3' （W链）3'-GGC|C=========GGC|C-5' （C链）
酶切后,3'端补平CG，再加A,然后加接头：(红色字体为人工碱基)
5'-P5adapter-T|CGG=========C|CGA-P7adapter-3' （W链）3'-P7adapter-AGC|C=========GGC|T-P5adapter-5' （C链）
BS处理：
5'-P5adapter-T|TGG=========T|TGA-P7adapter-3' (BSW链)3'-P7adapter-AGT|T=========GGT|T-P5adapter-5' (BSC链)
PCR：
5'-P5adapter-T|TGG=========T|TGA-P7adapter-3' (BSW链)---read13'-P5adapter-AAC|T=========AAC|T-P7adapter-5' (BSWR链)---read2
5'-P7adapter-T|CAA=========C|CAA-P5adapter-3' (BSCR链)---read23'-P7adapter-AGT|T=========GGT|T-P5adapter-5' (BSC链)---read1
由此可见，最后测序得到序列中，人工序列有(不计算末端A和T)：
read1：
3'TG
read2：
5'CA 3'CA
是不是都要去除？
实际上read2的3'CA虽然是人工序列，却并不影响原始链（W链和C链）中5'CGG中C碱基甲基化水平的计算（请仔细推敲）。这也是为什么酶切了这段DNA却不影响这几个位点甲基化水平计算的原因。
故，单酶切RRBS数据预处理，在去除接头污染之后，还要去除：
read1的3'TG
和
read2的5'CA
cutadapt实现：
cutadapt \
-a AGATCGGAAGAGC -A AGATCGGAAGAGC \
-a TG$ -G ^CA \
-n 2
-o trimmed.1.fastq.gz -p trimmed.2.fastq.gz \
reads.1.fastq.gz reads.2.fastq.gz
说明：
小写-a和-g表示切read1接头，大写-A和-G表示切read2接头
-a 表示切3'接头，-g表示切5'端接头
-n 表示对read最多需要剪切2次
前两个接是illumina通用接头的共同前缀部分，红色字体的A正是表示切补平后末端加的A
后面两个是上文说的RRBS需要切掉的碱基。
双酶切RRBS
有的同学可能没听过双酶切RRBS，说到这个技术小编不免瞬间自豪感爆棚，这是易基因现任实验研发部门负责人领衔开发的技术哦！它是一种基于常规RRBS开发的技术，相比单酶切RRBS，能覆盖C位点更多更全面，性价比更高！关于双酶切技术的论文可以到本文末尾查看。
双酶切RRBS采用的是MspI 和 ApeKI双酶切建库，这两种酶的酶切位点如下：
W表示A或T
我们分析一下整个建库过程中酶切位点周围的碱基变化：
双酶切可切三种片段：
1. 单MspI 酶切
5'-C|CGG=========C|CGG-3' （W链）3'-GGC|C=========GGC|C-5' （C链）
2.单ApeKI酶切
5'-G|CWGC=========G|CWGC-3' （W链）3'-CGWC|G=========CGWC|G-5' （C链）
3.MspI+ApeKI双酶切
5'-G|CWGC=========C|CGG-3' （W链）3'-CGWC|G=========GGC|C-5' （C链）
4.ApeKI+MspI双酶切
5'-C|CGG=========G|CWGC-3' （W链）3'-GGC|C=========CGWC|G-5' （C链）
不考虑加A的影响，
第二种的变化为：
补平：
5'-CWGC=========GCWG-3' （W链）3'-GWCG=========CGWC-5' （C链）
BS转化：
5'-TWGT=========GTWG-3' （BSW链）3'-GWTG=========TGWT-5' （BSC链）
PCR:
5'-TWGT=========GTWG-3' （BSW链） ---read13'-AWCA=========CAWC-5' （BSWR链）---read2
+
5'-CWAA=========ACWA-3' （BSCR链）---read23'-GWTG=========TGWT-5' （BSC链) ---read1
同样，人工加入的序列有：
read1的：3'TWG
read2的：5'CWA 和3'CWA
但是由于read2的3'CWA是原始非人工序列互补而来的，所以它不影响原始W和C链5'CWG中C碱基甲基化水平的计算。
所以需要去除的序列是：
read1的：3'TWG
read2的：5'CWA
第三种的变化为：
补平：
5'-CWGC=========CCG-3' （W链）3'-GWCG=========GGC-5' （C链）
BS转化：
5'-TWGT=========TTG-3' （BSW链）3'-GWTG=========GGT-5' （BSC链）
PCR:
5'-TWGT=========TTG-3' （BSW链） ---read13'-AWCA=========AAC-5' （BSWR链）---read2
+
5'-CWAC=========CCA-3' （BSCR链）---read23'-GWTG=========GGT-5' （BSC链) ---read1
同样，人工加入的序列有：
read1的：3'TG 和 3'TWG
read2的：5'CWA 和CA，3'CWA 和CA
但是由于read2的3'CWA和3'CA不影响原始W链5'CWG和C链5'CG中C碱基甲基化水平的计算。
所以需要去除的序列是：
read1的：3'TG 和 3'TWG
read2的：5'CWA和 5'CA
第四种PCR后得到的序列为：
5'-TGG=========GTWG-3' （BSW链） ---read13'-ACC=========CAWC-5' （BSWR链）---read2
+
5'-CAA=========ACWA-3' （BSCR链） ---read23'-GTT=========TGWT-5' （BSC链）---read1
同样，人工加入的序列有：
read1的：3'TG 和 3'TWG
read2的：5'CWA 和CA，3'CWA 和CA
但是由于read2的3'CWA和3'CA不影响原始W链5'CWG和C链5'CG中C碱基甲基化水平的计算。
所以需要去除的序列是：
read1的：3'TG 和 3'TWG
read2的：5'CWA和 5'CA
总结以上分析结果，得出结论：
双酶切RRBS需要去除的序列是：
read1的：3'TG 和 3'TWG
read2的：5'CWA和 5'CA
cutadapt实现：
cutadapt \
-a AGATCGGAAGAGC -A AGATCGGAAGAGC \
-a TG$ -a TAG$ -a TTG$\
-G ^CA -G ^CAA -G ^CTA\
-n 2
-o trimmed.1.fastq.gz -p trimmed.2.fastq.gz \
reads.1.fastq.gz reads.2.fastq.gz
说明：
小写-a和-g表示切read1接头，大写-A和-G表示切read2接头
-a 表示切3'接头，-g表示切5'端接头
-n 表示对read最多需要剪切2次
前两个接是illumina通用接头的共同前缀部分，红色字体的A正是表示切补平后末端加的A
后面一个是上文说的双酶切RRBS需要切掉的碱基，每一行所列要切的序列之间是“或”的关系，故最多切2次，n设置为2。
双酶切RRBS技术论文：Wang et al. BMC Genomics 2013, 14:11 http://www.biomedcentral.com/1471-2164/14/11