生信软件48 - 差异甲基化区域（DMR）识别的DNA 甲基化分析工具methylpy-优快云博客

本文链接：https://blog.youkuaiyun.com/LittleComputerRobot/article/details/149716507

methylpy 是一个用于分析亚硫酸氢盐测序（Bisulfite Sequencing, BS-seq）数据的 Python 库，专门针对哺乳动物基因组的 DNA 甲基化分析设计。它提供了从原始测序数据处理到差异甲基化区域（DMR）识别的完整流程，尤其适用于处理低输入量样本（如单细胞或胚胎）的甲基化数据。

可实现计算单个 CpG 位点或基因组区域的甲基化水平，识别低甲基化区域（LMRs）和高甲基化区域（HMRs）；差异甲基化分析（DMR calling），支持多样本比较；生成甲基化热图、箱线图和基因组浏览器视图，将甲基化区域注释到基因、启动子或功能元件等。

2. 安装

依赖要求cutadapt （>=1.9）用于原始读取修剪；picard （>=2.10.8）用于 PCR 重复去除；samtools （>=1.3）；wigToBigWig 用于格式转换。

####### conda安装（推荐）#######
conda env create --name methylpy
conda activate methylpy
conda install -y -c bioconda -c conda-forge methylpy    

wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig

####### github安装 #######
# clone项目
git clone https://github.com/yupenghe/methylpy.git
cd methylpy/
python setup.py install

帮助信息

3. 测试数据下载

wget http://neomorph.salk.edu/yupeng/share/methylpy_test.tar.gz
tar -xf methylpy_test.tar.gz
cd methylpy_test/
python run_test.py

4. 基本用法

4.1 基本分析流程

# 1. 构建转换后的基因组参考
methylpy build-reference \
	--input-files mm10_bt2/mm10.fa \
	--output-prefix mm10_bt2/mm10 \
	--bowtie2 True
	
	
# 2.1  处理亚硫酸氢盐测序和 NOMe-seq 数据
# SE数据
methylpy single-end-pipeline \
	--read-files raw/mESC_R1.fastq.gz \
	--sample mESC \
	--forward-ref mm10_bt2/mm10_f \
	--reverse-ref mm10_bt2/mm10_r \
	--ref-fasta mm10_bt2/mm10.fa \
	--num-procs 8 \
	--remove-clonal True \
	--path-to-picard="picard/"
	
# PE数据
methylpy paired-end-pipeline \
	--read1-files raw/mESC_R1.fastq.gz \
	--read2-files raw/mESC_R2.fastq.gz \
	--sample mESC \
	--forward-ref mm10_bt2/mm10_f \
	--reverse-ref mm10_bt2/mm10_r \
	--ref-fasta mm10_bt2/mm10.fa \
	--num-procs 8 \
	--remove-clonal True \
	--path-to-picard="picard/"
	
	
# 获取差异甲基化区域DMR	
# 接受压缩/未压缩的 allc 文件列表（来自 methylpy 管道的输出文件）作为输入并查找 DMR
methylpy DMRfind \
	--allc-files allc/allc_AD_HT.tsv.gz allc/allc_AD_IT.tsv.gz \
	--samples AD_HT AD_IT \
	--mc-type "CGN" \
	--chroms 1 2 3 4 5 \
	--num-procs 8 \
	--output-prefix DMR_HT_IT

输出文件：是 allc 格式的（压缩）制表符分隔文本文件。“allc”代表对于所有胞嘧啶（C）。allc 文件中的每一行对应于基因组中的一个胞嘧啶。一个 allc 文件包含 7 个必填列，没有标题。

4.2 从 BAM 文件中提取胞嘧啶甲基化状态

methylpy call-methylation-state \
	--input-file mESC_processed_reads_no_clonal.bam \
	--paired-end True \
	--sample mESC \
	--ref-fasta mm10_bt2/mm10.fa \
	--num-procs 8

4.3 获取基因组区域的甲基化水平

计算某些基因组区域的甲基化水平可以估计甲基化这些基因座的丰富性。

methylpy add-methylation-level \
	--input-tsv-file DMR_AD_IT.tsv \
	--output-file DMR_AD_IT_with_level.tsv \
	--allc-files allc/allc_AD_HT_1.tsv.gz allc/allc_AD_HT_2.tsv.gz \
		allc/allc_AD_IT_1.tsv.gz allc/allc_AD_IT_2.tsv.gz \
	--samples AD_HT_1 AD_HT_2 AD_IT_1 AD_IT_2 \
	--mc-type CGN \
	--num-procs 4

4.4 多个 allc 文件合并

将多个 allc 文件合并为一个 allc 文件。

methylpy merge-allc \
	--allc-files allc/allc_AD_HT_1.tsv.gz allc/allc_AD_HT_2.tsv.gz \
	--output-file allc/allc_AD_HT.tsv.gz \
	--num-procs 1 \
	--compress-output True

4.5 过滤 allc 文件

按胞嘧啶上下文、覆盖率等过滤位点。

methylpy filter-allc \
	--allc-file allc/allc_AD_HT_1.tsv.gz \
	--output-file allc/allCG_AD_HT_1.tsv.gz \
	--mc-type CGN \
	--min-cov 2 \
	--compress-output True

4.6 索引 allc 文件

允许为每个allc文件创建索引文件。索引文件可用于加快 allc 文件读取速度，类似于 .fasta 文件的 .fai 文件。

methylpy index-allc \
	--allc-files allc/allc_AD_HT_1.tsv.gz allc/allc_AD_HT_2.tsv.gz \
	--num-procs 2 \
	--no-reindex False

4.7 将allc文件转换为bigwig格式

从 allc 文件生成 bigwig 文件。甲基化水平将在平均分配的非重叠基因组箱中计算，输出将存储bigwig文件中。

methylpy allc-to-bigwig \
	--allc-file results/allc_mESC.tsv.gz \
	--output-file results/allc_mESC.bw \
	--ref-fasta mm10_bt2/mm10.fa \
	--mc-type CGN \
	--bin-size 100

4.8 用于亚硫酸氢盐测序读数的优质过滤器

过滤掉无法比对上或可能转化不足的 DNA 片段。

# 以下命令可用于过滤掉 MAPQ 分数小于 30 的读取 （对齐不良）和 mCH 水平大于 0.7（转换不足），如果读数包含 足够（至少 3 个）CH 位点
methylpy bam-quality-filter \
	--input-file mESC_processed_reads_no_clonal.bam \
	--output-file mESC_processed_reads_no_clonal.filtered.bam \
	--ref-fasta mm10_bt2/mm10.fa \
	--min-mapq 30 \
	--min-num-ch 3 \
	--max-mch-level 0.7 \
	--buffer-line-number 100

4.9 从现有结果中重新识别 DMR

根据之前 DMRfind 运行的结果重新识别 DMR。。

methylpy reidentify-DMR \
	--input-rms-file results/DMR_P0_FBvsHT_rms_results.tsv.gz \
	--output-file results/DMR_P0_FBvsHT_rms_results_recollapsed.tsv \
	--collapse-samples P0_FB_1 P0_FB_2 P0_HT_1 P0_HT_2 \
	--sample-category P0_FB P0_FB P0_HT P0_HT \
	--min-cluster 2