概述：

GWAS全稱“全基因組關聯分析”，使用統計模型找到與性狀關聯的位點，用于分子標記選擇（MAS）或者基因定位，全基因組關聯分析是對多個個體在全基因組范圍的遺傳變異多型性進行檢測，獲得基因型，進而將基因型與可觀測的性狀，即表型，進行群體水平的統計學分析，根據統計量或P值篩選出最有可能影響該性狀的遺傳變異，這次學習的教程是Plink做GWAS，

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：1430443

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個個體被移除，

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/296297.html

標籤：其他

上一篇：最詳細且易懂的Ajax、promise、axios知識

下一篇：Kubernetes 安裝手冊（Ubuntu非高可用版）

GWAS--Plink分析

概述：