在《WGBS 全基因组甲基化测序》)讲解了 WGBS 测序原理，简而言之是 BS 处理将未甲基化 C 转化为 U(T) 而甲基化的 C 不变，从而将两者做出区分。
Bismark 是知名的甲基化测序比对软件。测序 DNA BS 转化后经过 PCR 扩增发生了 C->T 转化，那么反向互补链则是 G->A 转化。Bismark 事先无法知道 Reads 的方向性，因此将所有 Reads 进行 C->T 和 G->A 转换，并且将参考基因组同样进行 2 种转换。组合起来每条 Reads 进行 4 次不同比对，从里面选择最佳比对作为结果。这样方案可行主要因为甲基化 C 是少数，大多数是未甲基化的。从最佳比对结果反推可以知道序列的链方向，根据错配的 CT 可以知道哪些碱基发生了甲基化。
Bismark比对
建立索引，使用目录路径，在目录下放要建索引的 fasta 文件会被自动识别。Bismark 支持 Bowtie 2 和 Hisat2，默认会用 Bowtie 2 建议将 Bowtie 2 添加到 PATH 就不用每次指定路径。
bismark_genome_preparation --genomic_composition /Example/openData/GRCh38_hg38/Bismark
比对，默认链特异建库所以只进行 OT 和 OB 比对就行，如果明确不是链特异建库，需要设置 --non_directional 参数进行完整的 4 次比对。注意多线程参数 --parallel 的设置，本身 bismark 的不少步骤(Bowtie 2, samtools...)已经是多线程的了。如果服务器内存大CPU核心多可以设置这个参数加快比对。小鼠基因组设置为 4 时约消耗 20 核心 40GB 内存。bismark 比对首先要将 reads 进行 C->T 和 G->A 转换然后进行 2/4 次比对，所以不设置多线程时耗时确实很多，根据自己硬件资源尽量跑满。
bismark --parallel 4 --output_dir /Example/WGBS/Bismark_1 --gzip --nucleotide_coverage \
/Example/openData/GRCh38_hg38/Bismark \
-1 /Example/WGBS/CleanData/BJBC001T1_R1.fq.gz \
-2 /Example/WGBS/CleanData/BJBC001T1_R2.fq.gz
WGBS 需要进行去重。RRBS 的话千万别去重。
deduplicate_bismark --paired --outfile BJBC001T1 --output_dir /Example/WGBS/Bismark_1 /Example/WGBS/Bismark_1/BJBC001T1_R1_bismark_bt2_pe.bam
然后是关键的提取甲基化数据，可以分 2 次进行。第一次用 --mbias_only 不提取数据，只产生 M-bias 质控图需要的数据，环境配置对的会自动生成图像，我没配置好就自己画图了。
下面的 M-bias 图左边显示不同碱基位置每一种甲基化比例，因为建库打断是随机的理论上不同碱基位置比例相同，所以好的数据大致是每种类型一条横线。右边部分是每种类型 C 总数目。上面 Read1 下面 Read2 。其中 Read2 前几个碱基的甲基化比例非常低，应该是建库末端修复步骤引入的非甲基化 C 碱基，应该移除。
结合 Fastqc 检测报告可以看出无论 Read1 还是 Read2 前10碱基占比波动很大。所以我这里第二次运行 bismark_methylation_extractor 进行提取甲基化数据时选择忽略前 10 碱基数据。有些流程为了方便默认将 2 端忽略 3-4 碱基，想省事的可以学这个，但我这自己手动跑还是依据质控图选择过滤参数。
M-bias
Read1 Fastq
Read2 Fastq
参数 --ignore 和 --ignore_r2 分别忽略 2 Reads 5' 端碱基。如果还需要配置 3' 端忽略，用 --ignore_3prime 和 --ignore_3prime_r2 参数。
bismark_methylation_extractor --paired-end --ignore 10 --ignore_r2 10 --comprehensive --output /datapool/pengguoyu/20200421BLCA_Omics/WGBS/Bismark_1 --parallel 4 --bedGraph --cytosine_report --genome_folder /share/database/openData/GRCh38_hg38/Bismark \
/Example/WGBS/Bismark_1/BJBC001T1.deduplicated.bam \
/Example/WGBS/Bismark_1/BJBC005T1.deduplicated.bam \
/Example/WGBS/Bismark_1/BJBC007T1.deduplicated.bam \
/Example/WGBS/Bismark_1/BJBC008T1.deduplicated.bam
运行后得到每一种甲基化的具体碱基位置和状态，以 CpG 甲基化为例。
Bismark methylation extractor version v0.22.3
CL100137746L1C001R001_234/1 - chr12 101373135 z
CL100137746L1C001R001_234/1 + chr12 101373168 Z
CL100137746L1C001R001_817/1 + chr20 37841302 Z
CL100137746L1C001R001_91/1 + chr8 75699871 Z
CL100137746L1C001R001_817/1 + chr20 37841355 Z
往后在这些结果目录用 bismark2report 就能得到 HTML 报告。包含比对情况，甲基化情况，碱基比例等。
比对概况
甲基化概况
碱基概况
Bismark 比对基本到这里结束，后面是更下游的分析以及可视化等。
网友评论