接上文GWAs——全基因组关联分析（质控1），此数据集模拟的是祖先来自欧洲西北部的犹他州居民，所以需要将没有欧洲背景的个体从数据集中剔除，即控制群体结构（Population Stratification，群体分层）。

一、控制群体结构

1、创建工作目录

不同于教程的目录结构，为了和文章整体结构对应，我将控制群体结构的工作目录放在了质控目录下（1_QC_GWAS）。由于此步骤产生文件数较多，将在子目录（pop_str）下再创建几个目录。

#创建工作目录
cd /{your directory}/GWAs/1_QC_GWAS/
#创建控制群体结构的主目录
mkdir pop_str

#将所需文件置入工作目录
cp Relatedness/HapMap_3_r3_12.* pop_str/
cp Relatedness/indepSNP.prune.in pop_str/
mv MDS_merged.R /pop_str/

我们将使用1000 Genomes Project（1KG）的数据检查人口分层，此数据包含有629个来自不同地域、民族的个体，数据集压缩包（ALL.2of4intersection.20100804.genotypes.vcf.gz）大小约61G。但是使用wget的方法下载速度过慢。
压缩包需要数据转换，将其下载到1KG_format目录下。

#创建存放1KG数据及进行格式转换的目录
cd pop_str
mkdir 1KG_format

方法1：使用wget下载

#使用wget下载
cd 1KG_format
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz

方法2：使用ascp加速下载

#使用加速下载器下载
cd 1KG_format
ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz ./

另外，也可以下载到本地，再挂载到服务器。

2、数据格式化

下载的1KG数据是**.vcf**格式的压缩文件，使用PLINK 1.9可以将数据转换为PLINK可处理的二进制文件（.bed，.fam和.bim），作者文章中说配套代码可在PLINK 1.7运行，但我在尝试时并没有成功。

#格式转换,这个过程会比较慢
/{your directory}/GWAs/plink-20220402/plink --vcf ALL.2of4intersection.20100804.genotypes.vcf.gz --make-bed --out ALL.2of4intersection.20100804.genotypes

在完成二进制转化的文件中，部分SNPs并没有以r开头的编号（.bim文件），取而代之的是“.”，可以使用以下代码查看，注意不要随便使用类似vim的工具打开数据文件，文件过大会导致故障发生。

#快速查看
zmore ALL.2of4intersection.20100804.genotypes.vcf.gz

因此使用以下命令为每个SNPs添加独立编号，参数格式为“@:#[b37]$1,$2”，其中@#不可替换，分别表示染色体编号（@）和位点（#），此例子中使用的编号格式为“染色体位置：位点[b37]allele1，allele2”。

#添加编号，注意，command行中的#不是忽略，#和后面的内容也是command的一部分
/nfs_fs/nfs1/sunwu/GWAs/plink-20220402/plink --bfile ALL.2of4intersection.20100804.genotypes --set-missing-var-ids @:#[b37]\$1,\$2 --make-bed --out ALL.2of4intersection.20100804.genotypes_no_missing_IDs

3、对1KG数据进行质控

作为background数据集的1KG数据同样需要质控，包括控制检出率和控制MAF两步，大致步骤同前。

（1）控制检出率

搭建工作环境。

#返回pop_str目录
cd /{your directory}/GWAs/1_QC_GWAS/pop_str/
#创建质控目录1KG_QC
mkdir 1KG_QC
cd 1KG_QC
#控制检出率目录
mkdir 1KG_missing
#置入数据集
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/1KG_format/ALL.2of4intersection.20100804.genotypes_no_missing_IDs.* 1KG_missing/

开始质控。切记，与之前的质控相同，必须先过滤缺失个体信息SNPs，再过滤缺失SNPs的个体。

cd 1KG_missing
#初检，检出missingness > 0.2
/{your directory}/GWAs/plink-20220402/plink --bfile ALL.2of4intersection.20100804.genotypes_no_missing_IDs --geno 0.2 --allow-no-sex --make-bed --out 1kG_MDS
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS --mind 0.2 --allow-no-sex --make-bed --out 1kG_MDS2

#二次检出，检出missingness > 0.02
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS2 --geno 0.02 --allow-no-sex --make-bed --out 1kG_MDS3
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS3 --mind 0.02 --allow-no-sex --make-bed --out 1kG_MDS4

最终数据集是以1kG_MDS4为名的三个二进制文件，下文使用1kG_MDS4.*表示。

（2）控制MAF

搭建工作环境。

cd /{your directory}/GWAs/1_QC_GWAS/pop_str/1KG_QC/
#创建控制MAF的工作目录
mkdir 1KG_MAF
#置入文件
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/1KG_QC/1KG_missing/1kG_MDS4.* 1KG_MAF/

开始质控。

cd 1KG_MAF
#移除MAF < 0.05的个体
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS4 --maf 0.05 --allow-no-sex --make-bed --out 1kG_MDS5

最终数据集1kG_MDS5.*

4、数据整合

进行群体结构控制需要将两个数据（1KG和HapMap）集进行整合，但是二者的数据内容和格式还存在很大的区别。为此，需要将二者的格式和内容做调整。

（1）提取相同的SNPs位点并更新位点信息

搭建工作目录。

cd /{your directory}/GWAs/1_QC_GWAS/pop_str/
#数据整合的主目录
mkdir merge
#提取SNP位点的工作目录
cd merge
mkdir extract

#置入数据集
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/1KG_QC/1KG_MAF/1kG_MDS5.* extract/
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/HapMap_3_r3_12.* extract/

从1KG数据集（1kG_MDS5）中提取在HapMap数据集（HapMap_3_r3_12.bim）中出现过的SNPs位点。

cd extract/
#提取HapMap_3_r3_12出现过的位点编号
awk '{print$2}' HapMap_3_r3_12.bim > HapMap_SNPs.txt
#从1KG数据集中提取位点
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

从HapMap数据集（HapMap_3_r3_12）中提取在1KG数据集（1kG_MDS6.bim）出现过的位点。注意由于后面要根据HapMap数据集统一每个SNPs的位点信息，此过程需要文本数据集中的file.map文件，在命令中添加“–recode”参数可以将数据集写入文本数据集文件。

#提取SNPs位点编号
awk '{print$2}' 1kG_MDS6.bim > 1kG_MDS6_SNPs.txt
#提取位点
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_3_r3_12 --extract 1kG_MDS6_SNPs.txt --recode --make-bed --out HapMap_MDS

通过位点信息的提取筛选，现在数据集HapMap_MDS和数据集1kG_MDS6包含有相同的SNPs位点集。但是还需要根据HapMap数据集的map文件（HapMap_MDS.map）更新1KG数据集中的SNPs位点信息。

#将HapMap数据集中map内记录的SNPs编号和位点信息写入txt文件
awk '{print$2,$4}' HapMap_MDS.map > buildhapmap.txt

#更新1KG数据集的位点信息
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS6 --update-map buildhapmap.txt --make-bed --out 1kG_MDS7

最终数据集：

（2）整合前检查

在数据整合之前需要确定两数据集是可以整合的。为此需要再一轮的检查。
搭建工作环境。

#创建工作目录
cd /{your directory}/GWAs/1_QC_GWAS/pop_str/merge/
mkdir prep_merge

#置入文件
mv extract/HapMap_MDS.* prep_merge/
mv extract/1kG_MDS7.* prep_merge/

数据集整合前，对数据的检查主要包括三个步骤：

第一步：Reference Genome

确定1KG数据集和HapMap数据集具有相似的Reference Genome。

cd prep_merge

#提取1KG里面的编号和对应的allele1的genotype
awk '{print$2,$5}' 1kG_MDS7.bim > 1kg_ref-list.txt
#强制将allele1的genotype设为参考Reference Genome
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_MDS --reference-allele 1kg_ref-list.txt --make-bed --out HapMap-adj
#此步骤会产生大量的warning信息，但无伤大雅

第二步：调整链方向

链问题（strand issue，序列倒置），即两个数据集之间的链方向不一致，需要翻转。

#分别提取量数据集的SNPs信息（编号，allele1，allele2）
awk '{print$2,$5,$6}' HapMap-adj.bim > HapMap-adj_tmp
awk '{print$2,$5,$6}' 1kG_MDS7.bim > 1kGMDS7_tmp

#排列两文件内的SNPs，并检出排列后只出现一次的SNPs
sort 1kGMDS7_tmp HapMap-adj_tmp |uniq -u > all_differences.txt

代码解析：sort将两文件内容全部合并排列；“|”是管道符，将排列后的结果传递给uniq -u将排列结果中只出现过一次的结果检出，将检出结果写入txt文件。
all_differences.txt中记录，两个数据集共有1624个SNPs信息不一致，即每个数据集有812个SNPs记录和另一个文件不一致。

#提取两个数集SNPs信息不一致的SNPs编号
awk '{print$1}' all_differences.txt | sort -u > flip_list.txt

导致这812个SNPs数据不一致的原因是其存在正-负链不一致的问题，为此，我们可以根据Reference Genome（1kg_ref-list.txt）的信息，把两个数据集的链统一，或全为正链，或全为负链。

#调整正负链
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap-adj --flip flip_list.txt --reference-allele 1kg_ref-list.txt --make-bed --out corrected_hapmap

再次检查：经过调整翻转后的两个数据集仍会存在一些SNPs信息不一致，这些SNPs极有可能是存在问题，需要检出。

#提取HapMap数据集调整后的SNPs信息
awk '{print$2,$5,$6}' corrected_hapmap.bim > corrected_hapmap_tmp

#检出不一致项（1KG和HapMap）
sort 1kGMDS7_tmp corrected_hapmap_tmp |uniq -u  > uncorresponding_SNPs.txt
awk '{print$1}' uncorresponding_SNPs.txt | sort -u > SNPs_for_exlusion.txt

从文件uncorresponding_SNPs.txt中可以发现，两个数据集共有84个SNPs信息仍然不一致，也就是42个相同编号但在两数据集信息不同的SNPs位点（SNPs_for_exlusion.txt）。在后续的步骤中我们需要将他们移除（根据SNPs_for_exlusion.txt文件）。

第三步：移除问题项

分别从1KG和HapMap数据集中移除存在问题的SNPs。

#从HapMap的数据集中移除
/{your directory}/GWAs/plink-20220402/plink --bfile corrected_hapmap --exclude SNPs_for_exlusion.txt --make-bed --out HapMap_MDS2

#从1KG数据集中移除
/{your directory}/GWAs/plink-20220402/plink --bfile 1kG_MDS7 --exclude SNPs_for_exlusion.txt --make-bed --out 1kG_MDS8

（3）数据整合

现在我们得到了整合前的最终数据集：

开始整合两个数据集。

#搭建工作目录
cd /{your directory}/GWAs/1_QC_GWAS/pop_str/merge/
mkdir mer
#置入数据集
mv /prep_merge/HapMap_MDS2.* mer
mv /prep_merge/1kG_MDS8.* mer

#整合数据集
cd mer
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_MDS2 --bmerge 1kG_MDS8.bed 1kG_MDS8.bim 1kG_MDS8.fam --allow-no-sex --make-bed --out MDS_merge2

最终整合数据集：MDS_merge2.*

5、多维尺度降维分析

MDS分析（Multi-dimensional Scaling，多维尺度降维），简单理解就是将研究的所有个体置入一个多维度的空间中，并依托数学和计算机手段，将这个多维尺度的空间压缩降维为一个二维可视化的面上，且同时能最大程度上反映各数据点之间的聚类关系（内在关系），从而达到聚类的目的。
本次示例中，对数据集涵盖的每个个体进行降维聚类，从而直观的反映群体结构。

（1）数据准备

搭建工作环境：

#创建工作目录
cd /{your directory}/GWSs/1_QC_GWAS/pop_str/
mkdir MDS

#置入数据集和文件
mv indepSNP.prune.in /MDS
mv MDS_merged.R /MDS
mv merge/mer/MDS_merge2.* /MDS
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/merge/prep_merge/HapMap_MDS.fam /MDS
mv /{your directory}/GWAs/1_QC_GWAS/pop_str/merge/extract/HapMap_3_r3_12.* /MDS

整合的数据集MDS_merge2.*里面还存在非独立的SNPs位点（主要来自1KG数据集），需要根据indepSNP.prune.in（独立或低关联SNPs目录）计算IBD和IBS相关数值（后续降维聚类需要）。

#计算IBD和IBS相关数值
/{your directory}/GWAs/plink-20220402/plink --bfile MDS_merge2 --extract indepSNP.prune.in --genome --out MDS_merge2

（2）进行降维聚类

现在使用整合数据集（MDS_merge2）和IBS和IBD相关数据（MDS_merge2.genome）进行降维聚类（10维，由–mds-plot参数设置）。

#开始聚类，10维
/{your directory}/GWAs/plink-20220402/plink --bfile MDS_merge2 --read-genome MDS_merge2.genome --cluster --mds-plot 10 --out MDS_merge2

此步骤共产出3大类文件，其中MDS_merge2.nosex记录了性别不明的个体ID，聚类结果记录在MDS_merge2.cluster1，MDS_merge2.cluster2，MDS_merge2.cluster3中，MDS结果记录在MDS_merge2.mds

（3）MDS可视化

为了获取1KG数据集中每个个体的背景信息，还需要下载20100804.ALL.panel，此文件主要由三到四列组成，第一列为个体ID，第二列为背景信息，第三、四列为测序的平台。下面是代码。

#下载1KG个体背景信息
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/20100804.ALL.panel

这些不同的个体来自于四个大的超群体：AFR、AMR、ASN和EUR。现在将属于这些超群体的群体内的个体归入到超群体下，得到1KG的racefile。

#提取个体ID和所属群体
awk '{print$1,$1,$2}' 20100804.ALL.panel > race_1kG.txt
#将群体内的个体划分到四个超群体下
sed 's/JPT/ASN/g' race_1kG.txt>race_1kG2.txt
sed 's/ASW/AFR/g' race_1kG2.txt>race_1kG3.txt
sed 's/CEU/EUR/g' race_1kG3.txt>race_1kG4.txt
sed 's/CHB/ASN/g' race_1kG4.txt>race_1kG5.txt
sed 's/CHD/ASN/g' race_1kG5.txt>race_1kG6.txt
sed 's/YRI/AFR/g' race_1kG6.txt>race_1kG7.txt
sed 's/LWK/AFR/g' race_1kG7.txt>race_1kG8.txt
sed 's/TSI/EUR/g' race_1kG8.txt>race_1kG9.txt
sed 's/MXL/AMR/g' race_1kG9.txt>race_1kG10.txt
sed 's/GBR/EUR/g' race_1kG10.txt>race_1kG11.txt
sed 's/FIN/EUR/g' race_1kG11.txt>race_1kG12.txt
sed 's/CHS/ASN/g' race_1kG12.txt>race_1kG13.txt
sed 's/PUR/AMR/g' race_1kG13.txt>race_1kG14.txt

提取我们自己数据集的个体ID、谱系信息，并在第三列加入OWN标签（对应1KG数据集race_1kG14.txt中的第三列），得到我们自己的racefile。

#创建我们自己的racefile
awk '{print$1,$2,"OWN"}' HapMap_MDS.fam>racefile_own.txt

最后将1KG和我们自己的racefile串联起来，并添加列名。

#串联两个racefile
cat race_1kG14.txt racefile_own.txt | sed -e '1i\FID IID race' > racefile.txt

代码解析：cat将两个racefile陆续弹出，sed接受弹出结果，-e参数指定其按自定义脚本’1i\FID IID race’ 运行，其中1i表示在第一行前插入，反斜线后内容为插入内容。

使用R脚本进行绘图。

#绘图
Rscript MDS_merged.R 

R脚本解析：MDS_merged.R

#读取MDS结果文件
data<- read.table(file="MDS_merge2.mds",header=TRUE)
#读取合并racefile文件
race<- read.table(file="racefile.txt",header=TRUE)
#数据合并
#将两个list以data的关键词IID进行merge
#并使用race的FID填充merge后的FID
datafile<- merge(data,race,by=c("IID","FID"))
#查看前6行
head(datafile)
#PDF绘图
pdf("MDS.pdf",width=7,height=7)
#以前两个主成分C1和C2为x和y坐标绘图
#不同超群体以不同颜色表示
for (i in 1:nrow(datafile))
{
if (datafile[i,14]=="EUR") {plot(datafile[i,4],datafile[i,5],type="p",xlim=c(-0.1,0.2),ylim=c(-0.15,0.1),xlab="MDS Component 1",ylab="MDS Component 2",pch=1,cex=0.5,col="green")}
par(new=T)
if (datafile[i,14]=="ASN") {plot(datafile[i,4],datafile[i,5],type="p",xlim=c(-0.1,0.2),ylim=c(-0.15,0.1),xlab="MDS Component 1",ylab="MDS Component 2",pch=1,cex=0.5,col="red")}
par(new=T)
if (datafile[i,14]=="AMR") {plot(datafile[i,4],datafile[i,5],type="p",xlim=c(-0.1,0.2),ylim=c(-0.15,0.1),xlab="MDS Component 1",ylab="MDS Component 2",pch=1,cex=0.5,col=470)}
par(new=T)
if (datafile[i,14]=="AFR") {plot(datafile[i,4],datafile[i,5],type="p",xlim=c(-0.1,0.2),ylim=c(-0.15,0.1),xlab="MDS Component 1",ylab="MDS Component 2",pch=1,cex=0.5,col="blue")}
par(new=T)
if (datafile[i,14]=="OWN") {plot(datafile[i,4],datafile[i,5],type="p",xlim=c(-0.1,0.2),ylim=c(-0.15,0.1),xlab="MDS Component 1",ylab="MDS Component 2",pch=3,cex=0.7,col="black")}
par(new=T)
}

#添加直线，lty指定线条类型
#h参数指定y坐标的垂直点
#v参数指定x坐标的垂直点
abline(v=-0.035,lty=3)
abline(h=0.035,lty=3)
#添加标签
legend("topright", pch=c(1,1,1,1,3),c("EUR","ASN","AMR","AFR","OWN"),col=c("green","red",470,"blue","black"),bty="o",cex=1)

可视化结果（Fig. 1）中显示，我们自己的选择实验个体与EUR（欧洲）个体聚到了一类，而我们的个体是祖先来自欧洲西北部的犹他州居民，这符合我们的实验预期，无需移除个体，但是为了展示移除离群个体的流程，下文会展示移除离群值的方法。

Fig. 1：MDS结果可视化

（4）移除离群值

剔除离群值是基于前文的MDS结果和可视化进行的，根据可视化结果，依据前两个成分（C1、C2），选择合适的阈值过滤掉“过分”离群的个体。

#将非离群值写入EUR_MDS_merge2
#筛选条件为C1 < -0.04且C2 > 0.03
awk '{ if ($4 <-0.04 && $5 >0.03) print $1,$2 }' MDS_merge2.mds > EUR_MDS_merge2

移除EUR_MDS_merge2中的符合实验预期的个体。（本步骤并没有移除任何个体，因为没有离群值，但是在实际研究中，如果存在阈值外的个体，那么他们会被EUR_MDS_merge2内的list排除在外）。

#提取合格个体
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_3_r3_12 --keep EUR_MDS_merge2 --make-bed --out HapMap_3_r3_13

至此，我们已经获得最终数据集：
HapMap_3_r3_13.bed
HapMap_3_r3_13.bim
HapMap_3_r3_13.fam

二、提取协变量

仅对经过所有质控步骤的数据进行MDS，生成的10个成分轴（C1至C10）就可以作为协变量。
搭建工作环境，注意：此目录存放的文件将是可进行关联分析的最终文件。

#创建工作目录
cd /{your directory}/GWAs/1_QC_GWAS/pop_str/
mkdir covariates

#置入数据集
mv MDS/HapMap_3_r3_13.* covariates/
#置入独立位点名录
mv MDS/indepSNP.prune.in covariates/

开始MDS，生成协变量。

#计算IBS、IBD值
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_3_r3_13 --extract indepSNP.prune.in --genome --out HapMap_3_r3_13
#进行MDS
/{your directory}/GWAs/plink-20220402/plink --bfile HapMap_3_r3_13 --read-genome HapMap_3_r3_13.genome --cluster --mds-plot 10 --out HapMap_3_r3_13_mds

但是此时的协变量还不能被直接读取利用，需要对其格式进行修改。

#提取协变量
awk '{print$1, $2, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13}' HapMap_3_r3_13_mds.mds > covar_mds.txt

至此，我们已经完成所有质控步骤，获得所有进行GWAs的文件和数据集：
HapMap_3_r3_13.bed
HapMap_3_r3_13.bim
HapMap_3_r3_13.fam
covar_mds.txt

Ending！