群体结构分析

目的：对群体结构和亲缘关系进行评估以确定使用的统计模型和获得相应的矩阵

评估内容（遗传上差异过大应剔除，相似性高的保留其一）

Structure是与PCA、进化树相似的方法，就是利用分子标记的基因型信息对一组样本进行分类，分子标记可以是SNP、indel、SSR… …当然，对于重测序应用的最多的还是SNP。进化树和PCA本质上都是计算样本序列间的差异程度，然后利用两两差异度聚类（进化树）或降维（PCA）来实现对样本的分类。structure本质上使用了与PCA、进化树完全不同的思路。

一、群体结构：构建系统发育树（必备）

同一物种内序列差异不大构建NJ树(mega)，序列差异较大，不同种构建ML树(RAxML),贝叶斯树(ExaBayes)

常用的编辑和显示树图的软件有TreeView、FigTree、MEGA、ITOL（http://itol.embl.de/）、R包（ggtree、APE）等。

进化树美化：https://itol.embl.de/

二、主成分分析

参考：「群体遗传学实战」第二课: 画出和文章几乎一样的PCA图

主成分分析(PCA)是一种线性降维方法，能从纷繁复杂的数据中抽离出关键因素，用来区分不同的样本。这里介绍以下三款：

Plink:https://www.cog-genomics.org/plink/

GCTA:https://cnsgenomics.com/software/gcta/

EIGENSOFT:https://github.com/DReichLab/EIG

无论使用哪款软件，始终都要记得它们最初为人类基因组设计，因此一定要关注和染色体有关的参数。例如Plink一定要加上 --allow-extra-chr允许非标准的染色体编号。

1、plink进行PCA分析

#将vcf文件转换为plink的ped和map格式 

plink --vcf watermelon_414acc_SNP2.vcf.gz --recode --out watermelon_414acc --const-fid --allow-extra-chr

#将ped和map转换为bed、bim、fam格式

plink --allow-extra-chr --file watermelon_414acc --noweb --make-bed --out watermelon_414acc  

进行PCA计算

plink --allow-extra-chr --threads 20 -bfile watermelon_414acc --pca 20 --out watermelon_414acc

这一步会得到两个文件，一个是以.eigenval结尾的文件，记录特征值，用来计算每个PC所占的比重。另一个是.eigenvec结尾的文件，记录特征向量，用于坐标轴。随后用ggplot2绘制pca图。

用plink的EIGENSOFT可以在分析时自动剔除离群值、做LD拟合等等。经过初步的分析，确定要剔除的离群值后，用plink的--keep和--remove参数剔除个体，值分析特定样本。它们要求输入文件为两列，一列是样本所在的群体编号，一列是样本编号。以剔除CN样本为例，它的样本编号为WM439，所在群体是0。代码如下：

echo '0\tWM439' > remove.txt

plink --remove remove.txt --allow-extra-chr -bfile watermelon_414acc --pca 20 --out watermelon_414acc_no_cn 

2、GCTA进行PCA分析

要求染色体的编号一定得是数字，不然在读取bim结尾的文件时，一定会报illegal chr number错误。

三、Admixture绘制群体结构图

参考：GWAS 学习之admixture、、Admixture：一款快速分析群体遗传结构的软件、、【一起学生信】群体结构图形—structure堆叠图

群体遗传学中测出很多个个体，得到了最终的SNP vcf文件，需要将其分成群体，看哪几个物种聚在一起。一般使用的软件是STRUCTURE，但是STREUTURE运行速度极慢，后面frappe软件提升了速度，但是也不是很快；admixture凭借其运算速度，成为了主流的分析软件。

admixture的下载：http://software.genetics.ucla.edu/admixture/download.html   不需要安装，解压出来即可使用

admixture输入文件：经plink处理后的.bed文件

1、创建PLINK 的-bed 文件

plink --noweb --file bw.chip --geno 0.05 --maf 0.05 --mind 0.05 --make-bed --out bw.chip  #这里包含了数据质控过滤

生成bed、bim、fam三种文件

2、用Admixture软件构建群体遗传结构和群体世系信息

K是样本所包含的亚群或者祖先数，若不知理想的K值，可以设定K=1,2,3,4,5，用admixture进行计算，命令如下：

for K in 1 2 3 4 5 6 7; do admixture --cv hapmap3.bed $K | tee log${K}.out; done

3、提取出CV值，确定最佳分群数

grep -h CV log*.out

提取CV值后，可以得到上一步得到的不同K值的错误率数据，一般认为CV error最小值为最佳K值。

4、使用R画图

tbl=read.table("hapmap3.3.Q")

pdf("/USER/XXX/Project/xj/reseq/result/11.admixture/Q7.pdf")

barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA) 

dev.off()

本例绘制的为K=3时的图，结果如下（供参考）

1种颜色代表该个体的血统是纯正的；多种颜色表明该个体很可能是从多个祖先亚群杂交而来