基于LD的动态窗口法

基于LD的动态窗口法(代码在附件中)原理通过计算SNP间的连锁不平衡系数r²根据实际遗传相关性动态调整block边界。重要特性允许孤立SNP不与其他SNP连锁的SNP不被强行加入任何block更符合生物学实际。算法流程初始化窗口大小500kbLD阈值r² 0.5对每条染色体当前位置 染色体起始 while 当前位置 染色体终止: 获取500kb窗口内的SNP if 窗口内SNP 2个: 跳过这些SNP不形成block continue 采样SNP最多100个避免计算量爆炸 计算窗口内SNP的平均r² if 平均r² 0.15: # 非常低的LD 跳过这些SNP孤立SNP不形成block continue elif 平均r² ≥ 0.5: # 高LD区域 扩展block至窗口末端 elif 平均r² 0.25: # 低LD区域 缩小block至窗口中点 else: # 中等LD 保持窗口大小 if block包含SNP数 ≥ 2: 形成有效block else: 跳过不形成block 移动到下一个位置LD r²计算r² [相关系数(SNP1, SNP2)]² 其中基因型编码 - 0 A2A2主等位基因纯合子 - 1 A1A2杂合子 - 2 A1A1次等位基因纯合子优化策略孤立SNP过滤新增判断标准窗口内平均r² 0.15阈值的30%最小block大小至少2个SNP生物学意义真正独立遗传的SNP不被强行分组示例某SNP与周围所有SNP的r²都0.15则不加入任何block采样加速每个窗口最多采样100个SNP均匀采样保证代表性计算复杂度O(n)而非O(n²)r²估计小窗口≤10 SNP计算所有配对大窗口10 SNP随机采样50对进度显示染色体级别进度条SNP读取进度条样本生成进度条 输出结果运行时统计信息划分了 2845 个LD blocks 包含 13,456,789 个SNP排除了 746,285 个孤立SNP说明包含的SNP存在于某个block中与其他SNP有LD关系孤立SNP不与周围任何SNP显著连锁r² 0.15独立遗传ld_blocks_info.csv包含每个block的详细信息列名说明示例block_idblock编号0, 1, 2…chr染色体1, 2, 3…start_pos起始位置(bp)1234567end_pos终止位置(bp)1734567length长度(bp)500000n_snpsSNP数量1523mean_ld_r²平均LD系数0.6234样本CSV文件每个样本一个文件如04K5672_ld_blocks.csvchrsplengthvarrankvartypeallele_freq119214983290.762300.854215002504998750.689110.3214 参数调整建议LD_WINDOW_KB窗口大小50-200kb精细划分适合高密度芯片500kb默认平衡精度和速度1000kb粗略划分加快计算LD_R2_THRESHOLDLD阈值0.2-0.3宽松标准更大block0.5默认中等标准0.7-0.8严格标准更小block⚡ 性能说明以477个样本、1420万SNP为例阶段耗时估计读取数据10秒划分blocksLD计算30-60分钟批量读取基因型5-10分钟生成477个CSV2-5分钟总计40-75分钟注意首次运行划分blocks较慢需计算LD但结果会保存在ld_blocks_info.csv中。 生物学意义什么是LD连锁不平衡两个SNP的等位基因在群体中非随机组合r² 1完全连锁r² 0完全独立r² 0.5强连锁为什么要基于LD划分生物学合理block内SNP趋向于共同遗传降维有效高LD的SNP携带冗余信息单倍型块符合人类基因组的block结构过滤孤立SNP新增不强制分组真正独立的SNP保持分析纯粹性 引用如果使用此方法建议引用Gabriel et al. (2002) The Structure of Haplotype BlocksPearson correlation-based LD calculationLD Block分析工具使用说明快速开始1. 安装依赖pipinstallnumpy pandas tqdm2. 运行脚本python generate_ld_block_csv.py3. 修改配置适配其他文件打开generate_ld_block_csv.py修改主函数中的配置# 配置区域 FILE_PREFIX477-merged# 改成你的文件前缀OUTPUT_DIRld_block_output# 输出目录LD_WINDOW_KB500# LD窗口大小kbLD_R2_THRESHOLD0.5# LD r²阈值0-1之间# LD Block划分方法基于真实LD关系的智能划分脚本使用LD-aware动态窗口法考虑实际连锁不平衡结构初始窗口: 500kb滑动窗口LD计算: 计算窗口内SNP间的平均r²值动态调整:高LD区域r² ≥ 0.5→ 扩展block边界低LD区域r² 0.25→ 缩小block边界中等LD区域 → 保持窗口大小参数说明:LD_WINDOW_KB 500: 基础窗口500kbLD_R2_THRESHOLD 0.5: r²高于0.5认为存在强LD与固定窗口法的区别方法优点缺点固定窗口快速简单忽略真实LD结构LD-aware符合生物学意义计算较慢输出文件格式每个样本生成一个CSV文件格式如下chr,sp,length,varrank,vartype,allele_freq 1,1921,998329,0.7623,0,0.8542 1,1000250,999875,0.6891,1,0.3214 2,50123,1000000,0.5432,2,0.5123 ...列说明chr: 染色体编号1, 2, 3, …sp: block起始位置碱基对位置length: block长度bpvarrank: 序列均匀度0-1之间越高表示基因型越多样化vartype: 该样本在此block的主要基因型0: 主等位基因纯合子A2A21: 次等位基因纯合子A1A12: 杂合子A1A2或其他allele_freq: block中主要等位基因A2的频率0-1之间文件要求必须包含以下3个文件文件前缀相同xxx.fam- 样本信息xxx.bim- SNP位点信息xxx.bed.gz或xxx.bed- 基因型数据性能说明477个样本约2000万SNP位点预计处理时间10-30分钟取决于CPU和磁盘速度内存占用约2-4GB常见问题Q: 如果BED文件未压缩怎么办A: 脚本会自动检测优先使用.bed.gz如果不存在则使用.bedQ: 如何更改block划分策略A: 修改LDBlockAnalyzer.define_blocks_simple()方法或实现基于LD r²的划分算法Q: 输出文件太多怎么办A: 可以修改脚本只处理特定样本# 只处理前10个样本forsample_idxinrange(10):analyzer.generate_sample_csv(sample_idx,output_dir)