GWAS--Plink分析

概述:

 GWAS全稱「全基因組關聯分析」,使用統計模型找到與性狀關聯的位點,用於分子標記選擇(MAS)或者基因定位。全基因組關聯分析是對多個個體在全基因組範圍的遺傳變異多態性進行檢測,得到基因型,進而將基因型與可觀測的性狀,即表型,進行羣體水平的統計學分析,根據統計量或P值篩選出最有可能影響該性狀的遺傳變異。此次學習的教程是Plink作GWAS。git

PLINK是一個免費的開源全基因組關聯分析工具集,在SNP數據統計,過濾,GWAS分析中均可以用得上,並且計算速度很快。直接去百度上搜索plink就能夠很容易就找到plink官網( http://www.cog-genomics.org/plink2)功能大概有如下幾種:
  • 數據管理: SNP數據格式的轉換,合併兩個或多個文件,提取SNP子集,以二進制文件格式壓縮數據等。
  • 質量控制的SNP數據統計: 計算丟失基因型率,等位基因,基因型頻率,HWE測試,個體和個體對的近親繁殖,IBS和IBD統計,LD區域計算等。
  • GWAS關聯分析
  • Meta分析

GWAS 操做流程1-1:數據下載和plink配置

GWAS分析的兩類性狀:
  • 分類性狀(閾值性狀,質量性狀):好比抗病性,顏色等等
  • 連續性狀(數量性狀):好比株高,體重,產量等等
GWAS的分析方法:
  • 分類性狀:logistic等等
  • 連續性狀:GLM,MLM模型等等
「通常線性模型(GLM):
這裏,SNP做爲固定因子,能夠考慮其它協變量(好比性別,PCA,羣體結構等等)
「混合線性模型(MLM):
  • 固定因子:SNP + 能夠考慮其它協變量(好比性別,PCA,羣體結構等等),這裏固定因子和前面的GLM同樣
  • 隨機因子:親緣關係矩陣(K矩陣或者A矩陣)

參考:

教程代碼和數據下載: https://github.com/MareesAT/GWA_tutorial/
這個教程很是的經典,我看網上不少人推薦。
教程中包括數據的過濾,SNP的過濾,樣本的過濾,質控的標準等等,介紹的很是清楚,看完這篇文章,感受plink的語法知識又增長了不少。
 
 
下載R語言和plink軟件
 
 
plink的調用
  1. 打開terminal,cd到解壓到的文件夾(把plink拖入終端便可)
  2. 以後每次使用須要cd到此路徑,輸入./plink再輸入參數,如:
./plink --help
./plink --vcf Arab.vcf --allow-extra-chr --maf 0.05 --recode --out Arab
 
 
若需在終端中直接調用plink,則須要把plink寫入全局路徑
 
  • vim .bash_profile
  • #而後按i鍵進入編輯模式,輸入如下語句,其中路徑爲解壓plink的文件夾路徑
  •  
  • export PATH=/Users/Downloads/plink_mac_20200428:$PATH
  •  
  • #輸完後回車,按ESC
  • :w #表示保存
  • :q #表示退出
  • source .bash_profile #更新便可 

GWAS 操做流程2-1:缺失質控

 
「主要介紹」
--geno篩選個體;--mind篩選SNP
GWAS分析時,拿到基因型數據,拿到表型數據,要首先作如下幾點:
  • 1,查看本身的表型數據,是否有問題
  • 2,查看本身的基因型數據,是否有問題
而後再進行建模,獲得顯著性SNP以及可視化結果。
「爲什麼對缺失數據進行篩選?」
不管是測序仍是芯片,獲得的基因型數據要進行質控,而對缺失數據進行篩選,能夠去掉低質量的數據。若是一個個體,共有50萬SNP數據,發現20%的SNP數據(10萬)都缺失,那這個個體咱們認爲質量不合格,若是加入分析中可能會對結果產生負面的影響,因此咱們能夠把它刪除。一樣的道理,若是某個SNP,在500個樣本中,缺失率爲20%(即該SNP在100個個體中都沒有分型結果),咱們也能夠認爲該SNP質量較差,將去刪除。固然,這裏的20%是過濾標準,能夠改變質控標準。下文中的質控標準是2%。
 

1. plink數據格式轉化

數據使用上一篇的數據,由於數據是plink的bfile格式,二進制不方便查看,咱們將其轉化爲文本map和ped的格式。
plink --bfile HapMap_3_r3_1 --recode --out test
結果生成:test.map  test.ped

2. 查看基因型個體和SNP數量

wc -l test.map test.ped
能夠看出,共有165個基因型個體,共有1447897個SNP數據。
 
「預覽一下ped文件:」
 
「預覽一下map文件:」

3. 查看一下個體缺失的位點數,每一個SNP缺失的個體數

看一下描述:
--missing: Sample missing data report written to plink.imiss, and variant-based
missing data report written to plink.lmiss.
結果生成兩個文件,分別是一個個體ID上SNP缺失的信息,另外一個是每一個SNP在個體ID中缺失的信息。
  • 個體缺失位點的統計在plink.imiss中
  • 單個SNP缺失的個體數在plink.lmiss.中
「個體缺失位點統計預覽:」第一列爲家系ID,第二列爲個體ID,第三列是否表型缺失,第四列缺失的SNP個數,第五列總SNP個數,第六列缺失率。
「SNP缺失的個體數文件預覽:」第一列爲染色體,第二列爲SNP名稱,第三列爲缺失個數,第四列爲總個數,第五列爲缺失率
 
 

4. 對個體及SNP缺失率進行篩選

  • 1, 若是一個SNP在個體中2%都是缺失的,那麼就刪掉該SNP,參數爲:--mind 0.02
  • 2,若是一個個體,有2%的SNP都是缺失的,那麼就刪掉該個體,參數爲:--geno 0.02

4.  1 對個體缺失率進行篩選

「先過濾個體缺失率高於2%的SNP」
plink --bfile HapMap_3_r3_1 --geno 0.02 --make-bed --out HapMap_3_r3_2
「轉化爲map和ped的形式:」
plink --bfile HapMap_3_r3_2 --recode --out test
查看一下過濾後的行數,
「以前的爲:」
1457897 test.map
165 test.ped
「如今的爲:」
1430443 test.map
165 test.ped
能夠看出,過濾了2萬多個位點。
從當時的log日誌裏也能夠看出這一點:
PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/
(

 

能夠看到--geno,過濾了27454個位點。

4.  2 對SNP缺失率進行篩選

「過濾SNP缺失率高於2%的個體」
plink --bfile HapMap_3_r3_2 --mind 0.02 --make-bed --out HapMap_3_r3_3
「查看日誌:」 

 

 「沒有過濾掉個體,剩餘:」個體:165 SNP:1430443github

4.  3  同時對個體和SNP的缺失率進行篩選

「兩步合在一塊兒,即過濾位點,又過濾個體:」
plink --bfile HapMap_3_r3_1 --geno 0.02 --mind 0.02 --make-bed --out HapMap_3_r3_5
plink日誌:
能夠看出,二者最終結果是同樣的。
 

GWAS 操做流程2-2:性別質控

人類性別的信息的質控,主要是根據性染色上SNP的比值,判斷性別,而後把性別錯誤的個體去掉或者更改性別信息。對其它物種參考意義不大,由於在動物中通常把性別信息的SNP去掉,植物中通常都是雌雄同體的,不涉及到這個問題,之因此會有這一篇,是由於原文中有這個信息,並且plink 也有--check-sex的參數,因此操做一下,留下筆記。
❞「原理:」檢查性別差別。先驗信息,女性的受試者的F值必須小於0.2,男性的受試者的F值必須大於0.8。這個F值是基於X染色體近交(純合子)估計。不符合這些要求的受試者被PLINK標記爲「PROBLEM」。
「上一步,去掉缺失信息後,如今有文件是過濾缺失後的文件:」
HapMap_3_r3_5.bed HapMap_3_r3_5.fam HapMap_3_r3_5.irem
HapMap_3_r3_5.bim HapMap_3_r3_5.hh HapMap_3_r3_5.log

1. 檢查性別衝突

plink --bfile HapMap_3_r3_5 --check-sex
結果文件:plink.sexcheck第一列爲家系ID,第二列爲個體ID,第三列爲系譜中的性別,第四列爲SNP推斷的性別,第五列是否正常,第六列爲F值。
 
 

2. 提取錯誤的ID

咱們使用grep過濾一下:根據STATUS列,若是有問題的話,爲「PROBLEM」,咱們能夠根據這個關鍵詞將有問題的行打印出來。
grep "PROBLEM" plink.sexcheck
1349 NA10854 2 1 PROBLEM 0.99
能夠看出,個體NA10854是有問題的。
將相關錯誤的ID提取出來(家系ID,個體ID),之因此提取家系ID和個體ID,由於plink有參數remove能夠根據ID進行篩選。
grep 'PROBLEM' plink.sexcheck | awk '{print $1,$2}' >sex_discrepancy.txt
咱們將結果保存在sex_discrepancy.txt。

3. 使用remove去掉個體

plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6
固然,你也能夠對個體進行斷定填充,這是用--impute-sex就能夠實現,這樣的話那個錯誤的個體會根據統計量更改性別信息。這裏咱們選擇的是刪掉這個個體。

4. 過濾的關鍵詞

去掉個體或者SNP,關鍵詞不同,容易混淆,這裏總結一下。
「保留或去掉個體:」
--keep <filename>
--remove <filename>
 
--keep-fam <filename>
--remove-fam <filename>
 
「保留或去掉SNP:」
--extract ['range'] <filename>
--exclude ['range'] <filename>
 

GWAS 操做流程2-3:MAF過濾

上一次咱們通過去掉缺失,去掉錯誤的性別信息,獲得的文件爲:
HapMap_3_r3_6.bed HapMap_3_r3_6.fam HapMap_3_r3_6.log
HapMap_3_r3_6.bim HapMap_3_r3_6.hh
這裏,咱們根據最小等位基因頻率(MAF)去篩選。
「爲何要根據MAF去篩選?」
最小等位基因頻率怎麼計算?好比一個位點有AA或者AT或者TT,那麼就能夠計算A的基因頻率和T的基因頻率,qA + qT = 1,這裏誰比較小,誰就是最小等位基因頻率,好比qA = 0.3, qT = 0.7, 那麼這個位點的MAF爲0.3. 之因此用這個過濾標準,是由於MAF若是很是小,好比低於0.02,那麼意味着大部分位點都是相同的基因型,這些位點貢獻的信息很是少,增長假陽性。更有甚者MAF爲0,那就是全部位點只有一種基因型,這些位點沒有貢獻信息,放在計算中增長計算量,沒有意義,因此要根據MAF進行過濾。

1. 去掉性染色體上的位點

「思路:」
  • 在map文件中選擇常染色體,提取snp信息
  • 根據snp信息進行提取
「提取常染色體上的位點名稱:」
由於這裏是人的數據,因此染色體只須要去1~22的常染色體,提取它的家系ID和個體ID,後面用於提取。
awk '{ if($1 >=1 && $1 <= 22) print $2}' HapMap_3_r3_6.bim > snp_1_22.txt
wc -l snp_1_22.txt
1399306 snp_1_22.txt
常染色體上共有1399306位點。

2. 提取常染色體上的位點

這裏,用到了位點提取參數--extract
plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7
 

3. 計算每一個SNP位點的基因頻率

首先,經過參數--freq,計算每一個SNP的MAF頻率,經過直方圖查看總體分佈。可視化會更加直接。
plink --bfile HapMap_3_r3_7 --freq --out MAF_check
結果文件:MAF_check.frq預覽:
 
 

4. 去掉MAF小於0.05的位點

plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8
 
日誌:
PLINK v1.90b6.5 64-bit (13 Sep 2018) www.cog-genomics.org/plink/1.9/

 

GWAS 操做流程2-4:哈溫平衡檢驗

「什麼是哈溫平衡?」
哈迪-溫伯格(Hardy-Weinberg)法則 哈迪-溫伯格(Hardy-Weinberg)法則是羣體遺傳中最重要的原理,它解釋了繁殖如何影響羣體的基因和基因型頻率。這個法則是用Hardy,G.H (英國數學家) 和Weinberg,W.(德國醫生)兩位學者的姓來命名的,他們於同一年(1908年)各自發現了這一法則。他們提出在一個不發生突變、遷移和選擇的無限大的隨機交配的羣體中,基因頻率和基因型頻率將逐代保持不變。---百度百科
「怎麼作哈溫平衡檢驗?」
「卡方適合性檢驗!」 ,一個羣體是否符合這種情況,即達到了遺傳平衡,也就是一對等位基因的3種基因型的比例分佈符合公式:p2+2pq+q2=1,p+q=1,(p+q)2=1.基因型MM的頻率爲p2,NN的頻率爲q2,MN的頻率爲2pq。MN:MN:NN=P2:2pq:q2。MN這對基因在羣體中達此狀態,就是達到了遺傳平衡。若是沒有達到這個狀態,就是一個遺傳不平衡的羣體。但隨着羣體中的隨機交配,將會保持這個基因頻率和基因型分佈比例,而較易達到遺傳平衡狀態。應用Hardy-Weinberg遺傳平衡吻合度檢驗方法,把計算獲得的基因頻率代入,計算基因型平衡頻率,再乘以總人數,求得預期值(e)。把觀察數(O)與預期值(e)做比較,進行χ2檢驗。病例組和對照組的基因型分佈的觀察值和預期值差別無顯著性(P>0.05),符合遺傳平衡定律.
「哈溫平衡過濾和MAF過濾的區別?」
以前,我對這兩個概念有點混淆,後來明白過來了。這兩個概念一個是對基因頻率進行的篩選,一個是對基因型頻率進行的篩選。對於一個位點「AA AT TT」,其中A的頻率爲基因頻率,AA爲基因型頻率。MAF直接是對基因頻率進行篩選,而哈溫平衡檢驗,則是根據基因型推斷出理想的(AA,AT,TT)的分佈,而後和實際觀察的進行適合性檢驗,而後獲得P值,根據P值進行篩選。即P值越小,說明該位點越不符合哈溫平衡。
「兩個目的:」
  • 計算全部位點的哈溫檢測結果
  • 刪除SNP中不符合哈溫平衡的位點

1. 計算全部位點的HWE的P值

plink --bfile HapMap_3_r3_8 --hardy
plink.hwe的數據格式:
  • CHR 染色體
  • SNP SNP的ID
  • TEST 類型
  • A1 minor 位點
  • A2 major 位點
  • GENO 基因型分佈:A1A1, A1A2, A2A2
  • O(HET) 觀測雜合度頻率
  • E(HET) 指望雜合度頻率
  • P 哈溫平衡的卡方檢驗P-value值
結果預覽:

2. 提取哈溫p值小於0.0001的位點

這裏咱們使用awk:
awk '{if($9 < 0.0001) print $0}' plink.hwe >plinkzoomhwe.hwe
共有123個位點,其中UNAFF爲45個位點。

3.  設定過濾標準1e-4

plink --bfile HapMap_3_r3_8 --hwe 1e-4 --make-bed --out HapMap_3_r3_9
日誌:

 

能夠看到,共有45個SNP根據哈溫的P值過濾掉了,和上面手動計算的同樣。
 
 
 

 GWAS 操做流程2-5:雜合率檢驗

通常天然羣體,基因型個體的雜合度太高或者太低,都不正常,咱們須要根據雜合度進行過濾。誤差可能代表樣品受到污染,近親繁殖。咱們建議刪除樣品雜合率平均值中偏離±3 SD的個體。
個人理解:非天然羣體中,好比自交系,雜交種F1,這些羣體不須要過濾雜合度。
「參數過濾和手動過濾」plink有個特色,全部的過濾標準,均可以生成過濾前的文件,而後能夠手動過濾,也能夠用參數進行過濾。
  • 好比:--missing生成結果,也能夠用--geno和--mind過濾。
  • 好比:--hardy生成結果,可使用--hwe過濾
  • 好比:--freq生成結果,能夠用--maf過濾 可是雜合度--het,沒有過濾的函數,只能經過編程去提取ID,而後用--remove去實現。

1. 計算雜合度

plink --bfile HapMap_3_r3_9 --het --out R_check
結果文件:
.het (method-of-moments F coefficient estimates)
Produced by --het.
 
A text file with a header line, and one line per sample with the following six fields:
 
FID Family ID # 家系ID
IID Within-family ID # 個體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 #
因此,雜合度 = (N-O)/N
 
 

GWAS 操做流程2-6:去掉親緣關係近的個體

「注意:」
這裏講親子關係的個體移除,不是必需要的,好比咱們分析的羣體裏面有親子關係的個體,想要進行分析,不須要作這一步的篩選。

1. 計算pihat > 0.2的組合

plink --bfile HapMap_3_r3_10 --genome --min 0.2 --out pihat_min0.2
說明文檔:
--genome invokes an IBS/IBD computation, and then writes a report with the following fields to plink.genome:
 
FID1 Family ID for first sample
IID1 Individual ID for first sample
FID2 Family ID for second sample
IID2 Individual ID for second sample
RT Relationship type inferred from .fam/.ped file
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)
PHE Pairwise phenotypic code (1, 0, -1 = AA, AU, and UU 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)
  

2. 提取Z1大於0.9的個體

awk '{if($8>0.9) print $0}' pihat_min0.2.genome > zoom_pihat.genome
過濾出91個組合:

 3. 刪除親子關係的個體

plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11
 
日誌:

能夠看出,51個個體被移除。編程

相關文章
相關標籤/搜索