5. GWAS：群体结构——Admixture
-
群体结构是指材料的亚群分化情况，会导致标记间的非连锁关联，进而导致关联分析结果出现假阳性。
-
地理隔离、人工选择、移民和遗传漂变等都可能导致群体分化。
-
是指遗传变异在物种或群体中的一种非随机分布；
-
将各材料归到每个亚群，计算每个材料基因组变异源于第K个亚群的可能性，用Q值表示，Q值越大，表明改材料来自这个亚群的可能性越大，一般可以用来推断祖先群，个体血缘组成，还有杂交事件;
-
常用软件：Admixture、Structure、Frappe等。
随着技术的发展，Structure速度较慢，无法满足大量分子标记计算的需求，因此，admixture逐渐成为群体结构分析的主流软件。本文将介绍如何通过admixture进行群体结构计算。
1.下载及安装
1.1 下载地址
http://dalexander.github.io/admixture/index.html
1.2 安装
cd /gss1/home/sunwei/software
tar xvf admixture_linux-1
cd /gss1/home/sunwei/software/admixture_linux-1.3.0
# 调用：./admixture
# 查看帮助：./admixture --help
2. 群体结构计算
2.1 整理成admixture所需的.ped(12recode)格式
在plink中将vcf文件转换成admixture所需的.ped或.bed格式，这里以.ped为例
cd /gss1/home/fzhang/sunwei/plink1.9
./plink --vcf root.id.vcf --allow-extra-chr --recode12 --out root12 --autosome-num 27
--vcf 输入文件名
--allow-extra-chr 允许其他格式染色体，如scaffold
--recode12 二进制编码
--out 输出文件名
--autosome-num 设置染色体数目，默认人类染色体数
2.2 Admixture
cd /gss1/home/fzhang/sunwei/admixture_linux-1.3.0
# 创建任务文件
vim adm.sh
# vim 文件名
# i 输入 左下角出现insert，可以输入
for K in 2 3 4 5 6 7 8 9 10; do ./admixture --cv root12.ped $K | tee log${K}.out; done
# ESC键 insert消失
# :wq
# 提交任务
bsub -n 4 -o log sh adm.sh
#查看任务
$ bjobs
JOBID USER STAT QUEUE FROM_HOST EXEC_HOST JOB_NAME SUBMIT_TIME
913421 fzhang RUN normal login 4*compute11 *f.gz;done Aug 24 01:14
每个K值都会生成两个文件，.P和.Q
P：储存推断的祖先种群的等位基因频率
Q：每个样本中各个祖先种群所占的百分比。
3. 最佳分群数确定及可视化
3.1 确定最佳分群数
查看cv值，K=3时，cv error最小，最佳分群为3。
$ grep -h CV root_log*.out
CV error (K=10): 0.4305
CV error (K=2): 0.40957
CV error (K=3): 0.39887
CV error (K=4): 0.39909
CV error (K=5): 0.40694
CV error (K=6): 0.40875
CV error (K=7): 0.40982
CV error (K=8): 0.41662
CV error (K=9): 0.41964
3.2 群体结构可视化
将CV结果复制粘贴至Excel中，绘制折线图。图中可看出最佳分群数为K=3。
在R中绘制群体结构图
提供几个我喜欢的配色：
K=3 "#FF4500","#9ACD32","#6495ED"
K=4 "#336666","darkred","steelblue","#CC9933"
K=5 "#FF4500","#5F7A61","#6495ED","#986D8E","#F6D167"
将K=3时的.Q文件拷贝至Windows中
> setwd("D:/研究生/数据/GWAS/群体结构")
> library("ggplot2")
> library("pophelper")
> tbl=read.table("root.3.Q")
> pdf("admixture.pdf",width = 9,height = 3)
> colorpal =c("#FF4500","#9ACD32","#6495ED")
> cols=rep(colorpal,700)
> barplot(t(as.matrix(tbl)), col=cols, xlab="", ylab="Ancestry",border = NA)
> dev.off()
3.3 确定样本属于哪个亚群
当确定最佳分群数是3时，打开K=3时的.Q文件，文件共包含三列，每行为一个样本，三列中哪一个数值最大，则这个样本属于哪一个亚群。
引用转载请注明出处，如有错误敬请指出。
