前言
因?yàn)榻M里面出了一批甲基化測(cè)序數(shù)據(jù)，使用的技術(shù)為BS-seq，處理的時(shí)候順帶記錄了學(xué)習(xí)過程，演示使用數(shù)據(jù)為官方提供的example.fastq。
DNA甲基化及CpG島定義了解
DNA甲基化（英語：DNA methylation）為DNA化學(xué)修飾的一種形式，能在不改變DNA序列的前提下，改變遺傳表現(xiàn)。為外遺傳編碼（epigenetic code）的一部分，是一種外遺傳機(jī)制。DNA甲基化過程會(huì)使甲基添加到DNA分子上，例如在胞嘧啶環(huán)的5'碳上：這種5'方向的DNA甲基化方式可見于所有脊椎動(dòng)物。
在人類細(xì)胞內(nèi)，大約有1%的DNA堿基受到了甲基化。在成熟體細(xì)胞組織中，DNA甲基化一般發(fā)生于CpG雙核苷酸（CpG dinucleotide）部位；而非CpG甲基化則于胚胎干細(xì)胞中較為常見[1] [2]。植物體內(nèi)胞嘧啶的甲基化則可分為對(duì)稱的CpG（或CpNpG），或是不對(duì)稱的CpNpNp形式（C與G是堿基；p是磷酸根；N指的是任意的核苷酸）。
特定胞嘧碇受甲基化的情形，可利用亞硫酸鹽定序（bisulfite sequencing）方式測(cè)定。DNA甲基化可能使基因沉默化，進(jìn)而使其失去功能。此外，也有一些生物體內(nèi)不存在DNA甲基化作用。
——維基百科
DNA甲基化作為基因組上的表觀修飾（區(qū)別于組蛋白修飾），存在于各種生物中。
簡(jiǎn)單來說，可以認(rèn)為甲基化代表著基因的失活，去甲基化則標(biāo)志基因的激活與表達(dá)
雖然CpG序列出現(xiàn)的頻率并不高，但是在某些基因區(qū)域內(nèi)，CpG的密度很高，俗稱CpG島。這些CpG島大多出現(xiàn)在基因的啟動(dòng)子區(qū)域（人類占到70%），長(zhǎng)度達(dá)300-3000bp。目前的研究表明，大多數(shù)的管家基因都含有CpG島，位于基因的5'端（其中的大多數(shù)CpG島都是未甲基化的）。
另外需要注意的是，目前的研究表明，腫瘤樣本與正常樣本的CpG島甲基化差異大多不是發(fā)生CpG島的內(nèi)部而是位于CpG島岸（CpG island shore）。
由于CpG位點(diǎn)的易甲基化導(dǎo)致胞嘧啶脫氨變成胸腺嘧啶，所以在漫長(zhǎng)的進(jìn)化過程中，CpG位點(diǎn)逐漸消失，但是又存在著對(duì)于基因表達(dá)的調(diào)控要求，所以CpG島的出現(xiàn)也被理解為抵抗甲基化經(jīng)常很，維持調(diào)控功能。
DNA甲基化測(cè)序原理與方法
此處略過，請(qǐng)自行了解(示例文件為WGBS單端測(cè)序文件)。
其實(shí)是因?yàn)槔斫馄饋肀容^累，知道大致原理就可以了，以后再來填坑
數(shù)據(jù)處理（使用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頁面需要搭梯子才能看到，簡(jiǎn)書不支持附件，所以有需要的可以留言，我email pdf文件
Dependencies
需要用戶已經(jīng)裝好bowtie1/bowtie2
測(cè)試數(shù)據(jù)獲取
此處使用測(cè)試數(shù)據(jù)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文件構(gòu)建
# under install folder
./bismark_genome_preparation --path_to_bowtie /PATH/to/bowtie2/ --verbose ~/PATH/to/GRCh38/
比對(duì)
# 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
甲基化位點(diǎn)提取
./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
生成處理報(bào)告
./bismark2report
所有結(jié)果文件
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
結(jié)果解讀
--cytosine_report參數(shù)會(huì)根據(jù)當(dāng)前目錄下的信息文件生成一個(gè)HTML格式的報(bào)告文件，即test_data_bismark_bt2_SE_report.html文件，它包括了比對(duì)信息，甲基化信息，M-bias等，可以對(duì)數(shù)據(jù)有一個(gè)大概的認(rèn)知（下圖只展示了一部分）：
同時(shí)因?yàn)槭褂昧?code>--comprehensive，所以結(jié)果合并正反鏈的數(shù)據(jù)后會(huì)輸出CpG/CHG/CHH三種類型的甲基化文件，包含了胞嘧啶所有的組合形式，但實(shí)際上我們自然最關(guān)注的是CpG位點(diǎn)的甲基化。其中
CpG_context_test_data_bismark_bt2.deduplicated.txt即CpG甲基化位點(diǎn)的文件。
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
# 第一列為測(cè)序信息
# 第二列為甲基化狀態(tài) + 代表甲基化 -代表未甲基化
# 第三列代表chromosome
# 第四列代表location
# 第五列代表methylation call，簡(jiǎn)單來說大寫的就是甲基化的（因?yàn)檫€有CHG，CHH的數(shù)據(jù)，分別對(duì)應(yīng)x, X , h, H）
test_data_bismark_bt2.deduplicated.bismark.cov文件則給了每個(gè)位點(diǎn)的甲基化比例，為下一步確定CpG島提供了基礎(chǔ)，其數(shù)據(jù)形式如下：
$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
# 第四列代表甲基化百分比
# 第五列代表甲基化數(shù)目
# 第六列代表未甲基化數(shù)目
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
# 第四，五列為甲基化和非甲基化的堿基數(shù)目
# 第六列為CG
# 第七列為背景（第一個(gè)C延伸兩個(gè)堿基）
其它參數(shù)
bam2nuc
bismark2summary
coverage2cytosine
NOMe_filtering
filter_non_conversion
# 有需要的可以自信探索 --help or manual
結(jié)語
此處根據(jù)測(cè)序數(shù)據(jù)得到了甲基化位點(diǎn)的信息，但是后續(xù)DML以及DMR的確定還需要R包的使用，以及后續(xù)的可視化還以探索以下包：
Methy-Pipe: An Integrated Bioinformatics Pipeline for Whole Genome Bisulfite Sequencing Data Analysis
2014年出的一個(gè)pipeline，分析作圖一條龍，有興趣的同學(xué)安排一下。
