1. 程式人生 > 其它 >【翻譯】全基因組關聯分析教程:質量控制和統計分析【第二部分:軟體介紹&質量控制】

【翻譯】全基因組關聯分析教程:質量控制和統計分析【第二部分:軟體介紹&質量控制】

原文標題:A tutorial on conducting genome-wide association studies: Quality control and statistical analysis
原文連結:https://pubmed.ncbi.nlm.nih.gov/29484742/

2 | 軟體

質量控制流程和統計分析將使用免費的、開源的全基因組管理分析工具集PLINK 1.07版(Purcell等人,2007)進行說明,該工具集可以從 http://zzz.bwh.harward.edu/plink 下載。PLINK 1.9測試版包含相同的選項,同時執行的速度要快得多 https://www.cog-genomics.org/plink/1.9/

。由於PLINK 1.9目前是測試版本,我們在本教程中使用了官方版本的PLINK。但是,也可以使用PLINK 1.9完成教程的所有內容。儘管本文中討論的某些步驟可以在R等傳統統計軟體包中執行,但專門用於分析遺傳資料的軟體包使用起來要方便得多。除了PLINK,還有許多其他不錯的選擇可用於分析SNP資料,例如Genebel(Aulchenko等人,2007)和 SNPTEST(Marchini等人,2007)。 此外,允許在基於家庭的GWAS中測試關聯的方法的軟體也被開發出來了(Chen和Yang,2010,Ott等)。我們建議使用基於GNU/Linux的計算機環境,儘管許多選項也可以通過Windows版本的PLINK獲得。可以在
http://www.ee.surrey.ac/Teaching/Unix
上找到shell和命令列的基本介紹。Github上的示例指令碼生成的所有圖形都將使用免費的開源程式語言R( https://www.r-project.org/ )獲得。

2.1 | 資料格式

PLINK可以讀取文字格式的檔案,也可以讀取二進位制格式的檔案。由於讀取大型的文字檔案可能耗費大量的時間,於是我們推薦使用二進位制格式的檔案。文字形式的PLINK資料由兩個檔案組成:一個包含了個體及其基因型的資訊(*.ped),另一個包含了遺傳標記資訊(*.map;詳情見圖1)。相應的,二進位制形式的PLINK資料由三個檔案組成,其中一個二進位制檔案包含了個體的識別符號(IDs)和基因型(*.bed

),另兩個文字檔案各自包含了個體資訊(*.fam)和遺傳標記資訊(*.bim,詳情見圖1)。例如,在雙向情感障礙的研究中,*.bed檔案包含了所有患者和健康對照的基因分型結果;*.fam檔案包含了與受試者相關的資料(與研究中其他參與者的家庭關係、性別和臨床診斷);而*.bim檔案包含了有關SNP物理位置的資訊。使用協變數(covariates)的分析通常需要第四個檔案,其中包含了每個人的相關協變數的值(參見圖1)。

圖1 各種常用的PLINK檔案。SNP即單核苷酸多型性

2.2 | PLINK基礎命令

PLINK是一款命令列程式,因此,使用PLINK需要一個活動的shell來等待命令。這可以通過游標前的提示($或者>)來識別。通常情況下,當前目錄的路徑會顯示在提示符之前,如圖2所示。當前目錄是PLINK使用中的一個核心概念,因為在預設情況下,PLINK將從當前目錄載入資料檔案,並將執行產生的結果檔案儲存在當前目錄中。當前目錄可以通過傳統的Unix命令(通常是cd)來改變為任何目錄。在提示之後,通過輸入plink關鍵字表示我們需要使用PLINK。如果PLINK沒有安裝在標準目錄中,則必須在命令前面加上PLINK的安裝目錄路徑,例如,/usr/local/bin/plink

在輸入的plink關鍵字之後,控制PLINK的工作流程的其他選項將緊隨其後,並以空格分隔。這些選項全都以兩個橫線(--)開頭。需要提供的第一個選項之一是資料檔案的格式和名稱:對於文字檔案,使用--file {your_file},二進位制檔案則使用--bfile {your_file}。之後,可以新增所有其他必須的選項,例如,--assoc選項表示執行關聯分析,如圖2所示;此選項將告訴PLINK為每個SNP對感興趣的表型執行\(X^{2}\)測試。可以在單個命令列中結合使用多個選項。在PLINK中已經對命令選項實現了預設的順序,無論執行的命令列中的命令順序如何,它都能正常工作。一個有用的、有時也是強制性的選項是--out {outfile},它為PLINK提供了輸出檔案的名稱(PLINK會根據需要新增檔案字尾名)。請注意,PLINK會在靜默狀態下(不通知使用者)刪除任何與之同名的現有檔案(所以要注意當前資料夾下的現有檔案是否存在衝突)。請注意,PLINK中的命令選項超出了本文所討論的範圍;關於全部的可用選項,請參見 http://zzz.bwh.harvard.edu/plink/

圖2 PLINK命令列的結構。*並非所有shell都像這樣顯示。**如果PLINK不在當前目錄中,請提供安裝PLINK的目錄的路徑(例如,\usr\local\bin\plink)。請注意,此示例命令是使用Putty(一個免費的SSH和Telnet客戶端)生成的。使用其他資源(ssh連線軟體或者終端軟體)進行操作時,可能在圖形介面上有略微的變化;但是,PLINK命令的基本結構將是相同的。

3 | 基因(遺傳?)資料的質量控制

任何GWAS中都應該存在的一個重要的步驟就是使用適當的質量控制(Quality Control),如果沒有廣泛的質量控制,GWAS產生的結果將不會是可靠的,因為原始的基因型資料本質上是不完美的。資料中的錯誤可能有多種原因,例如,由於DNA樣本的質量很差、DNA與晶片(array)的雜交效果不佳、基因型探針效能不佳以及樣本的混淆或汙染(contamination)。例如,因為未能徹底控制這些資料上的問題,導致Sebastiani等人發表的一篇文章被撤回(2010) in Science (Sebastiani 等人, 2010, 2011; Sebastiani等人, 2012; Sebastiani等人, 2013)。被撤回的文章的結果受到了Illumina 610晶片技術錯誤和質量控制不足的影響。儘管經過適當的質量控制後,(文章中的)主要科學發現仍然能得到足夠的支援,但新的分析得到的結果偏差很大,以至於作者決定撤回該文章。

3.1 | 使用HapMap資料進行資料模擬

為了能夠使用真實的遺傳資料說明所有分析步驟,我們使用了來自國際HapMap專案(http://hapmap.ncbi.nlm.nih.gov/downloads/genotypes/2010-05_phaseIII/plink_format/ ;Gibbs 等人,2003 年)的公開資料模擬了一個數據集(N=207),其中包含了二進位制形式的結果度量。在本教程中,為了建立種族同質的資料集,(在資料集中)我們只包括了來自北歐和西歐(CEU)血統的猶他州(Utah)居民。由於HapMap資料的樣本量相對較小,這些模擬中的遺傳效應大小設定為比通常在複雜性狀遺傳研究中觀察到的值要大。需要注意到的是,檢測複雜性狀的遺傳風險因素需要更大的樣本量(例如,至少在數千,甚至可能是數萬或數十萬)。可以在 https://github.com/MareesAT/GWA_tutorial/ (1_QC_GWAS.zip) 找到具有模擬表型特徵的 HapMap 資料。

3.2 | 質量控制步驟概述

由於GWAS所面臨的挑戰,我們旨在說明基本的質量控制步驟並提供示例的指令碼。表1總結了七個質量控制的步驟,幷包括有關特定閥值的建議。但是,閥值可能會根據研究的具體特徵而有所不同。七個質量控制的步驟根據以下的內容來過濾掉SNP和個體:(1)個體和SNP缺失,(2)受試者的分配和遺傳性別不一致(見性別差異,sex discrepancy),(3)次要等位基因頻率(MAP,minor allele frequency),(4)與哈代-溫伯格平衡(Hardy-Weinberg equilibrium, HWE)的偏差,(5)雜合率,(6)相關性,以及(7)種族異常值(見種群分層,population stratification)。

通過遵循我們的線上教程(https://github.com/MareesAT/GWA_tutorial/ 上的1_QC_GWAS.zip + 2_Population_stratification.zip)中概述的所有步驟,可以獲得有關質量控制步驟1-7效能表現的實踐經驗。教程中提供了用於資料質量控制和潛在偏差來源視覺化的指令碼。這些指令碼對HapMap資料的CEU組執行了質量控制,可以應用於其他資料集,不過基於家庭的資料集和涉及多個不同種族的資料集除外。通常來說,如果樣本包括多個種族(例如,非洲人、亞洲人和歐洲人),建議分別對每個種族進行關聯測試並使用適當的方法,例如meta分析(Willer等人,2010),來將結果結合起來。如果您的樣本中包括來自單一族群的受試者,則可以通過以下討論的方法來校正剩餘的種群分層。

表1 在遺傳關聯分析之前應進行的七個質量控制步驟

步驟 命令 功能 閾值和相關解釋
1:個體和SNP缺失 (1)--geno

(2)--mind
(1)排除大部分受試者中缺失的SNP。在這一步中,刪除具有低基因型呼叫(low genotype calls)的SNP。
(2)排除基因型缺失率高的個體。在此步驟中,刪除具有低基因型呼叫的個體。
我們建議首先根據寬鬆的閾值(0.2;>20%)過濾SNP和個體,因為這將過濾掉SNP和具有非常高水平缺失的個體。然後可以應用具有更嚴格閾值的過濾器(0.02)。
注意,SNP過濾應該在個體過濾之前進行。
2:性別差異 -check-sex 根據X染色體雜合率/純合率檢查資料集中,記錄的個體的性別與其性別之間的差異。 可以揭示樣品混淆。如果許多受試者都有這種差異,則應仔細檢查資料。男性的X染色體純合度估計值應大於0.8,而女性的值應小於0.2。
3:次要等位基因頻率(MAF) --maf 僅包括高於設定的MAF閾值的SNP。 具有低MAF的SNP很少見,因此缺乏檢測SNP表型關聯的能力。
這些SNP也更容易出現基因分型錯誤。
MAF閾值應取決於您的樣本大小,較大的樣本可以使用較低的MAF閾值。
對於大樣本(N=100,000)與中等樣本(N=10000),通常使用0.01和0.05作為MAF閾值。
4:哈代-溫伯格平衡(HWE) --hwe 排除偏離哈代-溫伯格平衡的標記。 基因分型錯誤的常見指標,也可能表明進化選擇。
對於二元性狀,我們建議的用於排除的閾值:在cases中的HWE p值小於1e-10和在controls中p值小於1e-6。不太嚴格的case閾值可避免在選擇時丟棄與疾病相關的SNP(參見我們的線上教程)。
對於數量性狀,我們建議HWE的p值小於1e-6。
5:雜合率 見線上教程中的示例指令碼 排除雜合率較高或較低的個體。 (雜合率的)偏差(deviation)可以揭示樣品汙染和近親繁殖的存在。
我們建議刪除偏離樣本雜合率平均值±3 SD的個體。
6:相關性 (1)--genome

(2)--min
(1)計算所有樣本對的IBD(identity of descent)。
(2)設定閾值,並建立相關性高於所選閾值的個體列表。這意味著可以檢測到相關物件(個體),例如pi-hat>0.2(即二級親屬,second degree relatives)
在此分析中使用獨立的SNP(修剪,pruning)並將其限制為僅常染色體。
隱藏的相關性會干擾關聯分析。如果您有一個基於家庭的樣本(例如,父母-後代),則不需要刪除相關對,但統計分析應考慮家庭相關性。然而,對於基於種群的樣本,我們建議使用0.2的pi-hat閾值,這與文獻中提到的一致(Anderson 等人,2010 年;Guo 等人,2014 年)。
7:種群分層 (1)--genome
(2)--cluster --mds-plot k
(1)計算所有樣本對的IBD(identity of descent)。
(2)基於IBS生成資料中任何子結構的k維表示。
在此分析中使用獨立的SNP(修剪,pruning)並將其限制為僅常染色體。
K是維數,需要定義(通常是10)。這是有多個程式組成的質量控制中的重要步驟,但出於完整性考慮,我們在表中簡要地提及此步驟。此步驟將在“控制種群分層”一節中更詳細地描述。