基于基础Python实现多样本SNP密度滑动窗口绘图的问题咨询

滑动窗口步长核心逻辑实现 #

你已经完成了VCF读取、基因型清洗、突变位点提取，首先将输出整理为如下结构的标准数据：
sample_snp_pos = {"样本1": {"chr1": [100, 230, 450, ...], "chr2": [120, 360, ...]}, "样本2": ...}
其中每个染色体对应的SNP坐标列表需按从小到大排序（VCF输出默认已排序，可直接使用）。

使用Python标准库bisect实现高效的滑动窗口计算，无需引入任何生物相关第三方库，步长逻辑直接通过窗口起始位置累加步长值实现：

import bisect
from collections import defaultdict

def calc_snp_density(sample_snp_pos, window_size, step_size, chrom_lengths=None):
    """
    计算各样本各染色体的滑动窗口SNP密度
    参数：
        sample_snp_pos: 整理好的样本突变位点坐标
        window_size: 滑动窗口大小（单位bp）
        step_size: 步长（单位bp）
        chrom_lengths: 可选，各染色体长度字典，无传入时默认取染色体最大SNP位置为计算终点
    返回：
        {样本名: {染色体: [(窗口起始, 窗口结束, SNP密度), ...]}}
    """
    res = defaultdict(lambda: defaultdict(list))
    for sample, chrom_map in sample_snp_pos.items():
        for chrom, pos_list in chrom_map.items():
            if not pos_list:
                continue
            # 确定染色体计算终点
            max_pos = chrom_lengths[chrom] if chrom_lengths else pos_list[-1]
            current_win_start = 1
            # 滑动窗口遍历
            while current_win_start <= max_pos:
                current_win_end = current_win_start + window_size - 1
                # bisect快速统计区间内SNP数量，时间复杂度O(logN)
                left = bisect.bisect_left(pos_list, current_win_start)
                right = bisect.bisect_right(pos_list, current_win_end)
                snp_count = right - left
                density = snp_count / window_size
                res[sample][chrom].append((current_win_start, current_win_end, density))
                # 步长移动窗口：核心逻辑就是起始位置累加步长值
                current_win_start += step_size
    return res

步长参数控制规则：

• 步长=窗口大小：窗口无重叠，适合全基因组全局密度统计
• 步长<窗口大小：相邻窗口存在重叠区域，适合高密度精细分析
• 步长>窗口大小：相邻窗口存在间隔，适合低覆盖度数据的快速统计

多样本折线图生成 #

使用通用可视化库matplotlib实现多样本折线图输出，不属于你禁止的生物相关工具范畴：

import matplotlib.pyplot as plt

def draw_density_line(density_res, target_chrom, window_size, step_size, save_path="snp_density.png"):
    plt.figure(figsize=(12, 6))
    for sample, chrom_data in density_res.items():
        if target_chrom not in chrom_data:
            continue
        win_list = chrom_data[target_chrom]
        # 取窗口中点为X轴坐标，折线展示更连贯
        x_axis = [(win[0] + win[1])/2 for win in win_list]
        y_axis = [win[2] for win in win_list]
        plt.plot(x_axis, y_axis, label=sample, linewidth=1)
    plt.xlabel(f"基因组位置（{target_chrom}，单位bp）")
    plt.ylabel("SNP密度（单位：个/bp）")
    plt.title(f"滑动窗口SNP密度分布（窗口大小{window_size}bp，步长{step_size}bp）")
    plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
    plt.tight_layout()
    plt.savefig(save_path, dpi=300)
    plt.close()

如果你不希望使用任何第三方库，可以将calc_snp_density输出的结果直接写入csv文件，再用系统自带的表格工具生成折线图即可。

注意事项 #

• 滑动窗口计算默认按单条染色体独立进行，不跨染色体统计，符合群体遗传分析的通用规则
• 如果需要生成全基因组拼接的密度图，可对不同染色体的X轴坐标做固定偏移后再拼接绘图
• VCF中的INDEL位点如果不需要统计，你之前的清洗步骤过滤即可，不影响后续计算逻辑

内容的提问来源于stack exchange，提问作者user8769986