本系列文章采用的數據集與代碼來自https://github.com/MareesAT/GWA_tutorial。
該教程獲得了許多人的推薦，是一份很詳細的step-by-step guide。
本文將介紹該教程中的QC部分（1_QC_GWAS.zip），后續或將繼續添加有關QC的其他細節。
0. 目錄
1. 准備
首先，使用下述命令即可將該Github項目下載到本地：
git clone https://github.com/MareesAT/GWA_tutorial.git
下載后，將文件1_QC_GWAS.zip解壓縮即可得到該部分的教程文件與數據。
其中，教程文件為1_Main_script_QC_GWAS.txt。
由於本教程還需使用plink軟件，plink 1.9版本的下載頁面見https://www.cog-genomics.org/plink2/。只需在頁面選擇對應系統版本的二進制軟件壓縮包下載，解壓后即可直接使用。
此外，環境中還應裝有R語言。
2. 數據簡介
該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
3. 分析步驟
第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”。
作者提供的R腳本（gender_check.R）可以繪制所有樣本F值的分布情況，並分別給出了數據集中標記為男性、女性的個體的F值分布。
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-xd6SOjgc-1596263025956)(https://upload-images.jianshu.io/upload_images/15241933-651738542a5ce8fd.png?imageMogr2/auto-orient/strip|imageView2/2/w/1240)]
接下來，作者給出了兩種刪除異常個體的方法。此處介紹第一種。
#從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<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
檢查結果中，每一列分別代表：
- 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
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-EFz5VZ4L-1596263025964)(https://upload-images.jianshu.io/upload_images/15241933-fb39879b597bb5da.png?imageMogr2/auto-orient/strip|imageView2/2/w/1240)]
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
通常，應使用基於特定家庭的方法來分析基於家庭的數據。在本教程中，出於說明目的，將無關樣本中的相關性視為隱含的相關性。
在本步驟中，我們旨在從數據集中刪除所有“關聯性”。為了證明大多數相關性是親子關系所致，我們將僅在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
4. 結語
以上就是該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