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 允许非标准染色体编号
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
然后在进行和上述一样的分析即可
网友评论