tructure是与PCA、进化树相似的方法,就是利用分子标记的基因型信息对一组样本进行分类,分子标记可以是SNP、indel、SSR。相比于PCA,进化树,群体结构分析可明确各个群之间是否存在交流及交流程度
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.
第一步:将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位点
## 对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
然后在进行和上述一样的分析即可
手机扫一扫
移动阅读更方便
你可能感兴趣的文章