我最近基于nanopore的数据进行全基因组5mC甲基化分析,具体是基于王顺老师编写的脚本进行分析(https://gitee.com/wangshun1121/nanopolish-methyl),最后得到了频率文件,我想把它转换成划分窗口的全基因组甲基化水平文件,以便进行后续的可视化探索。因此,我写了一个python脚本,具体代码如下:
import pandas as pd
# 读取结果文件
df = pd.read_csv("results.meth_freq.tsv.gz", sep="\t", compression="gzip")
# 定义窗口大小
window_size = 40000
# 结果列表
results = []
# 按染色体和窗口分组计算甲基化水平
for chrom in df["chromosome"].unique():
chrom_data = df[df["chromosome"] == chrom]
# 获取染色体范围并按窗口分段
max_pos = chrom_data["end"].max()
for start in range(1, max_pos, window_size):
end = start + window_size - 1
# 筛选窗口内的数据
window_data = chrom_data[(chrom_data["start"] >= start) & (chrom_data["end"] <= end)]
# 计算平均甲基化水平
if not window_data.empty:
mean_methylation = window_data["methylated_frequency"].mean()
else:
mean_methylation = None
# 添加结果
results.append([chrom, start, end, mean_methylation])
# 转换为 DataFrame 并保存
output_df = pd.DataFrame(results, columns=["Chromosome", "Start", "End", "Methylation_Level"])
output_df.to_csv("genome_methylation_levels.tsv", sep="\t", index=False)
参考资料:https://www.jianshu.com/p/60d620a166fa