如何使用MethylDackel：BS-seq实验的终极甲基化提取工具

MethylDackel是一款专为BS-seq（bisulfite sequencing）实验设计的强大且通用的甲基化分析软件，能够从坐标排序并带有索引的BAM或CRAM文件中提取每碱基的甲基化指标，为表观遗传学研究提供精准数据支持。

项目核心功能与优势

MethylDackel前身为PileOMeth，采用MIT许可协议发布，以其灵活性和广泛适用性在生物信息学领域备受青睐。该工具默认计算CpG位点的甲基化指标，同时支持CHG和CHH上下文中的甲基化分析，满足不同实验需求。

在处理数据时，MethylDackel会将所有胞嘧啶分为CpG、CHG和CHH三种序列上下文。其中，H代表除G以外的任何核苷酸，若参考序列中遇到N，则会根据情况归类为CHG或CHH上下文。对于染色体/重叠群末端的胞嘧啶，因无法推断上下文，默认归类为CHH上下文。

快速安装指南

Conda一键安装

确保系统中已安装Anaconda，通过以下命令可快速安装MethylDackel：

conda install -c bioconda methyldackel

源码编译安装

若需从源码编译，可按以下步骤操作：

git clone https://gitcode.com/gh_mirrors/me/MethylDackel
cd MethylDackel
make LIBBIGWIG="/some/path/to/libBigWig.a"
make install prefix=/some/installation/path

若链接器无法找到htslib头文件和库，可使用CFLAGS和LIBS指定路径：

make install CFLAGS="-O3 -Wall -I/some/path/include " LIBS="-L/some/path/lib" prefix=/some/installation/path LIBBIGWIG="/some/path/to/libBigWig.a"

基础使用方法

基本命令格式

MethylDackel最基本的使用方式如下：

MethylDackel extract reference_genome.fa alignments.bam

该命令会计算CpG位点的甲基化指标，并输出到alignments_CpG.bedGraph文件。文件中第4列表示甲基化C的读数/读对数量，第5列表示未甲基化C的对应数量。使用-o选项可指定输出文件名前缀。

高级选项设置

• 上下文选择：默认仅计算CpG上下文的甲基化指标，使用--CHH和--CHG选项可分别分析CHH和CHG上下文，--noCpG选项则忽略CpG中的胞嘧啶。
• 过滤参数：可通过-q和-p选项调整MAPQ和Phred分数的最低阈值，默认分别为10和5。
• 甲基化偏差校正：使用--OT、--OB、--CTOT和--CTOB选项可考虑甲基化偏差。

methylation bias分析与校正

在理想实验中，读段任意位置观察到甲基化C的概率应恒定，但实际中读段末端常出现甲基化率的增减，即甲基化偏差。MethylDackel提供了 methylation bias绘图功能，帮助用户识别和校正偏差：

MethylDackel mbias reference_genome.fa alignments.sorted.bam output_prefix

该命令会为每个有有效比对的链创建methylation bias图，输出SVG格式文件。下图展示了一个示例的methylation bias图，图中线条表示特定位置的平均甲基化百分比，阴影区域为99.9%置信区间：

从图中可观察到读段2末端甲基化的峰值和读段1起始处的谷值，这些区域可通过建议的修剪边界忽略，以提高数据准确性。

输出文件解析

bedGraph文件格式

MethylDackel extract生成的bedGraph文件包含6列制表符分隔的数据：

1. 染色体/重叠群/支架名称
2. 起始坐标
3. 结束坐标
4. 四舍五入为整数的甲基化百分比
5. 报告甲基化碱基的比对数/对数
6. 报告未甲基化碱基的比对数/对数

所有坐标均为0-based半开放格式，符合bedGraph定义。以下是输出示例：

track type="bedGraph" description="SRR1182519.sorted CpG methylation levels"
1	25115	25116	100	3	0
1	29336	29337	50	1	1

合并上下文输出

使用--mergeContext选项可将单个CpG或CHG中的胞嘧啶指标合并，生成每CpG或每CHG的指标。例如，将：

track type="bedGraph" description="SRR1182519.sorted CpG methylation levels"
1	25114	25115	66	2	1
1	25115	25116	100	3	0

转换为：

track type="bedGraph" description="SRR1182519.sorted merged CpG methylation levels"
1	25114	25116	83	5	1

若已存在每胞嘧啶指标的bedGraph文件，可使用MethylDackel mergeContext命令转换为每CpG/CHG指标。

高级功能与最佳实践

排除低覆盖区域

使用--minDepth选项可设置最小覆盖度阈值，默认值为1。例如，--minDepth 10表示仅输出覆盖度至少为10x的位点。该选项可与--mergeContext结合使用，对合并后的CpG或CHG进行过滤。

排除部分转化的读段

虽然一般不建议，但可使用--minConversionEfficiency选项过滤具有不完全亚硫酸氢盐转化证据的读段。该选项取值范围为0到1.0，默认值为0，表示忽略CHH和CHG上下文胞嘧啶的转化状态。

排除可能的变异位点

对于遗传异质样本，可使用--maxVariantFrac选项排除可能的变异位点。该选项通过跟踪参考序列中C的互补链上非G的数量，若其比例超过阈值，则排除该位置。--minOppositeDepth选项可设置互补链的最小覆盖深度，默认值为0，表示跳过变异位点排除过程。

典型应用案例

全基因组甲基化分析

对于大规模BS-seq数据分析，MethylDackel可并行处理多个染色体或分批次处理大文件，优化资源利用。建议预先进行MAPQ和Phred质量过滤，提高分析结果可信度。

与其他工具协同工作

在生物信息学生态系统中，MethylDackel常与BWA、samtools等数据预处理工具，以及minfi、ChAMP等R包协同工作。用户还可结合Bedops、BEDTools等工具进行区域操作，或结合GATK进行变异检测，进一步拓展其在基因组研究中的应用范围。

通过遵循上述指南，您可以充分利用MethylDackel的强大功能，为BS-seq实验的甲基化分析提供高效、准确的支持，推动表观遗传学研究的深入开展。