GWAS学习笔记（一）：质量控制（QC）

本系列文章采用的数据集与代码来自https://github.com/MareesAT/GWA_tutorial。
该教程获得了许多人的推荐，是一份很详细的step-by-step guide。

本文将介绍该教程中的QC部分（1_QC_GWAS.zip），后续或将继续添加有关QC的其他细节。

准备

首先，使用下述命令即可将该Github项目下载到本地：

git clone https://github.com/MareesAT/GWA_tutorial.git

下载后，将文件1_QC_GWAS.zip解压缩即可得到该部分的教程文件与数据。

1_QC_GWAS.zip解压缩后得到的文件

其中，教程文件为1_Main_script_QC_GWAS.txt。

由于本教程还需使用plink软件，plink 1.9版本的下载页面见https://www.cog-genomics.org/plink2/。只需在页面选择对应系统版本的二进制软件压缩包下载，解压后即可直接使用。

此外，环境中还应装有R语言。

数据简介

该Github教程使用了可免费获取的HapMap数据：hapmap3_r3_b36_fwd.consensus.qc。编写者模拟了一组二进制表型特征，并将其添加到该数据集中，并命名为HapMap_3_r3_1。未经添加的原始HapMap数据可见http://hapmap.ncbi.nlm.nih.gov/downloads/genotypes/2010-05_phaseIII/plink_format/

典型的plink数据集包括三个文件：.bed文件、.bim文件和.fam文件

bed文件：二进制文件，主要是存储等位基因信息。它开头前三个字节永远是0x6c, 0x1b, 和0x01，接下来就是V组N/4个字节的序列，这里V是指遗传变异的个数，N是指样本数，假如N无法被4整除，那么将N/4的结果取整后加1作为各组的字节数，编码信息如下：

00：基因型是bim文件中allele 1的纯合子

01：基因型缺失

10：基因型是杂合子

11：基因型是bim文件中allele 2的纯合子

bim文件：文本文件。包含染色体编号（默认可用1-22、X、Y。可采用扩展编号）、SNP编号、位点的摩尔距离（可用0代表不知道）、物理位置、allele 1（常为次等位基因）、allele 2（常为主等位基因）。其中allele用0代表缺失。

fam文件：文本文件。包含Family ID、Individual ID、Paternal ID（父本ID）、Maternal ID（母本ID）、Sex（雄性为1，雌性2，未知为0）、Phenotype。这里的Phenotype取值可为1（对照组）、2（实验组/病例）、-9/0（表示实验组/对照组表型缺失）。如果出现了 {-9, 0, 1, 2}之外的值，则表型会被读取为数量性状。

有公众号文章（小麦穗粒数转录组分析(四)----使用plink进行关联分析）指出，对于小麦，可将染色体改名为数字（如1A->1，1B->2类推）来适应plink的染色体编号规则。

plink 1.9支持将染色体用字符表示。在由VCF格式向plink转换时，使用参数“--allow-extra-chr”即可正常转换形如“chr1A”的染色体名称。

该GitHub项目中提供的脚本可在简单修改后适用于其他数据集的研究。但由于脚本是针对二进制表型（binary outcome measure）开发的，因此并不适用于数量性状的研究（需要进一步修改脚本）。

Note, most GWAS studies are performed on an ethnic homogenous population, in which population outliers are removed. The HapMap data, used for this tutorial, contains multiple distinct ethnic groups, which makes it problematic for analysis.
Therefore, we have selected only the EUR individuals of the complete HapMap sample for the tutorials 1-3. This selection is already performed in the HapMap_3_r3_1 file from our GitHub page.

有关质量控制的进一步细节，可参见项目作者发表的文章A Tutorial on Conducting Genome-Wide Association Studies: Quality Control and Statistical Analysis

第1步：SNP缺失率过滤

计算数据集缺失率：

plink --bfile HapMap_3_r3_1 --missing

该步会输出文件plink.imiss和plink.lmiss，分别包含每个样本个体的SNP缺失率，和对每个SNP而言存在缺失的样本数量的比例。

随后，可用项目提供的R脚本绘制直方图

Rscript --no-save hist_miss.R

随后，可以按照阈值删除缺失率过高的SNP位点和样本个体：

# 删除缺失率>0.2的SNP
plink --bfile HapMap_3_r3_1 --geno 0.2 --make-bed --out HapMap_3_r3_2

# 删除缺失率>0.2的个体
plink --bfile HapMap_3_r3_2 --mind 0.2 --make-bed --out HapMap_3_r3_3

# 删除缺失率>0.02的SNP
plink --bfile HapMap_3_r3_3 --geno 0.02 --make-bed --out HapMap_3_r3_4

# 删除缺失率>0.02的个体
plink --bfile HapMap_3_r3_4 --mind 0.02 --make-bed --out HapMap_3_r3_5

在过滤环节，项目作者指出，先从不严格的阈值开始过滤是个好习惯。

第2步：检测、过滤性别差异

该步利用X染色体的近交（纯合度）估计，并依据先验确定的女性的受试者的F值必须<0.2、男性的受试者的F值必须> 0.8进行检测。

plink --bfile HapMap_3_r3_5 --check-sex 

Rscript --no-save gender_check.R

该步骤会生成文件plink.sexcheck、plink.log和plink.hh。在plink.sexcheck中记录了每个样本个体的状态（STATUS）和F值等。其中，不满足前述要求个体的状态会被标记为“PROBLEM”。

plink.sexcheck

作者提供的R脚本（gender_check.R）可以绘制所有样本F值的分布情况，并分别给出了数据集中标记为男性、女性的个体的F值分布。

数据集中全部样本的F值分布

接下来，作者给出了两种删除异常个体的方法。此处介绍第一种。

#从plink.sexcheck中筛选带有字符“PROBLEM”的行->筛选结果中提取第1、2列写入sex_discrepancy.txt文件
grep "PROBLEM" plink.sexcheck| awk '{print$1,$2}'> sex_discrepancy.txt

#删除样本
#plink使用Family ID和Individual ID定位到特定个体
plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6

第3步：次等位基因频率（MAF）过滤

该步骤将删除MAF过低的SNP，即删除变异过于罕见的位点。项目作者首先筛选出了位于22对常染色体上的SNP。

awk '{ if ($1 >= 1 && $1 <= 22) print $2 }' HapMap_3_r3_6.bim > snp_1_22.txt
plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7

随后用plink生成MAF分布数据，并绘制成直方图。

plink --bfile HapMap_3_r3_7 --freq --out MAF_check
Rscript --no-save MAF_check.R

MAF分布

最后，删除MAF<0.05的位点。通常，GWAS项目采用的MAF阈值在0.01-0.05之间，取决于样本大小。

plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8

第4步：删除不符合Hardy-Weinberg平衡的SNP

首先，检查所有SNP的HWE p值分布

plink --bfile HapMap_3_r3_8 --hardy

HWE检查结果

检查结果中，每一列分别代表：

snp 所在染色体
snp 名称
test的名称
Minor allele code
Major allele code
具体数据也就是 AA Aa aa 的个数
观察到的2pq 的值
期望的2pq的值
对这个数据进行卡方检验，看显不显著

选择HWE p-value低于0.00001的SNP绘图，以放大严重偏离的SNP

awk '{ if ($9 <0.00001) print $0 }' plink.hwe>plinkzoomhwe.hwe
Rscript --no-save hwe.R

HWE分布与放大之后的分布

HWE的过滤部分，作者分了两步。首先是默认参数的HWE过滤，plink仅会过滤对照组。添加参数“--hwe-all”后，plink将过滤全部样本。

#默认参数，仅过滤对照组
plink --bfile HapMap_3_r3_8 --hwe 1e-6 --make-bed --out HapMap_hwe_filter_step1

#过滤全部样本
plink --bfile HapMap_hwe_filter_step1 --hwe 1e-10 --hwe-all --make-bed --out HapMap_3_r3_9

更多理论背景可见项目作者的文章A Tutorial on Conducting Genome-Wide Association Studies: Quality Control and Statistical Analysis

第5步：控制杂合率偏差

本步骤将生成受试者杂合率分布图，并删除杂合率偏离均值超过3倍标准差的个体。

首先需要对不相关的SNP进行杂合率检查。项目作者提供了文件inversion.txt，其中记录了在高LD区域中的SNP编号。检查时，排除了这些SNP，并使用参数--indep-pairwise来精简SNP。参数参数“ 50 5 0.2”分别代表：窗口大小，每步移动窗口经过的SNP数量（步长），以及一个SNP在所有其他SNP上同时回归的多重相关系数（multiple correlation coefficient for a SNP being regressed on all other SNPs simultaneously）(即r^2 threshold)。

plink --bfile HapMap_3_r3_9 --exclude inversion.txt --range --indep-pairwise 50 5 0.2 --out indepSNP

该命令生成两个文件indepSNP.prune.in和indepSNP.prune.out。prune.in文件中包含的就是通过筛选条件的、独立的SNP位点，而prune.out则记录具有LD的SNP位点。

plink --bfile HapMap_3_r3_9 --extract indepSNP.prune.in --het --out R_check

该命令将生成文件R_check.het和R_check.log。其中，“.het”文件格式见下表。

列名	含义
FID	Family ID
IID	Within-family ID
O(HOM)	Observed number of homozygotes，实际观测到的纯合子数量
E(HOM)	Expected number of homozygotes，预测的纯合子数量
N(NM)	Number of non-missing autosomal genotypes，常染色体上非缺失的基因型数量
F	Method-of-moments F coefficient estimate，矩量法估计的近交系数F

下述脚本读取了R_check.het文件，并绘制了杂合率分布图

Rscript --no-save check_heterozygosity_rate.R

杂合率分布直方图

项目作者提供了脚本heterozygosity_outliers_list.R以生成杂合率与平均值相差超过三倍标准差的样本。结果保存为fail-het-qc.txt，并对该文件进行处理，处理结果交由plink进行删除。

#提取杂合率与平均值相差超过三倍标准差的样本所在行
Rscript --no-save heterozygosity_outliers_list.R

#去掉文件中的引号并提取前两行，以匹配plink的格式
sed 's/"// g' fail-het-qc.txt | awk '{print$1, $2}'> het_fail_ind.txt

#plink根据提取结果删除对应样本
plink --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10

第6步：清除隐含相关的个体

该步骤进行样本间的亲缘关系估计，并将清除未标记明显亲子关系的个体亲缘关系过近的个体。

下述命令将进行样本间亲缘关系的估计，生成pihat_min0.2.genome和pihat_min0.2.log。

plink --bfile HapMap_3_r3_10 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2

“.genome”格式文件介绍如下表

列名	含义
FID1	First sample's family ID
IID1	First sample's within-family ID
FID2	Second sample's family ID
IID2	Second sample's within-family ID
RT	Relationship type inferred from .fam/.ped file，从fam/ped文件推断的样本关系
EZ	IBD sharing expected value, based on just .fam/.ped relationship
Z0	P(IBD=0)
Z1	P(IBD=1)
Z2	P(IBD=2)
PI_HAT	Proportion IBD, i.e. P(IBD=2) + 0.5*P(IBD=1)，IBD比例
PHE	Pairwise phenotypic code (1, 0, -1 = case-case, case-ctrl, and ctrl-ctrl pairs, respectively)
DST	IBS distance, i.e. (IBS2 + 0.5*IBS1) / (IBS0 + IBS1 + IBS2)
PPC	IBS binomial test
RATIO	HETHET : IBS0 SNP ratio (expected value 2)

IBD全称Identity By Descent, 又叫做血缘同源，指的是两个个体中共有的等位基因来源于共同祖先
IBS全称Identity By State, 又叫做状态同源，指的是两个个体中共有的等位基因序列相同。
对于某个等位基因，IBS state只要求allel的个数相同即可，而IBD state则进一步要求相同的allel来自于共同祖先。
此文章介绍了IBS与IBD的更多细节：IBD血缘同源简介

其中，列RT的取值可为：

列名	含义
FS	full siblings，全同胞，指有着相同父母的个体之间的关系
HS	half siblings，半同胞，指同父异母或同母异父的个体之间的关系
PO	parent-offspring，亲子代关系
OT	other，其他
UN	不相关的个体

已知HapMap包含亲子关系，下述代码将使用z值可视化这些亲子关系。

#如果第8列即Z1>0.9，输出该行
awk '{ if ($8 >0.9) print $0 }' pihat_min0.2.genome > zoom_pihat.genome

Rscript --no-save Relatedness.R

Pihat分布与样本亲子关系在z0和z1上的可视化

通常，应使用基于特定家庭的方法来分析基于家庭的数据。在本教程中，出于说明目的，将无关样本中的相关性视为隐含的相关性。

在本步骤中，我们旨在从数据集中删除所有“关联性”。为了证明大多数相关性是亲子关系所致，我们将仅在founders（数据集中没有父母的个人，家系创建人）中进行pihat筛选。

#筛选出仅含founders的数据集
plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11

#再次进行亲缘关系估计
plink --bfile HapMap_3_r3_11 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2_in_founders

文件pihat_min0.2_in_founders.genome显示，在排除所有非founders之后，HapMap数据中仅剩下1对pihat大于0.2的个体对。根据Z值，这可能是完整的同胞或DZ双胞胎对。值得注意的是，他们在HapMap数据中没有获得相同的家庭身份（FID）。

对于每对pihat> 0.2的“相关”个体，项目作者建议删除call rate最低的个体。

#对于筛选出的仅含founders的数据集，再次计算缺失率
plink --bfile HapMap_3_r3_11 --missing

在缺失率计算结果中，手动找到“相关”个体对中缺失率更高的个体（在本例中，是13291 NA07045），并构建文件0.2_low_call_rate_pihat.txt，在文件中输入以下内容：

13291  NA07045

保存文件，并用如下命令删除对应样本。（如有多个“相关”对，则寻找每对中缺失率更高的个体记录在本文件中，每个个体占一行）

plink --bfile HapMap_3_r3_11 --remove 0.2_low_call_rate_pihat.txt --make-bed --out HapMap_3_r3_12

结语

以上就是该GWAS教程的质控部分了。本部分生成的文件中，下述文件将在下一部分中继续使用：

The bfile HapMap_3_r3_12 (i.e., HapMap_3_r3_12.fam,HapMap_3_r3_12.bed, and HapMap_3_r3_12.bim）
indepSNP.prune.in