低覆盖度无填充全基因组测序数据：ngsLD Perl修剪脚本替代方案问询

针对你180例低覆盖度全基因组样本的LD修剪需求，结合HPC7天任务时限的约束，以下是几个适配低覆盖数据且速度远优于Perl脚本的方案：

1. PLINK 2.0（首选） #

PLINK 2.0是C语言重写的版本，计算效率比Perl脚本高一个数量级，同时针对低覆盖数据做了优化：

• 核心优势：支持基于似然的LD计算（而非硬基因型计数），更适配低覆盖下的基因型不确定性；原生支持多线程，可充分利用HPC节点资源；支持灵活的区域拆分并行。
• 操作示例：

  # 按染色体拆分+多线程运行，同时做预过滤适配低覆盖
  plink2 --bfile chr${chr}_input \
         --indep-pairwise 50 5 0.2 \  # 滑动窗口LD修剪参数
         --geno 0.2 \  # 保留缺失率≤20%的SNP（适配低覆盖）
         --maf 0.05 \  # 过滤低频SNP减少计算量
         --out chr${chr}_ld_pruned \
         --threads 32  # 利用节点全部核心
  

• 额外优化：将染色体拆分为10Mbp左右的子区域并行提交任务，单个任务运行时间可控制在7天内。

2. GCTA #

GCTA的LD修剪模块同样适配低覆盖数据，且计算速度快，支持多线程：

• 核心优势：可通过--impute-avg参数用等位基因频率填充低覆盖导致的缺失基因型，提升LD计算的可靠性；支持按染色体或区域并行。
• 操作示例：

  gcta64 --bfile chr${chr}_input \
         --ld-prune 0.2 \  # LD阈值0.2
         --ld-window 50 \  # 窗口大小50个SNP
         --ld-window-step 5 \  # 滑动步长5个SNP
         --impute-avg \  # 填充缺失基因型
         --out chr${chr}_ld_pruned \
         --thread-num 32
  

3. SNPrelate（R包） #

适合熟悉R环境的用户，支持基于基因型剂量的LD计算，完美适配低覆盖数据：

• 核心优势：直接处理VCF/PLINK格式的剂量数据（低覆盖样本常用输出格式），LD计算更准确；支持多线程并行。
• 操作示例：

  library(SNPrelate)
  # 读取PLINK文件
  genofile <- snpgdsOpen("input.bed")
  # LD修剪，设置多线程
  pruned <- snpgdsLDpruning(genofile, 
                            method="r2", 
                            ld.threshold=0.2,
                            slide.max.bp=50000,
                            ncores=32)
  # 提取修剪后的SNP列表
  pruned_snps <- unlist(pruned)
  snpgdsClose(genofile)
  

通用优化策略 #

• 预过滤升级：先过滤缺失率>20%、MAF<0.03的SNP，减少后续计算量；
• 任务拆分细化：将染色体拆分为5-10Mbp的子区域，每个子区域作为独立任务提交，确保单任务运行时长不超过7天；
• 节点资源申请：优先申请高核心数（32/64核）、大内存的HPC节点，最大化并行效率。

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