用 spike-in DNA 计算 BS 处理的转化率

WGBS 实验 BS 处理 C->U(T) 转化效率不是 100% 的，转化失败的部分 C 会成为假阳性信号，所以 BS 转化率也是很重要的指标。计算转化率需要建库实验参入已知不含甲基化的 DNA ，理论上所有该 DNA 的 C 都会被转化为 U(T)，用 Bismark 软件分析会得到甲基化率为 0. 实际肯定不为 0，这部分甲基化就是转化失败的，所以 1 减去甲基化率就是转化率。
在这里以参入的 DNA 是 Enterobacteria phage lambda() 为例子演示如何计算 BS 处理转化率。参考基因组序列 NC_001416.1 可以在 NCBI 下载。

下载参考基因组 NC_001416.1 后先用 bismark_genome_preparation 建立索引。

bismark_genome_preparation --genomic_composition /Example/Enterobacteria_phage_lambda

然后 bismark 命令将 fastq 文件比对到 NC_001416.1 参考基因组。

bismark  --parallel 4 --output_dir /Example/Bismark_lambda --gzip --nucleotide_coverage /Example/Enterobacteria_phage_lambda -1 /Example/CleanData/${sample}_R1.fq.gz -2 /Example/CleanData/${sample}_R2.fq.gz

最后从比对报告里手动计算转化率。

Final Alignment report
======================
Sequence pairs analysed in total:       552427505
Number of paired-end alignments with a unique best hit: 2788190
Mapping efficiency:     0.5% 
Sequence pairs with no alignments under any condition:  549639315
Sequence pairs did not map uniquely:    0
Sequence pairs which were discarded because genomic sequence could not be extracted:    1737

Number of sequence pairs with unique best (first) alignment came from the bowtie output:
CT/GA/CT:       1443997 ((converted) top strand)
GA/CT/CT:       0       (complementary to (converted) top strand)
GA/CT/GA:       0       (complementary to (converted) bottom strand)
CT/GA/GA:       1342456 ((converted) bottom strand)

Number of alignments to (merely theoretical) complementary strands being rejected in total:     0

参入的 DNA 是非常少量的，所以比对率非常低(这里 0.5%)如果比对率高需要检查有没有出问题。

Final Cytosine Methylation Report
=================================
Total number of C's analysed:   121038850

Total methylated C's in CpG context:    99903
Total methylated C's in CHG context:    128558
Total methylated C's in CHH context:    299809
Total methylated C's in Unknown context:        154

Total unmethylated C's in CpG context:  27624300
Total unmethylated C's in CHG context:  29106450
Total unmethylated C's in CHH context:  63779830
Total unmethylated C's in Unknown context:      10902

C methylated in CpG context:    0.4%
C methylated in CHG context:    0.4%
C methylated in CHH context:    0.5%
C methylated in unknown context (CN or CHN):    1.4%

这里总的 C 碱基数目 121038850 跟下面甲基化和非甲基化的和是对不上的，如果不计算 Unknown context 就对的上，也许是概念有些重复。

> 99903 + 128558 + 299809 + 27624300 + 29106450 + 63779830
[1] 121038850

把甲基化的 C 相加除以总的 C 就是甲基化率(其实是假阳性)，用 1 减去这个甲基化率就是 BS 处理转化率。

> 99903 + 128558 + 299809
[1] 528270
> 528270 / 121038850
[1] 0.004364466
> 1 - 0.004364466
[1] 0.9956355

所以这里 BS 处理转化率为 99.57% 。

另一个方案是将参入 DNA 基因组与人基因组一起建立索引并正常进行 Bismark 流程，后面将属于参入 DNA 的甲基化输出数据提取计算甲基化率。未曾尝试，仅作介绍。