GWAs——全基因组关联分析二（质控2）

接上文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