GWAS 總體流程理解版

自己找了一些文章和視頻挪鹏，先總結(jié)了一部分复局，后面再做補充和實操

一. 相關(guān)概念理解

（1）GWAS：

全稱“全基因組關(guān)聯(lián)分析”嗜憔，使用統(tǒng)計模型找到與性狀關(guān)聯(lián)的位點涣旨，用于分子標(biāo)記選擇（MAS）或者基因定位

（2）GWAS分析的兩類性狀：

分類性狀（閾值性狀，質(zhì)量性狀）：比如抗病性蹂喻，顏色等等

質(zhì)量性狀指相對性狀的變異呈不連續(xù)性葱椭，呈現(xiàn)質(zhì)的中斷性變化的性狀。由1對或少數(shù)幾對主基因控制口四。如雞羽的蘆花斑紋和非蘆花斑紋孵运、角的有無、毛色蔓彩、血型等都屬于質(zhì)量性狀治笨。

連續(xù)性狀（數(shù)量性狀）：比如株高，體重赤嚼，產(chǎn)量等等（一般是呈現(xiàn)正態(tài)分布）

數(shù)量性狀指相對性狀的變異呈連續(xù)性旷赖，個體之間的差異不明顯，很難明確分組更卒。受微效多基因控制等孵，控制數(shù)量性狀的基因稱為數(shù)量性狀位點(quantitative trait loci, QTLs)。在 QTLs 中, 基因的效應(yīng)也有大有小蹂空。其中, 效應(yīng)較大的稱為主效QTL, 效應(yīng)較小的稱為微效QTL(或微效多基因)俯萌。動植物的許多重要經(jīng)濟性狀都是數(shù)量性狀，如作物的產(chǎn)量上枕、成熟期咐熙，奶牛的泌乳量，棉花的纖維長度辨萍、細度等等棋恼。

（3）GWAS的分析方法：

分類性狀：logistic回歸模型等等
連續(xù)性狀：GLM，MLM模型等等

二锈玉、準備工作

（1）準備文件

全基因組參考序列
全基因組注釋文件
樣本重測序文件（雙端測序）200個樣本左右或以上

（2）各類軟件和包

性狀分析
R 和 Rstudio
psych R包
lme4 R包
pheatmap R包
reshape2 R包
CMplot R包（繪制 snp 密度圖）

處理初始數(shù)據(jù)
fstqc （質(zhì)控）
bowtie2 （構(gòu)建基因組 index 爪飘，及序列比對）
samtools （構(gòu)建 samtools 基因組 index ， sam , bam 文件轉(zhuǎn)換）
bcftools （生成 vcf 原始文件）

標(biāo)記過濾
vcftools linux版
plink linux版
Tassel linux版

三拉背、實驗流程

（1）表型數(shù)據(jù)分析處理

將得到的表型數(shù)據(jù)（一般是數(shù)量性狀數(shù)據(jù)）進行分析處理

對數(shù)據(jù)進行描述性統(tǒng)計悦施，繪制直方圖觀察數(shù)據(jù)是否存在不合適的樣本數(shù)據(jù)

檢測正態(tài)性

剔除不合適的樣本

（2）原始數(shù)據(jù)處理（識別樣本中 snp 位點）

根據(jù)全基因組參考序列使用 bowtie2 建立 index 進行比對準備

使用 fastqc 對樣本重測序文件進行測序質(zhì)量檢測

使用 fastx_tolltik（或者 trimmomatic ）去除不良序列

使用 cutadaptor 去掉 read 兩端的 adaptor

使用 bowtie2 比對生成 .sam文件

使用 samtools view 將 sam 文件轉(zhuǎn)化為 bam文件

使用 samtools sort 將 .bam 文件進行排序

在 java 環(huán)境下使用 picard MarkDuplicates 去除 PCR 冗余（超聲波打斷的叫rad需要用Picard處理，酶切的叫ddrad不用處理可以省略這一步）

使用 bcftools mpileup 獲得原始 vcf 文件（這里相當(dāng)于將所有的序列文件進行合并去团，里面就含有 snp 的信息）

補充抡诞，在做完這些后穷蛹，可以在 R 中 CMplot 包繪制繪制SNP密度圖

（3）基因型過濾（網(wǎng)上教程大多數(shù)是從這一步開始，有 vcf 文件昼汗，或者將 vcf 文件轉(zhuǎn)換后的 plink 格式的文件肴熏，bed,bim,fam 文件）

如果是vfc文件

使用 vcftools 過濾（關(guān)于其參數(shù)設(shè)置的問題有待研究，一般是過濾掉低于缺失率高于50%的位點顷窒，次等位基因深度低于3蛙吏。在實際中，要更嚴格鞋吉，篩除第二等位基因率（次等位基因）頻率小于5%（在國際人類基因組單體型圖計劃（HapMap）中鸦做，MAF大于0.05 （5%）的SNP都被作為了調(diào)查目標(biāo)），缺失率大于10%,等位基因的個數(shù)要有兩個——這個可以調(diào)整）
使用 vcftools 將過濾后的文件轉(zhuǎn)換為 plink 格式的文件（或者使用 plink 也可以直接轉(zhuǎn)換）谓着，主要得到的是 bed 文件泼诱。

plink格式的文件主要有兩組五種
ped 和 map 是一組的
bed fam bim 是一組的
兩組可通過 plink -recode 參數(shù)相互轉(zhuǎn)換

轉(zhuǎn)換后可以使用plink再次過濾（對于計算不同的東西，如群體結(jié)構(gòu)Structure要求位點要少一些篩選的條件就不一樣）

（4）群體結(jié)構(gòu)分析

分析之前需要進行第三步的標(biāo)記篩選赊锚，然后根據(jù)以下條件去再次篩選：（一般只要按照LD去篩選就可以）

一定的物理距離取一個標(biāo)記作為代表進行分析
全基因組上抽取一部分標(biāo)記進行群體結(jié)構(gòu)的分析
按 LD 篩選治筒，LD強度大于一定閾值的標(biāo)記只保留其中一個用于分析

數(shù)據(jù)過濾，使用 plink 進行缺失和 maf 篩選
LD 篩選使用 plink 按照 LD 進行篩選
格式轉(zhuǎn)換舷蒲，然后使用 recode 參數(shù)進行轉(zhuǎn)換并得到 str 相關(guān)矩陣文件（后續(xù)就用該文件進行群體結(jié)構(gòu)分析）（可以根據(jù)需求轉(zhuǎn)換成 structure 或者 admixture 格式耸袜，structure比較麻煩一些）
利用得出的structure 或者 admixture 格式文件計算出最適 K 值，并在 R 中繪制不同 K 值時最低交叉驗證誤差的變化
繪制 structure 結(jié)構(gòu)圖（有些文章把這個省略掉了牲平，根據(jù)文章的側(cè)重不同可選擇保留）
PCA分析堤框，計算 PCA 矩陣（還可以通過EIGENSTRAT軟件計算主成分，計算各個主成分是否有顯著的統(tǒng)計學(xué)意義纵柿，將P值小于0.05的主成分計算在內(nèi)）蜈抓，然后繪制 PCA 圖

（5）親緣關(guān)系分析

可用于親緣關(guān)系分析的工具有很多，如：GCTA藐窄，LDAK，SPAGeDi酬土，EIGENSOFT荆忍，TASSEL，現(xiàn)在使用 TASSEL 比較多

GCTA撤缴，LDAK 常用于 snp 遺傳力估計或者性狀遺傳力的估計

同樣需要前期使用 plink 進行合理篩選
使用相應(yīng)軟件進行親緣關(guān)系矩陣計算（TASSEL 可以使用界面版也可以使用命令行版本）
計算PCA矩陣（還可以通過EIGENSTRAT軟件計算主成分刹枉，計算各個主成分是否有顯著的統(tǒng)計學(xué)意義，將P值小于0.05的主成分計算在內(nèi)）屈呕，繪制 PCA 圖
結(jié)果可視化（ kinship 值的分布圖和 kinship 熱圖）

（6）關(guān)聯(lián)分析

使用 tassel 進行 GLM/MLM/CMLM 分析(分為界面版和 Linux 版本微宝，界面版操作比較方便，但是用慣了 Linux 系統(tǒng)的肯定不會選界面版)
（這里要根據(jù)實驗?zāi)康幕⒄＃热珞w色蟋软，生長镶摘，低氧等）

界面版可以參考 tassel 使用說明書，下面主要講 Linux 操作

首先要對 vcf 文件進行排序岳守，這里用到的是一個 perl 腳本凄敢，不排序后面操作會報錯

轉(zhuǎn)換成 hapmp 格式，也可以不轉(zhuǎn)換直接操作湿痢，注意后面的參數(shù)設(shè)置就行

進行關(guān)聯(lián)分析（ GLM 和 MLM ）

對 tassel 計算的 GLM 或 MLM 結(jié)果進行校正涝缝，校正方法 Bonferroni 和 FDR （前者比較絕對，但篩選的結(jié)果一定是正確的譬重，后者可能會保留一些有意義的實驗結(jié)果拒逮，看情況使用）

FDR 校正，在 R 中使用 p.adjust 函數(shù)進行
Bonferroni 校正比 FDR 嚴格臀规，得到的有效 SNP 位點會少一些
這里可以參考之前我寫的關(guān)于兩者的區(qū)別

使用 CMplot 包進行結(jié)果可視化滩援，曼哈頓圖，QQplot（QQplot應(yīng)該是在校正前做出來還是校正后以现？）
篩選有意義的 SNP 位點（包括潛在有意義的位點）
計算膨脹系數(shù)lambda

基因組膨脹因子λ定義為經(jīng)驗觀察到的檢驗統(tǒng)計分布與預(yù)期中位數(shù)的中值之比狠怨，從而量化了因大量膨脹而造成結(jié)果的假陽性率。換句話說邑遏，λ定義為得到的卡方檢驗統(tǒng)計量的中值除以卡方分布的預(yù)期中值佣赖。預(yù)期的P值膨脹系數(shù)為1，當(dāng)實際膨脹系數(shù)越偏離1记盒，說明存在群體分層的現(xiàn)象越嚴重憎蛤，容易有假陽性結(jié)果，需要重新矯正群體分層纪吮。

（7）根據(jù) GWAS 結(jié)果找到 QTL區(qū)域（這個后面的操作就了解的比較少了,后面學(xué)到了會再補充）

利用 R/qtl 軟件 MapQTL 軟件等定位 QTL
根據(jù) QTL 定位找到相關(guān)性狀的基因（這個是用什么軟件俩檬？）
根據(jù)基因的位置和功能來反向驗證結(jié)果（這里是不是要用富集分析之類的？）
后面還有一連串對 QTL 的分析

（8）驗證實驗

驗證實驗也有很多種碾盟，敲除棚辽，抑制基因表達，定量PCR等