全基因组甲基化分析简述

DNA甲基化是一种非常基础且重要的表观修饰，在调控基因表达、转录因子结合和抑制转座子元件中起到关键的作用。

目前，DNA甲基化检测的技术已经比较成熟，例如高通量的WGBS、RRBS、MeDIP-seq、MBD-seq，低通量的BSP、MSP等，其中以WGBS（Whole genome bisulfite sequencing）最为经典。WGBS可以单碱基分辨率尺度下在全基因组范围内鉴定和定量甲基化状态。然而，如何分析WGBS数据通常是一种挑战，涉及众多的分析步骤和注意事项。

以下简要介绍WGBS数据的分析步骤，主要包括：

数据清洗、质量检查和比对

DNA甲基化水平评估

差异甲基化

数据可视化

一、第1步 - trim reads

数据比对之前需要对接头和低质量碱基序列进行去除。

处理工具非常的多，这里以cutadapt为例：

单端数据：

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \

-m 50 -q 10,10 reads.fastq  \

> trimmed_reads.fastq

双端数据：

cutadapt -a AGATCGGAAGA -A CTCTTCCGATC -q 10,10 \

-o trimmed_read1.fastq -p trimmed_read2.fastq \

read1.fastq read2.fastq

压缩数据以节省硬盘存储：

gzip trimmed_reads.fastq

gzip trimmed_read1.fastq trimmed_read2.fastq

详细使用说明参考：http://cutadapt./en/stable/

二、第2步 - 数据质量评估

为了评估测序reads的质量，可以使用fastQC来进行，可以轻松的得到reads的质量评估报告。

使用同样非常简单且网页版质量评估简单易懂。

单端数据：

fastqc trimmed_reads.fastq.gz

双端数据：

fastqc trimmed_read1.fastq.gz trimmed_read2.fastq.gz

详细的说明参考：http://www.bioinformatics./projects/fastqc/

fastQC

三、第3步 - 基因组比对

下载物种的基因组序列，并建立基因组索引，通常WGBS建库需要加入内参lambda DNA序列。

首先，需要合并2个基因组序列：

cat genome.fa lambda.fa > genome_lambda.fa

对于基因组比对，这里推荐BS-Seeker2 [https:///10.1186/1471-2164-14-774]

其与经典的bismark相比，在比对率、精确率等方面均略胜一筹，同时包含甲基化检测套件，使用非常方便，具体信息各位可以阅读原文。

BS-seeker2

建立参考基因组索引：

python bs_seeker2-build.py -f genome_lambda.fa \

--aligner=bowtie2

随后，reads同索引后的参考基因组进行比对：

单端数据：

python bs_seeker2-align.py -g genome_lambda.fa \

--aligner=bowtie2 \

-u unmapped.fa \

-o mapped.bam \

--bt2-p \

-i trimmed_reads.fastq.gz

双端数据：

python bs_seeker2-align.py -g genome_lambda.fa \

--aligner=bowtie2 \

-u unmapped.fa \

-o mapped.bam \

--bt2-p \

-1 trimmed_read1.fastq.gz \

-2 trimmed_read2.fastq.gz

比对完成后对生成的BAM文件进行排序：

samtools sort -@ -T temp\

-O bam mapped.bam sorted

PCR重复reads是必须要去除的：

这里使用picard tools，samtools也可以但是不是很准确

java -jar picard.jar MarkDuplicates I=sorted.bam \

O=filtered.bam \

M=duplicate_stats.txt \

REMOVE_DUPLICATES=true \

AS=true

四、第4步 - DNA甲基化鉴定

基于全基因组胞嘧啶位置的reads的C/T比对情况进行检测，得到每个胞嘧啶的覆盖深度、序列背景等信息。

python bs_seeker2-

call_methylation.py -i filtered.bam \

--sorted -o \

--db

为了保证甲基化水平评估的准确性，需要对内参DNA进行重亚硫酸氢盐处理的转化率进行计算，需要保证转化率在98%以上。

得到全基因组范围内单个胞嘧啶的甲基化比率后即可对全基因组、染色体、功能区域、重复序列、基因等进行甲基化水平的计算，以进一步研究该物种的甲基化模式。

五、第5步 - 差异甲基化

当进行多个样本比较时，需要进行差异甲基化胞嘧啶（DMC）、差异甲基化区域（DMR）等分析，以寻找样本间特异、特有的甲基化模式。

这里推荐DSS进行差异位点和差异区域的寻找。

首先，需要准备其所需要的文件

格式如下：

例如，只对CpG类型的胞嘧啶进行差异分析：

awk 'BEGIN {FS=OFS=”\t”} {if ($4 == “CG”) print $1,

$3, $7, $8-$7}’ sample.CGmap > cg_positions.tsv

R软件环境下进行运行：

library(DSS)

column_names <- c(“chr”, “pos”, “N”, “X”)

 condition1 <- read.table('condition1_cg.tsv’,

header=column_names)

 condition2 <- read.table('condition2_cg.tsv’,

header=column_names)

 experiment <- makeBSseqData(list(condition1,

condition2), c('C1’, 'C1’))

dmlTest <- DMLtest(experiment, group1=c('C1’),

+group2=c('C2’)

dml <- callDML(dmlTest, delta=0.1,

+p.threshold=0.001)

如此，即可得到差异甲基化的CpG位点。

上述dml对象可以进一步用于更大区域的差异分析，即DMR分析：

dmrs <- callDMR(dmlTest, delta=0.1, minCG=3

+p.threshold=0.001, minlen=50,

+dist.merge=100)

最后，输出上述分析结果到文件：

write.table(dmrs, “cg_dmrs.bed”, sep=”\t”,

+ row.names=FALSE, quote=FALSE)

六、第6步 - 可视化

DMR注释、差异碱基位点、差异甲基化区域等信息均可进行可视化，以直观展示甲基化水平的模式。

常用的可视化工具有IGV、UCSC基因组浏览器、deeptools等。

可视化示例

参考资料：

https:///10.1038/nature06745

https:///10.14806/ej.17.1.200

https:///10.1038/nmeth.1923

https:///10.1186/1471-2164-14-774

https:///10.1093/bioinformatics/btp352

https:///10.1093/nar/gku154

https://github.com/xie186/ViewBS

https:///10.1093/nar/gku365

微信公众号（生信笔记）同步更新：https://mp.weixin.qq.com/s/Yw75_43h3afLs8UC4HlbKA