使用AdmixTools(admixr)做D-statistic

2021-04-17 本文已影响0人 DumplingLucky

ADMIXTOOLS是David Reich实验室开发的软件包，在群体遗传学中，用于计算和测试群体混合统计，验证杂交事件。R包admixr是实现ADMIXTOOLS的一种直接方法，它本质上是ADMIXTOOLS软件的接口，要在ADMIXTOOLS安装好的环境下运行。
admixr分析案例：

1. 安装ADMIXTOOLS

wget https://github.com/DReichLab/AdmixTools/archive/refs/heads/master.zip
unzip master.zip
cd AdmixTools-master/src/
make clobber
make all
make install
#可执行程序位于../bin路径下

2. 配置R环境

#为了确保R可以找到ADMIXTOOLS程序包，在R环境文件(.Renviron)中添加ADMIXTOOLS程序包的路径
PATH=<path to the location of ADMIXTOOLS programs>
#在R环境中
install.packages("admixr")
install.packages("tidyverse")
library(admixr)
library(tidyverse)

3. 产生输入文件

ADMIXTOOLS需要特征(eigenstrat)文件，要将vcf文件转换为eigenstrat。

（1）eigenstrat格式包含三种文件：

ind文件 - 指定每个样本的唯一名称，性别(可选-可以用“U”表示“undefined”)和样本标签；
snp文件 - 指定SNP，REF/ALT等位基因的位置；

1: SNP string ID
2: chromosome
3: genetic distance
4: position along a chromosome
5: reference allele
6: alternative allele
geno文件 - SNP数据(一行一个位点，每个样本一个字符):
0：个体是纯合子ALT
1：个体是杂合子
2：个体是纯合子REF
9：数据丢失
上述三种文件信息均包含在VCF文件中。

ind.file

snp.file

geno.file

（2）转换方法

#vcf格式转换成eigenstrat格式，Joana脚本
sh convertVCFtoEigenstrat.sh test.vcf.gz
#这个脚本链接我还没有找到......随后更新

4. 实战分析

（1）输入文件准备

#使用admixr的download_data()方法下载测试数据
(prefix <- download_data(dirname = tempdir()))
[1] "/var/folders/t7/9gjtb6m92flbnp930618vt3r0000gn/T//RtmpveArl0/snps/snps"
list.files(path = dirname(prefix), pattern = basename(prefix), full.names = TRUE)
[1] "/var/folders/t7/9gjtb6m92flbnp930618vt3r0000gn/T//RtmpveArl0/snps/snps.geno"
[2] "/var/folders/t7/9gjtb6m92flbnp930618vt3r0000gn/T//RtmpveArl0/snps/snps.ind" 
[3] "/var/folders/t7/9gjtb6m92flbnp930618vt3r0000gn/T//RtmpveArl0/snps/snps.snp"

（2）D统计

D统计（或称为ABBA-BABA检验）是检测杂交事件的简单方法。比如W，X，Y，Z具有 ((W，X), Y), Z) 的系统发育关系，其中Z为外群，检测Y是否有向X的基因渗入。 ADMIXTOOLS中实现的用于计算D统计量的公式基于等位基因频率。

#确定W群体
pops <- c("French", "Sardinian", "Han", "Papuan", "Khomani_San", "Mbuti", "Dinka")
#D检验
result <- d(W = pops, X = "Yoruba", Y = "Vindija", Z = "Chimp", data = snps)
head(result)
W   X   Y   Z   D   stderr  Zscore  BABA    ABBA    nsnps
French  Yoruba  Vindija Chimp   0.0313  0.006933    4.510   15802   14844   487753
Sardinian   Yoruba  Vindija Chimp   0.0287  0.006792    4.222   15729   14852   487646
Han Yoruba  Vindija Chimp   0.0278  0.006609    4.199   15780   14928   487925
Papuan  Yoruba  Vindija Chimp   0.0457  0.006571    6.953   16131   14721   487694

D - 𝐷统计值
stderr - 𝐷统计量的标准误差
Zscore - 𝐷从0开始的标准误数，即我们否定无杂交的零假设的强烈程度
BABA/ABBA - 观察到的位点模式的计数
nsnps - 用于给定计算的SNP数

（3）ggplot2绘制结果

ggplot(result, aes(fct_reorder(W, D), D, color = abs(Zscore) > 2)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_errorbar(aes(ymin = D - 2 * stderr, ymax = D + 2 * stderr))

结果解读：非洲人的𝐷值与0并无显著差异，这意味着该数据中非洲人没有尼安德特人的基因渗入，与零假设相符。另一方面，所有非洲人以外的人群拒绝零假设，这表明尼安德特人与当今非洲人以外人的祖先存在基因交流。

（4）其他检验方法

admixr也可以做f4，f4-ratio，f3检验等。

#f4
result <- f4(W = pops, X = "Yoruba", Y = "Vindija", Z = "Chimp", data = snps)
#f4-ratio
result <- f4ratio(X = pops, A = "Altai", B = "Vindija", C = "Yoruba", O = "Chimp", data = snps)
ggplot(result, aes(fct_reorder(X, alpha), alpha, color = abs(Zscore) > 2)) +
  geom_point() +
  geom_errorbar(aes(ymin = alpha - 2 * stderr, ymax = alpha + 2 * stderr)) +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(y = "Neandertal ancestry proportion", x = "present-day individual")

结果解读：当今的非非洲人占尼安德特人血统的2-3％，巴布亚新几内亚人的尼安德特人血统比例更高-超过4％。这与早期的研究相一致，早期的研究表明在当今的巴布亚人的祖先中还有其他的古混合事件。

#f3
pops <- c("French", "Sardinian", "Han", "Papuan", "Mbuti", "Dinka", "Yoruba")
result <- f3(A = pops, B = pops, C = "Khomani_San", data = snps)
#根据f3值对人口标签进行排序
ordered <- filter(result, A == "Mbuti", B != "Mbuti") %>% arrange(f3) %>% .[["B"]] %>% c("Mbuti")
#成对f3热图绘制
result %>% 
filter(A != B) %>%
  mutate(A = factor(A, levels = ordered),
         B = factor(B, levels = ordered)) %>%
  ggplot(aes(A, B)) + geom_tile(aes(fill = f3))

当我们基于成对的𝑓3统计值对热图标签进行排序时，展示出了已知的人口分化顺序（即Sar先分开，然后是Mbuti等）。
参考：
https://speciationgenomics.github.io/ADMIXTOOLS_admixr/
https://bodkan.net/admixr/articles/tutorial.html