5. GWAS：群体结构

#5. GWAS：群体结构| 来源: 网络整理| 查看: 265

群体结构是指材料的亚群分化情况，会导致标记间的非连锁关联，进而导致关联分析结果出现假阳性。

地理隔离、人工选择、移民和遗传漂变等都可能导致群体分化。

是指遗传变异在物种或群体中的一种非随机分布；

将各材料归到每个亚群，计算每个材料基因组变异源于第K个亚群的可能性，用Q值表示，Q值越大，表明改材料来自这个亚群的可能性越大，一般可以用来推断祖先群，个体血缘组成，还有杂交事件;

常用软件：Admixture、Structure、Frappe等。

随着技术的发展，Structure速度较慢，无法满足大量分子标记计算的需求，因此，admixture逐渐成为群体结构分析的主流软件。本文将介绍如何通过admixture进行群体结构计算。

1.下载及安装 1.1 下载地址

http://dalexander.github.io/admixture/index.html

1.2 安装 $ tar xvf admixture_linux-1.3.0.tar.gz $ cd your/path/admixture_linux-1.3.0 # 调用：./admixture # 查看帮助：./admixture --help 2. 群体结构计算 2.1 整理成admixture所需的.ped(12recode)格式

在plink中将vcf文件转换成admixture所需的.ped或.bed格式:

$ cd your/path/plink1.9 $ ./plink --vcf genotype.vcf --allow-extra-chr --recode12 --out genotype12 --autosome-num 27

--vcf 输入文件名 --allow-extra-chr 允许其他格式染色体，如scaffold --recode12 二进制编码 --out 输出文件名 --autosome-num 设置染色体数目，默认人类染色体数

2.2 Admixture $ cd your/path/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 xxx RUN normal login 4*compute11 sh adm.sh Aug 24 01:14

每个K值都会生成两个文件，.P和.Q P：储存推断的祖先种群的等位基因频率 Q：每个样本中各个祖先种群所占的百分比。

3. 最佳分群数确定及可视化 3.1 确定最佳分群数

查看cv值，cv error最小的K值为最佳分群数。

$ grep -h CV log*.out CV error (K=10): 0.65873 CV error (K=2): 0.71095 CV error (K=3): 0.63424 CV error (K=4): 0.68598 CV error (K=5): 0.67584 CV error (K=6): 0.66818 CV error (K=7): 0.66301 CV error (K=8): 0.66083 CV error (K=9): 0.65919 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") > install.packages(c("ggplot2","gridExtra","label.switching","tidyr","remotes"),repos="https://cloud.r-project.org") > remotes::install_github('royfrancis/pophelper') > library("pophelper") > tbl=read.table("genotype.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文件，文件共包含三列，每行为一个样本，三列中哪一个数值最大，则这个样本属于哪一个亚群。

引用转载请注明出处，如有错误敬请指出。

【本文地址】

公司简介

联系我们