admixture 群体结构分析
阅读原文时间:2023年07月10日阅读:9

tructure是与PCA、进化树相似的方法,就是利用分子标记的基因型信息对一组样本进行分类,分子标记可以是SNP、indel、SSR。相比于PCA,进化树,群体结构分析可明确各个群之间是否存在交流及交流程度

1 软件安装

conda install -c bioconda admixture

admixture
****                   ADMIXTURE Version 1.3.0                  ****
****                    Copyright 2008-2015                     ****
****           David Alexander, Suyash Shringarpure,            ****
****                John  Novembre, Ken Lange                   ****
****                                                            ****
****                 Please cite our paper!                     ****
****   Information at www.genetics.ucla.edu/software/admixture  ****

Usage: admixture <input file> <K>
See --help or manual for more advanced usage.

2 简单使用

第一步:将VCF变为plink格式

## vcftools
vcftools --gzvcf SNP.vcf.gz --plink --out test
## 没有压缩,则为--vcf

## 或者plink 直接转
plink --vcf SNP.vcf.gz--recode --out test--const-fid --allow-extra-chr

 # --vcf vcf 或者vcf.gz
 # --recode 输出格式ped(默认bed)
 # --out 输入前缀
 # --const-fid  添加群体信息familyID sampleID)(将family设置为0);
 # --allow-extra-chr 允许非标准染色体编号

$\color{red}{**}$plink转化,map文件,SNP那列为点,vcftools 则不是,但是ref为0

# 如果用vcftools转换, 重新添加染色体
paste <(cut -f2 Test.map |awk -F ":" '{print $1}') <(cut -f2-4 Apple.map)  >map
## 如果用plink 转换, 重新添加SNP编号
awk '{x+=1}{print $1"\tSNP"x"\t"$3"\t"$4}' Test.map >map
#重命名
mv map Test.map

因为plink本身是针对人类进行开发的,所以在运行时,加上--allow-extr-chr。此外对于vcf中的sampleID(familyID_sampleID), plink 默认为下划线分隔(也可以通过参数--id-delim进行修改),分别作为family ID和sampleID。但是一般我们的样本并不是那样命名的,所以可以添加--double-id参数,将familyID和sampleID命名为相同,或者--const-fid,将familyID命名为0,表明-9

第二步 plink进行过滤,得到bed文件

plink --allow-extra-chr --noweb -file test --geno 0.05 --maf 0.05 --hwe 0.0001 --make-bed --out test1
# --noweb 不显示网页

第三步:寻找合适的K值

for i in {1..7};do admixture --cv test1.bed $i |tee log${i}.out;done
## 根据CV error 确定K值
grep -h 'CV'  log*.out
CV error (K=1): 0.22317
CV error (K=2): 0.15018
CV error (K=3): 0.12804
CV error (K=4): 0.12109
CV error (K=5): 0.12656

当k=4时,cv error 最小,选择4

此外,

如果SNP数据集非常大,则可以随机选择SNP进行K值选择分析,比如随机选取20000个SNP进行分析,每个K值跑20次,确定最终的k值,然后分析

当利用plink转的格式中,在运行上述命令是出现以下报错

Invalid chromosome code!  Use integers

将*.bim中的第一列改为数值就可以了

第四步:多线程

admixture test1.bed 4 -j 5
## j, 线程数

第五步:作图

ta1 = read.table("test1.4.Q")
head(ta1)
barplot(t(as.matrix(ta1)),col = rainbow(3),
        xlab = "Individual",
        ylab = "Ancestry",
        border = NA)

根据LD进行筛选SNP

首先 筛选合格的LD位点

## 对vcf进行
plink --vcf SNP.vcf.gz --indep-pairwise 100 50 0.2 -out Test-id-maf0.05-LD --allow-extra-chr

## 或者对bed进行操作也可以
plink --bfile  test1 --indep-pairwise 100 50 0.2

## --indep-pairwise 100 50 0.2; 100Kb窗口,50步长,0.2LD强度

然后进行提取

plink  --bfile test1 --extract plink.prune.in --make-bed --out test2

然后在进行和上述一样的分析即可

欢迎扫码交流

参考

--Admixture使用说明文档cookbook

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章