2022-05-25甲基化转化率统计

Bisulfitesequencing，实验 BS 处理 C->U(T) 转化效率并不是 100% 的，转化失败的部分 C 会成为假阳性信号，所以 BS 转化率也是很重要的指标。计算转化率需要建库实验掺入已知不含甲基化的 DNA（由于线粒体 DNA 可以认为是无甲基化修饰的 DNA，也可使用线粒体 DNA 来评估），理论上所有该 DNA 的C 都会被转化为 U(T)，通过比对和甲基化计算分析，会得到该 DNA 的甲基化率为 0. 实际肯定不为 0，这部分甲基化就是转化失败的，所以 1 减去甲基化率就是转化率。BS 转化率，我们一般要求要在 99%以上。掺入已知不含甲基化的 DNA 一般常用是 Enterobacteria phage lambda (lambda DNA)。参考基因组序列 NC_001416.1 可以在 NCBI 下载：https://www.ncbi.nlm.nih.gov/nuccore/NC_001416.1如果使用样本本身的线粒体 DNA 来评估，建议使用线粒体 DNA 上非 CG 位点/或者 CC 位点的甲基化水平来评估。

1.比对

bsmap -a 1.fq -b 2.fq -d lambda.fa/chrMt.fa -o BSconvert.sam -p 8 -z 33 -r 1 -m 50 -x 1000

2.sam2bam

samtools view -Sb BSconvert.sam > BSconvert.bam

3.提取甲基化位点

python methratio.py -d lambda.fa/chrMt.fa -m 1 -o BSconvert.met.txt BSconvert.bam

4.统计甲基化转化率

perl -alne 'next if $.==1; next if $F[3] eq "CG"; $m+=$F[6]; $n+=$F[7]; END{ print "“BS_convert_ratio\t".sprintf("%.5f", $n/($m+$n));}' BSconvert.met.txt