群体遗传分析方法：LD,FST,eQTL

Posted 2023-05-02

参考技术A

LD(连锁不平衡):计算使用plink,
FST(遗传分化指数):计算使用vcftools,可视化分为箱线图和散点图，单组比较使用在染色体上的散点图，多组比较使用箱线图。
FST的原理，计算方法，可视化的方法 https://www.jianshu.com/p/bb0beec0ed63
haploPS和XP-EHH
平均测序深度：
等位基因频率：
vcftools使用说明：1. https://www.dazhuanlan.com/2020/03/04/5e5f850672ddf/

Fst值的取值范围是[0,1]，最大值为1表明两个群体完全分化，最小值为0表明群体间无分化。
实际使用FST<0--0.05,表示群体分化很小；0.05--0.15，中等程度的分化；0.15--0.25，较大的分化；0.25以上，分化很大。
其实Fst分析就是看两个群体之间分化程度的一种方法，Fst值越大（越接近1）表明两个群体间分化程度越高，亲缘关系越远；Fst值越小（越接近0）表明群体间分化程度越低，亲缘关系越近。
分为按照 SNP单点计算 和 滑动窗口模式计算 ，一般使用的是滑动窗模式
####### SNP单点计算

参数说明：
--vcf 输入vcf文件
--weir-fst-pop 输入群体的群体ID名，该文件必须是txt格式，每个ID占一行。
--fst-window-size 500000 --fst-window-step 50000 ，这里窗口设置为500kb，步长设置为50kb。根据情况调整窗口大小。
计算haploPS、XP-EHH、 Fst，正向选择分析方法寻找性状相关的位点（拿几个群体的重测序数据，轻松发核心中文或者<3的SCIE）

π，核苷酸多样性，越大说明核苷酸多样性越高，越低说明两个座位DNA序列差异越小。
vcftools

计算和可视化， 参考

找到重测序的数据，基因分型，找到单倍型，call SNP,过滤，注释snp,可视化，分别计算haploPS,XP-EHH,Fst.求交集，公共基因进行注释，富集分析。

https://www.jianshu.com/p/6e6d54d7483e

https://www.genedenovo.com/news/364.html

LDSC的安装使用

vcftools --vcf F2_3.vcf --plink --out F2_3 #把vcf格式转为plink格式

群体遗传分析—LD连锁不平衡

参考技术A 当位于某一座位的特定等位基因与另一座位的某一等位基因同时出现的概率大于群体中因随机分布的两个等位基因同时出现的概率时，就称这两个座位处于连锁不平衡状态（linkage disequilibrium）。

D 是 LD（连锁不平衡） 的基本单位，度量观察到的单倍型频率与平衡状态下期望频率的偏差。虽然D能够很好的表达LD的基本含义，但是由于其严格依赖于等位基因频率（allele frequency），故不适合应用于表述实际的LD强度，尤其是进行不同研究的LD值的相互比较。几个常用于度量LD的符号中，最重要的是D'和r2，两者都是基于D，各有各的特点及用途。

LD计算方法如下：
1、设有两个位点（A、B），等位基因分别是A、a、B、b，在群体中对应频率f(A)、f(a)、f(B)和f(b)

2、两个位点共有四种单倍型AB、Ab、aB、ab，对应频率f(AB)、f(Ab)、f(aB)和f(ab)

3、计算：Dab=f(AB)-f(A)*f(B)

当Dab=0时，处于连锁平衡状态；

当Dab≠0时，处于连锁不平衡状态。

LD度量：

当Dab>0，|D'|=(Dab)2/min(f(AB), f(ab))；

当Dab<0，|D'|=(Dab)2/min(f(Ab), f(aB))；

r2=(Dab)2/(f(A) f(a) f(B)*f(b))；

D'=0, r2=0时处于完全连锁平衡状态；

D'=1,r2=1时处于完全连锁不平衡状态；

从0-1度量越高，LD越高，如果两个位点连锁，连锁程度也越强。

r2和D'反映了LD的不同方面。r2包括了重组和突变，而D'只包括重组史。D'能更准确地估测重组差异，但样本较小时，低频率等位基因组合可能无法观测到，导致LD强度被高估，所以D'不适合小样本群体研究；

LD衰减作图中通常采用r2来表示群体的LD水平；Haplotype Block中通常采用D'来定义Block；迁移、突变、选择、有限的群体大小以及其他引起等位基因频率改变的因素，这些都会引起LD的改变。

plink2 ：

( https://www.cog-genomics.org/plink2 )

haploview ：

( https://www.broadinstitute.org/haploview/haploview )

plink计算R2值的命令行(基于vcf)：

--vcf 指定输入的文件为vcf格式，如果是bed格式文件，使用--bfile接文件前缀，如果数据是ped 、map格式，使用 --map接.map文件，--ped接.ped文件

--allow-no-sex 表示允许没有性别信息

--maf 指定maf阈值

--geno 指定缺失率阈值，与我们的完整度意思相反，0.2对应的完整度为0.8

--r2表示计算r2值

--ld-window 表示计算LD的区间，表示距离小于这个值的标记对都要进行LD的计算

--ld-window-r2

这个参数只能和--r2参数搭配使用，默认值为0.2, 对输出结果进行过滤，只输出R2大于该参数值的LD分析结果。

haploview 计算R2值的命令行：

java -jar Haploview.jar -nogui -memory 10240 -info test.hapmap.info -pedfile test.hapmap.ped -out test -maxdistance 500 -minGeno 0.5 -minMAF 0.05 -missingCutoff 0.5 -hwcutoff 0 -dprime

LD的衰减指位点间由连锁不平衡到连锁平衡的演变过程；LD衰减的速度在不同物种间或同物种的不同亚群间，往往差异非常大。所以，通常会使用1个标准——“LD衰减距离”来描述LD衰减速度的快慢。

LD衰减距离通常指的是：当平均LD系数r2 衰减到一定大小的时候，对应的物理距离。“一定大小”是这个定义的关键点，但没有特别统一的标准，在不同文章中标准不同。常见的标准包括：

a）LD系数降低到最大值的一半；

b）LD系数降低到0.5以下；

c）LD系数降低到0.1以下；

d）LD系数降低到基线水平（注意，不同物种的基线值是不同的）。

所以，下次你在文章中看到“LDdecay distance is XXkb”的时候，不要忘了看看文章使用的标准是什么。

Nature Biotechnology 30, 105–111 (2012) doi:10.1038/nbt.2050

值的获取：成对计算指定距离范围内的所有SNP的r2 值，按区间取平均

1、判断GWAS所需标记量，决定GWAS的检测效力以及精度；

GWAS标记量 = 基因组大小/LD衰减距离

2、辅助分析进化与选择

在同一个连锁群上，LD衰减的慢说明该群体受到选择。一般来说，野生群体比驯化改良群体LD衰减快，异花授粉植物比自花授粉植物LD衰减快。比如玉米：地方品种1kb，自交系2kb，商用自交系100kb。

单体型块，即连锁不平衡区域，是指同一条染色体上处于连锁不平衡状态的一段连续的区域。单体型块分析可以用于筛选tag SNP、确定候选基因的范围等。

如果GWAS检测到显著关联的区间，可以通过进一步绘制局部的LD单体型块图，来进一步判断显著相关的SNP和目标基因间是否存在强LD关系。