DNA甲基化测序数据处理（一）：数据比对

前言

因为组里面出了一批甲基化测序数据，使用的技术为BS-seq，处理的时候顺带记录了学习过程，演示使用数据为官方提供的example.fastq。

DNA甲基化及CpG岛定义了解

DNA甲基化（英语：DNA methylation）为DNA化学修饰的一种形式，能在不改变DNA序列的前提下，改变遗传表现。为外遗传编码（epigenetic code）的一部分，是一种外遗传机制。DNA甲基化过程会使甲基添加到DNA分子上，例如在胞嘧啶环的5'碳上：这种5'方向的DNA甲基化方式可见于所有脊椎动物。
在人类细胞内，大约有1%的DNA碱基受到了甲基化。在成熟体细胞组织中，DNA甲基化一般发生于CpG双核苷酸（CpG dinucleotide）部位；而非CpG甲基化则于胚胎干细胞中较为常见[1] [2]。植物体内胞嘧啶的甲基化则可分为对称的CpG（或CpNpG），或是不对称的CpNpNp形式（C与G是碱基；p是磷酸根；N指的是任意的核苷酸）。
特定胞嘧碇受甲基化的情形，可利用亚硫酸盐定序（bisulfite sequencing）方式测定。DNA甲基化可能使基因沉默化，进而使其失去功能。此外，也有一些生物体内不存在DNA甲基化作用。
——维基百科

DNA甲基化作为基因组上的表观修饰（区别于组蛋白修饰），存在于各种生物中。

简单来说，可以认为甲基化代表着基因的失活，去甲基化则标志基因的激活与表达

虽然CpG序列出现的频率并不高，但是在某些基因区域内，CpG的密度很高，俗称CpG岛。这些CpG岛大多出现在基因的启动子区域（人类占到70%），长度达300-3000bp。目前的研究表明，大多数的管家基因都含有CpG岛，位于基因的5'端（其中的大多数CpG岛都是未甲基化的）。

另外需要注意的是，目前的研究表明，肿瘤样本与正常样本的CpG岛甲基化差异大多不是发生CpG岛的内部而是位于CpG岛岸（CpG island shore）。

由于CpG位点的易甲基化导致胞嘧啶脱氨变成胸腺嘧啶，所以在漫长的进化过程中，CpG位点逐渐消失，但是又存在着对于基因表达的调控要求，所以CpG岛的出现也被理解为抵抗甲基化经常很，维持调控功能。

DNA甲基化测序原理与方法

此处略过，请自行了解(示例文件为WGBS单端测序文件)。

其实是因为理解起来比较累，知道大致原理就可以了，以后再来填坑

数据处理（使用Bismark软件处理）

Bismark下载

Bismark官网

The less the people know about how sausages and our code are made, the better they sleep at night (untracable author)
bismark_v0.19.0.tar.gz
解压即用tar xzf bismark_v0.X.Y.tar.gz，常用的话将路径写入PATH:PATH=/PATH/TO/Bismark/:$PATH
manual页面需要搭梯子才能看到，简书不支持附件，所以有需要的可以留言，我email pdf文件

Dependencies

需要用户已经装好bowtie1/bowtie2

测试数据获取

此处使用测试数据test.fastq
（from SRR020138, Lister et al., 2009; trimmed to 50 bp; base call qualities are Sanger encoded Phred values (Phred33)).

$ls -l total 2.1M -rw-r--r-- 1 sxj users 2.1M Apr 17 2018 test_data.fastq

INDEX文件构建

# under install folder
./bismark_genome_preparation --path_to_bowtie /PATH/to/bowtie2/ --verbose ~/PATH/to/GRCh38/

比对

# paired-end data ./bismark --genome ~/PATH/to/GRCh38/ -1 read1.fastq.gz -2 read2.fastq.gz -p 4 -o ./ 2>test.log # single-end data ./bismark --genome ~/PATH/to/GRCh38/ test_data.fastq -p 4 -o ./ 2>test.log

甲基化位点提取

./bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder ~/PATH/to/GRCh38/ test_data_bismark_bt2.bam 2>extracor.log

生成处理报告

./bismark2report

所有结果文件

ls -l
-rw-r--r-- 1 sxj users  82K Apr 18 16:11 CHG_context_test_data_bismark_bt2.deduplicated.txt.gz
-rw-r--r-- 1 sxj users 154K Apr 18 16:11 CHH_context_test_data_bismark_bt2.deduplicated.txt.gz
-rw-r--r-- 1 sxj users 138K Apr 18 16:11 CpG_context_test_data_bismark_bt2.deduplicated.txt.gz
-rw-r--r-- 1 sxj users 105K Apr 18 17:30 extracor.log
-rw------- 1 sxj users  20K Apr 17 15:56 nohup.out
-rw-r--r-- 1 sxj users 450K Apr 17 15:56 test_data_bismark_bt2.bam
-rw-r--r-- 1 sxj users 449K Apr 18 10:34 test_data_bismark_bt2.deduplicated.bam
-rw-r--r-- 1 sxj users  48K Apr 18 16:11 test_data_bismark_bt2.deduplicated.bedGraph
-rw-r--r-- 1 sxj users  55K Apr 18 16:11 test_data_bismark_bt2.deduplicated.bismark.cov
-rw-r--r-- 1 sxj users 217M Apr 18 17:30 test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt.gz
-rw-r--r-- 1 sxj users 2.9K Apr 18 16:11 test_data_bismark_bt2.deduplicated.M-bias.txt
-rw-r--r-- 1 sxj users  766 Apr 18 16:11 test_data_bismark_bt2.deduplicated_splitting_report.txt
-rw-r--r-- 1 sxj users  260 Apr 18 10:34 test_data_bismark_bt2.deduplication_report.txt
-rw-r--r-- 1 sxj users 347K Apr 19 23:18 test_data_bismark_bt2_SE_report.html
-rw-r--r-- 1 sxj users 1.8K Apr 17 15:56 test_data_bismark_bt2_SE_report.txt
-rw-r--r-- 1 sxj users 2.1M Apr 17 15:14 test_data.fastq
-rw-r--r-- 1 sxj users 6.2K Apr 17 15:56 text.log

结果解读

--cytosine_report参数会根据当前目录下的信息文件生成一个HTML格式的报告文件，即test_data_bismark_bt2_SE_report.html文件，它包括了比对信息，甲基化信息，M-bias等，可以对数据有一个大概的认知（下图只展示了一部分）：

比对信息

M-bias

同时因为使用了--comprehensive，所以结果合并正反链的数据后会输出CpG/CHG/CHH三种类型的甲基化文件，包含了胞嘧啶所有的组合形式，但实际上我们自然最关注的是CpG位点的甲基化。其中

CpG_context_test_data_bismark_bt2.deduplicated.txt即CpG甲基化位点的文件。

head CpG_context_test_data_bismark_bt2.deduplicated.txt Bismark methylation extractor version v0.19.0 SRR020138.15024317_SALK_2029:7:100:1672:902_length=86 - 1 57333019 z SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 + 2 10026473 Z SRR020138.15024331_SALK_2029:7:100:1673:1282_length=86 + 11 78025243 Z SRR020138.15024343_SALK_2029:7:100:1673:202_length=86 + 10 121617231 Z SRR020138.15024357_SALK_2029:7:100:1673:879_length=86 - 4 75173715 z SRR020138.15024361_SALK_2029:7:100:1673:235_length=86 - 2 130768889 z SRR020138.15024368_SALK_2029:7:100:1673:123_length=86 + 10 106402850 Z SRR020138.15024376_SALK_2029:7:100:1673:1908_length=86 - 12 124097382 z SRR020138.15024380_SALK_2029:7:100:1673:397_length=86 + 8 95038627 Z # 第一列为测序信息 # 第二列为甲基化状态 + 代表甲基化 -代表未甲基化 # 第三列代表chromosome # 第四列代表location # 第五列代表methylation call，简单来说大写的就是甲基化的（因为还有CHG，CHH的数据，分别对应x, X , h, H）

test_data_bismark_bt2.deduplicated.bismark.cov文件则给了每个位点的甲基化比例，为下一步确定CpG岛提供了基础，其数据形式如下：

$head test_data_bismark_bt2.deduplicated.bismark.cov
1   975476  975476  100 1   0
1   975488  975488  100 1   0
1   975490  975490  100 1   0
1   2224487 2224487 100 1   0
1   2224489 2224489 100 1   0
1   2224514 2224514 100 1   0
1   2224520 2224520 100 1   0
1   2313220 2313220 0   0   1
1   9313902 9313902 100 1   0
1   9313914 9313914 100 1   0

# 第一列代表chromosome
# 第二，三列代表location
# 第四列代表甲基化百分比
# 第五列代表甲基化数目
# 第六列代表未甲基化数目

test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt文件则是背景信息：

$head test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt 1 10469 + 0 0 CG CGC 1 10470 - 0 0 CG CGA 1 10471 + 0 0 CG CGG 1 10472 - 0 0 CG CGC 1 10484 + 0 0 CG CGG 1 10485 - 0 0 CG CGG 1 10489 + 0 0 CG CGC 1 10490 - 0 0 CG CGG 1 10493 + 0 0 CG CGC 1 10494 - 0 0 CG CGG # 第一列为chromosome # 第二列为location # 第三列为strand # 第四，五列为甲基化和非甲基化的碱基数目 # 第六列为CG # 第七列为背景（第一个C延伸两个碱基）

其它参数

bam2nuc
bismark2summary
coverage2cytosine
NOMe_filtering
filter_non_conversion
# 有需要的可以自信探索 --help or manual

结语

此处根据测序数据得到了甲基化位点的信息，但是后续DML以及DMR的确定还需要R包的使用，以及后续的可视化还以探索以下包：

挖坑待填

Methy-Pipe: An Integrated Bioinformatics Pipeline for Whole Genome Bisulfite Sequencing Data Analysis
2014年出的一个pipeline，分析作图一条龙，有兴趣的同学安排一下。