GWAS--Plink分析

時間 2021-08-13

標籤 git github 編程 vim bash 函數工具學習測試 spa 欄目 Git 简体版

原文原文鏈接

概述：

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/

❞

這個教程很是的經典，我看網上不少人推薦。

❝

相關的文章： https://onlinelibrary.wiley.com/doi/full/10.1002/mpr.1608

❞

教程中包括數據的過濾，SNP的過濾，樣本的過濾，質控的標準等等，介紹的很是清楚，看完這篇文章，感受plink的語法知識又增長了不少。

下載R語言和plink軟件

R： https://www.r-project.org/
plink： http://zzz.bwh.harvard.edu/plink/ https://www.cog-genomics.org/plink2

plink的調用

打開terminal，cd到解壓到的文件夾（把plink拖入終端便可）
以後每次使用須要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值越小，說明該位點越不符合哈溫平衡。

❞

「兩個目的：」