ld

Posted on 2024-08-04 Views: Waline: Word count in article: 6.2k Reading time ≈ 11 mins.

使用PopLDdecay进行ld分析

PopLDdecay是华大开发一款快速进行连锁不平衡衰减分析的工具，只需输入vcf文件就可完成所有分析，计算速度非常快。

连锁不平衡（LD）：当两个或多个位点的等位基因频率组合偏离随机组合时，这些位点之间被称为连锁不平衡。LD 的度量通常使用 r² 或 D’ 值。r² 值介于 0 和 1 之间，表示连锁不平衡的强度；D’ 值也介于 0 和 1 之间，表示连锁不平衡的比例。

计算

等位基因频率：

pAp_ApA 和 pBp_BpB 分别表示等位基因 A 和 B 在各自位点的频率。
pap_apa 和 pbp_bpb 分别表示等位基因 a 和 b 在各自位点的频率。

单倍型频率：

fABf_{AB}fAB：同时携带 A 和 B 的单倍型频率。
fAbf_{Ab}fAb：同时携带 A 和 b 的单倍型频率。
faBf_{aB}faB：同时携带 a 和 B 的单倍型频率。
fabf_{ab}fab：同时携带 a 和 b 的单倍型频率。

1. 计算 D 值

D 值表示两个等位基因的实际联合频率与预期联合频率之间的差异：

1	D=fAB−pA⋅pBD = f_{AB} - p_A \cdot p_BD=fAB−pA⋅pB

2. 计算 D’ 值

D’ 是对 D 进行标准化以去除等位基因频率的影响，其计算公式为：

1	D′=DDmaxD' = \frac{D}{D_{\text{max}}}D′=DmaxD

其中 DmaxD_{\text{max}}Dmax 是在等位基因频率给定的情况下 D 值可能达到的最大值。

1
2

如果 D>0D > 0D>0：Dmax=min⁡(pA⋅(1−pB),(1−pA)⋅pB)D_{\text{max}} = \min(p_A \cdot (1 - p_B), (1 - p_A) \cdot p_B)Dmax=min(pA⋅(1−pB),(1−pA)⋅pB)
如果 D<0D < 0D<0：Dmax=min⁡(pA⋅pB,(1−pA)⋅(1−pB))D_{\text{max}} = \min(p_A \cdot p_B, (1 - p_A) \cdot (1 - p_B))Dmax=min(pA⋅pB,(1−pA)⋅(1−pB))

3. 计算 r² 值

r² 表示两个位点之间的 LD 强度，计算公式为：

1	r2=D2pA⋅(1−pA)⋅pB⋅(1−pB)r^2 = \frac{D^2}{p_A \cdot (1 - p_A) \cdot p_B \cdot (1 - p_B)}r2=pA⋅(1−pA)⋅pB⋅(1−pB)D2

安装

可以参考华大github的安装方式：BGI-shenzhen/PopLDdecay: PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format(VCF) files (github.com)

这里演示conda安装

conda install hcc::poplddecay 
#值得注意的是如果主环境zlib不匹配，环境solving 会fail
#这个时候新建一个环境就可以
#conda create -n poplddecay  python=3.9
#在这个环境内安装就可以了

安装好测试

参数说明

1. 输入和输出文件
-InVCF <str>: 输入的 SNP VCF 文件（压缩格式，通常为 .vcf.gz）。
-InGenotype <str>: 输入的 SNP 基因型文件（另一个输入格式选项，通常为纯文本格式）。
-OutStat <str>: 输出的统计结果文件，包含 LD 值（例如，r^2 或 D'）。

2.可选参数
-SubPop <str>: 子群体样本文件列表（默认是所有样本）。这个文件应列出感兴趣的样本名称，一行一个样本。
-MaxDist <int>: 两个 SNP 之间的最大距离（以千碱基为单位），超过这个距离的 SNP 将不计算 LD 值。默认值为 300 kb。
-MAF <float>: 最小次常见等位基因频率（MAF）过滤阈值。默认值为 0.005，即次常见等位基因频率必须至少为 0.5% 才会被包含在分析中。
-Het <float>: 杂合子等位基因比例的最大阈值。默认值为 0.88，即如果杂合子比例超过 88%，该 SNP 将被过滤掉。
-Miss <float>: 缺失等位基因比例的最大阈值。默认值为 0.25，即如果缺失比例超过 25%，该 SNP 将被过滤掉。
-EHH <str>: 如果需要运行 EHH（Extended Haplotype Homozygosity）区域衰减分析，设置起始位点（默认值为 `NA`，即不运行 EHH 分析）。
-OutFilterSNP: 输出最终用于计算的 SNP 列表。
-OutType <int>: 输出类型。可以选择以下值：
 - 1: 仅输出 r^2值（默认）。
 - 2: 输出 r^2 和 D' 值。
 - 3: 输出成对 LD 结果。
 - 更多输出类型请参阅帮助文档，参数范围为 1-8。

- -help: 显示更多帮助信息和版本号（当前版本为 hewm2008 v3.42）。
1
2
3
4
5
6
7
8
9
10
11
12
13

   

### 计算举例

```sh
# 直接处理 GATK VCF 文件
./bin/PopLDdecay -InVCF SNP.vcf.gz -OutStat LDdecay
# 对于 PLINK 格式的 [.ped .map] 文件，先将 PLINK 文件转换为 PopLDdecay 所需的基因型格式再用 PopLDdecay 进行 LD 衰减计算：
perl bin/mis/plink2genotype.pl -inPED in.ped -inMAP in.map -outGenotype out.genotype
./bin/PopLDdecay -InGenotype out.genotype -OutStat LDdecay
# 计算 VCF 文件中子群体的 LD 衰减将子群体样本名称放入 GroupA_sample.list 文件中，然后运行：
./bin/PopLDdecay -InVCF -OutStat -SubPop GroupA_sample.list

可视化举例

# 单个群体
perl bin/Plot_OnePop.pl -inFile LDdecay.stat.gz -output Fig
# 单个群体，多个染色体
perl bin/Plot_OnePop.pl -inList Chr.ResultPath.List -output Fig
# 多个群体
perl bin/Plot_MutiPop.pl -inList Pop.ResultPath.list -output Fig

部分参考文档

分群体分析

这个是官方说明文档告诉我们的群体list文件包含这个群体样本ID，一个ID一行。

在绘图时multi群体list文件每行有两列，第一列是LD计算结果文件绝对路径，第二列是分组类名。目测是tab分隔的。

这个时候就有人疑惑了，conda安装的plot_pop脚本去哪里了，坑人的家伙出坏教程，还是用官方github编译教程靠谱，plot_pop脚本就在src目录下，这个conda安装教程莫非要多此一举绕去github找源码？

conda和docker给我们这些没有root权限无法编译的带来的便利是巨大的。不要质疑conda，plot_pop就在软件安装的路径上级目录。

1	which PopLDdecay

分染色体分析

先介绍单一群体，官方举例的情况是vcf文件已经按照chr拆分，可以使用基本命令分析。推测不支持拆分染色体功能，每条染色体单独分析效果好一些，我推荐用vcftools拆分染色体。

#计算
PopLDdecay -InVCF chr1.vcf.gz -OutStat chr1.stat.gz
for i in {1..23}; do (PopLDdecay -InVCF chr$i.vcf.gz -OutStat chr$i.stat.gz)& done
#绘图
perl  Plot*pl -inList chr.list -output out
#list文件是分析的结果文件一个一行。

举个例子

多群体按照多染色体计算

好戏登场了，这个方案我觉得很不错，跟上面的类似先按照单染色体计算，不过要注意加入群体分组list，不要用单一群体vcf再指定不是这个群体内的ID哦。我觉得你也可以按照这两个类先用vcftools拆分两遍在计算，但是比较消耗存储资源。

在可视化时是先把每个群体的单独绘制，再使用multiplot整合到一起

#计算
PopLDdecay -InVCF chr1.vcf.gz -OutStat chr1.stat.gz -SubPop GroupA_sample.list
PopLDdecay -InVCF chr1.vcf.gz -OutStat chr1.stat.gz -SubPop GroupB_sample.list
#可视化
perl Plot_OnePop.pl -inList  GroupA.chr.list -output ga.chr.cat  #chr.list上面已经演示
perl Plot_OnePop.pl -inList  GroupB.chr.list -output gb.chr.cat  #chr.list上面已经演示
perl Plot_MultiPop.pl -inList  multi.list -output ab.chr.out  #multi.list 按照分群体计算那一步的办法自己改改

最后给出我的方案（未拆分chr），不清楚可以咨询

1 2	bash ld.sh --vcf input.vcf.gz --pop_dir populations_dir --out_dir output_dir # 选择性修改 --maxdist 300 --maf 0.005 --het 0.88 --miss 0.25

ld.sh代码如下：

#!/bin/bash
#create by lx
function usage() {
    echo "Usage: bash $0 --vcf <vcf_file> --pop_dir <pop_dir> --out_dir <output_dir> [--maxdist <int>] [--maf <float>] [--het <float>] [--miss <float>]"
    echo "  --vcf        <str>   Input VCF file"
    echo "  --pop_dir    <str>   Directory containing population sample files"
    echo "  --out_dir    <str>   Output directory"
    echo "Optional parameters:"
    echo "  --maxdist    <int>   Max distance (kb) between two SNPs [default: 300]"
    echo "  --maf        <float> Min minor allele frequency filter [default: 0.005]"
    echo "  --het        <float> Max ratio of het allele filter [default: 0.88]"
    echo "  --miss       <float> Max ratio of miss allele filter [default: 0.25]"
    exit 1
}

maxdist=300
maf=0.005
het=0.88
miss=0.25
vcf_file=""
pop_dir=""
output_dir=""

while [[ $# -gt 0 ]]; do
  case "$1" in
    --vcf)
      vcf_file="$2"
      shift 2
      ;;
    --pop_dir)
      pop_dir="$2"
      shift 2
      ;;
    --out_dir)
      output_dir="$2"
      shift 2
      ;;
    --maxdist)
      maxdist="$2"
      shift 2
      ;;
    --maf)
      maf="$2"
      shift 2
      ;;
    --het)
      het="$2"
      shift 2
      ;;
    --miss)
      miss="$2"
      shift 2
      ;;
    *)
      echo "Unknown option: $1"
      usage
      ;;
  esac
done

if [[ -z "$vcf_file" || -z "$pop_dir" || -z "$output_dir" ]]; then
    echo "Error: Missing required parameters."
    usage
fi

mkdir -p "$output_dir"
total_files=$(ls -1 "$pop_dir"/*.txt | wc -l)
current_file=0

for pop_file in "$pop_dir"/*.txt; do
    pop_name=$(basename "$pop_file" .txt)
    output_stat="$output_dir/$pop_name.stat"

    (
        PopLDdecay -InVCF "$vcf_file" -OutStat "$output_stat" -SubPop "$pop_file" -MaxDist "$maxdist" -MAF "$maf" -Het "$het" -Miss "$miss"
        current_file=$((current_file + 1))
        echo "Completed $current_file of $total_files: $pop_name"
    ) &
done

wait
echo "程序执行完毕，输出文件保存在 $output_dir"

plot不再展示。