参考来源1：https://www.jianshu.com/p/5d7e550abc1a
参考来源2：https://www.jianshu.com/p/971ff2ee533a
参考来源3：https://www.plob.org/article/10843.html
第二部分 实战篇之数据比对
1. 数据处理（使用Bismark软件处理）
1.1 软件安装
- Bismark
- Bowtie2
1.2 测序数据下载
推荐使用aspera软件。
1.3 SRA文件处理
fastq-dump --gzip --split-3 SRR3503136.sra
fastq-dump --gzip --split-3 SRR3503134.sra
fastq-dump --gzip --split-3 SRR3503142.sra
1.4 准备基因组（构建基因组索引文件）
将rice基因组文件放在index目录下，后缀为.fa或者.fasta，然后运行以下命令:
# under install folder
./bismark_genome_preparation --verbose --bowtie2 ./index/
-
–bowtie2 选定bowtie2，如果bowtie不在环境变量路径内，还要另外指明
-
./index/ 选定参考基因组位置
运行成功后，会在index所在目录生成如下的目录结构：
可以看到，分别对CT和GA转换的基因组建立了bowtie2的索引。
1.5 测序文件比对
# paired-end data
./bismark --bowtie2 --genome ./index/ -1 ./data/SRR3503136_1.fastq.gz -2 ./data/SRR3503136_2.fastq.gz -p 8 -o ./Cr-DJ 2>Cr-DJ_mapping.log
./bismark --bowtie2 --genome ./index/ -1 ./data/SRR3503134_1.fastq.gz -2 ./data/SRR3503134_2.fastq.gz -p 8 -o ./osddm1a1b 2>osddm1a1b_mapping.log
./bismark --bowtie2 --genome ./index/ -1 ./data/SRR3503142_1.fastq.gz -2 ./data/SRR3503142_2.fastq.gz -p 8 -o ./osdrm2 2>osdrm2_mapping.log
# single-end data
./bismark --genome ~/PATH/to/index/ test_data.fastq -p 8 -o ./ 2>test.log
1.6 甲基化位点提取
./bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder ./index ./Cr-DJ/SRR3503136_1_bismark_bt2_pe.bam -o ./Cr-DJ_extractor 2>Cr-DJ_extractor.log
./bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder ./index ./osddm1a1b/SRR3503134_1_bismark_bt2_pe.bam -o ./osddm1a1b 2>osddm1a1b_extractor.log
./bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder ./index ./osdrm2/SRR3503142_1_bismark_bt2_pe.bam -o ./osdrm2_extractor 2>osdrm2_extractor.log
1.7 甲基化胞嘧啶提取
./coverage2cytosine --dir ./Cr-DJ -o Cr-DJ-cytosine_genome --genome_folder ~/Bismark_v0.20.0/index --CX ./Cr-DJ/SRR3503136_1_bismark_bt2_pe.bismark.cov.gz
./coverage2cytosine --dir ./osddm1a1b -o osddm1a1b-cytosine_genome --genome_folder ~/Bismark_v0.20.0/index --CX ./osddm1a1b/SRR3503134_1_bismark_bt2_pe.bismark.cov.gz
./coverage2cytosine --dir ./osdrm2 -o osdrm2-cytosine_genome --genome_folder ~/Bismark_v0.20.0/index --CX ./osdrm2/SRR3503142_1_bismark_bt2_pe.bismark.cov.gz
**提取出MC后，如需在IGV上查看reads在基因组上的分布情况，则需要将其进一步转换为wig文件：
perl make.wig.files.pl Cr-DJ-cytosine_genome.CX_report.txt Cr-DJ Cr-DJ 0.05
perl make.wig.files.pl osddm1a1b-cytosine_genome.CX_report.txt osddm1a1b osddm1a1b 0.05
perl make.wig.files.pl osdrm2-cytosine_genome.CX_report.txt osdrm2 osdrm2 0.05
# 过滤值可考虑用0.05，以输出更多位点。
From wig to TDF file: igvtools toTDF *.wig *.wig.tdf 5 90%
1.8 生成处理报告
./bismark2report
需要用root权限在上一步统计信息的文件夹下面运行，生成HTML格式的图形统计结果。
1.9 所有结果文件
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
1.10 结果解读（以软件自带的test数据为例）
–cytosine_report参数会根据当前目录下的信息文件生成一个HTML格式的报告文件，即test_data_bismark_bt2_SE_report.html文件，它包括了比对信息，甲基化信息，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
NOMe_filtering
filter_non_conversion
# 有需要的可以自信探索 --help or manual
2. 结语
此处根据测序数据得到了甲基化位点的信息，但是后续DML以及DMR的确定还需要R包的使用，以及后续的可视化还以探索以下包：
