通過前面的幾篇推文，相信大家對GWAS的質(zhì)控和分析原理已經(jīng)初步了解。接下來，我們需要將之前學(xué)習(xí)到的知識運(yùn)用到實(shí)踐中。為了更好的進(jìn)行實(shí)踐，選擇一個(gè)好的參考數(shù)據(jù)集就就成為了一個(gè)必不可少的環(huán)節(jié)。千人基因組數(shù)據(jù)集，作為一個(gè)旨在繪制最詳盡的、最有醫(yī)學(xué)應(yīng)用價(jià)值的人類基因組遺傳多態(tài)性圖譜數(shù)據(jù)集，就成為了我們這次分析的首選參考數(shù)據(jù)集。本篇推文，將通過對千人基因組數(shù)據(jù)集進(jìn)行深入而細(xì)致的GWAS分析，挖掘出不同種屬之間的進(jìn)化關(guān)系。那么，準(zhǔn)備好了么，讓我們開始GWAS應(yīng)用分析之旅?。?！

一?GWAS分析之旅

1.1?前人基因組的下載與數(shù)據(jù)整理

實(shí)際上，本教程并未下載人體所有染色體上的snp數(shù)據(jù)，僅下載了第21號染色體上的snp數(shù)據(jù)。這是因?yàn)槿绻腥旧w上的snp數(shù)據(jù)將近有T級別的數(shù)據(jù)，下載和分析都需要耗費(fèi)大量時(shí)間。但是，我們分析所用21號染色體數(shù)據(jù)僅8個(gè)G的數(shù)據(jù)左右，大大加快了分析相關(guān)的進(jìn)程。這也對應(yīng)之前的推文中<GWAS分析<二>之?dāng)?shù)據(jù)質(zhì)控的原理>所提到的如何大家如何使用最少的數(shù)據(jù)完成本次項(xiàng)目的分析的策略。

下文的兩行代碼分別為將千人基因組文件（vcf）轉(zhuǎn)為plink格式（bed、fam、bim），并為缺失rs標(biāo)識符的SNP分配唯一的標(biāo)識符。雖然，對于千人基因組而言，數(shù)據(jù)中缺失的rs標(biāo)識符不是問題。

1.2?移除不符合標(biāo)準(zhǔn)的snp類型

分別通過‘--geno 0.2’、‘--mind 0.2’、‘--geno 0.02’、‘--mind 0.02’和‘--maf 0.05’選項(xiàng)完成對符合標(biāo)準(zhǔn)的snp進(jìn)行過濾。過濾的標(biāo)準(zhǔn)詳見之間的推文。

1.3?提取snp變異

一方面，本教程將從千人基因組數(shù)據(jù)集中提取HapMap數(shù)據(jù)集中存在的snp變異，另一方面，本教程也將從HapMap數(shù)據(jù)集中提取千人基因組數(shù)據(jù)集中存在的snp變異。

1.4 hapmap文件的構(gòu)建

利用上一步提取的HapMap_MDS.map文件構(gòu)建hapmap.txt，該文件每行包含一個(gè)SNP id和物理位置。隨后，利用plink軟件將hapmap.txt中snp信息更新1kG_MDS6中，輸出到1kG_MDS7文件中。

1.5?等位基因的校準(zhǔn)

本步驟是確保HapMap和1000基因組項(xiàng)目數(shù)據(jù)集中的參考等位基因被重新指定或調(diào)整。

1.6?鏈校正問題

本步驟將檢查潛在的鏈校正問題，并為不應(yīng)該存在的A1等位基因分配生成警告輸出。

1.7 snp翻轉(zhuǎn)

本步驟翻轉(zhuǎn)（正義或反義鏈上）SNP以解決鏈校正問題。本步驟將輸出SNP標(biāo)識符并刪除重復(fù)項(xiàng)，再生成一個(gè)包含812個(gè)SNP的文件。文件中的這些snp來自于兩個(gè)文件之間不對應(yīng)的SNP。

1.8 snp再次檢查

檢查翻轉(zhuǎn)后仍有問題的SNP。

1.9?去除有問題的snp

本步驟將從HapMap和1000個(gè)基因組中刪除有問題的SNP。上一步驟生成了一個(gè)42個(gè)SNP的列表，這些SNP在翻轉(zhuǎn)和參考等位基因被重新指定或調(diào)整后，在HapMap和1000基因組數(shù)據(jù)集之間造成了84處差異。因此，本步驟將分別從兩個(gè)數(shù)據(jù)集中刪除42個(gè)有問題的SNPs，并合并HapMap與千人基因組數(shù)據(jù)。值得注意，雖然HapMap和千人基因組的數(shù)據(jù)集之間存在樣本重疊，但是這對本次分析并不重要。在由千人基因組數(shù)據(jù)上錨定的HapMap-CEU數(shù)據(jù)上執(zhí)行MDS。

1.10?種群異常個(gè)體去除

本步驟將人口標(biāo)識轉(zhuǎn)換為更高級的人口標(biāo)識（即AFR、AMR、ASN和EUR），并創(chuàng)建一個(gè)自己的數(shù)據(jù)文件。通過對結(jié)果進(jìn)行可視化，顯示我們的數(shù)據(jù)屬于歐洲1000個(gè)基因組的數(shù)據(jù)組（輸出文件是MDS.pdf）。并且，出于教程訓(xùn)練的目的，本步驟給出腳本也將過濾出人口分層異常值，從而為下一步驟生成適當(dāng)?shù)奈募?，用以排除種群中異常的個(gè)體。本步驟將在HapMap數(shù)據(jù)中選擇低于過濾閾值的個(gè)體。需要注意的是，異常水平不是固定的閾值，而是基于前兩個(gè)維度的可視化來確定（讀者可自行確定）。

1.11 plink協(xié)變量文件創(chuàng)建

本步驟將在單倍體圖譜數(shù)據(jù)中提取這些個(gè)體（排除異常值的個(gè)體）。注意，由于我們的單體型圖數(shù)據(jù)確實(shí)包含任何種族異常值，因此在這一步中沒有刪除任何個(gè)體（如果我們的數(shù)據(jù)包括超出我們設(shè)定的閾值的個(gè)人，那么這些個(gè)體將被刪除）。隨后，我們將基于MDS創(chuàng)建協(xié)變量。注意，僅對沒有種族異常值的HapMap數(shù)據(jù)執(zhí)行MDS分析。更改文件的格式。將mds文件轉(zhuǎn)換為plink協(xié)變量文件。注意：10個(gè)MDS維度的值隨后被用作關(guān)聯(lián)分析中的協(xié)變量。

進(jìn)行到這里，本公眾號也要恭喜讀者，實(shí)現(xiàn)了人口分層數(shù)據(jù)的質(zhì)控工作?。。?/p>

1.12?關(guān)聯(lián)分析

對于關(guān)聯(lián)分析，本步驟使用的文件來自于上一步驟中（人口分層）中生成的文件（HapMap_3_r3_13.bed、HapMap_3_r3_13.bim、HapMap_3_r3_13.fam和covar_mds.txt）。需要注意的是，本教程中的--assoc選項(xiàng)不允許校正主成分（PC）/MDS成分等協(xié)變量，這使得它不太適合關(guān)聯(lián)分析。因此，本步驟將使用10個(gè)主成分（covar_MDS. txt）作為協(xié)變量進(jìn)行l(wèi)ogistic分析。

本步驟使用ide-covar選項(xiàng)僅在輸出文件中顯示SNP的相加結(jié)果。本步驟也將刪除NA值，因?yàn)檫@些值會在后續(xù)生成圖片時(shí)引發(fā)問題。需要注意的是，本步驟GWAS分析中所產(chǎn)生的結(jié)果將在隨后的步驟用于可視化分析，顯示數(shù)據(jù)集中是否包含基因組范圍內(nèi)的重要SNP。另外，如果是定量結(jié)果測定，則本步驟的選項(xiàng)--logistic應(yīng)替換為--linear。使用--assoc選項(xiàng)也可以用于定量結(jié)果測量。除了常規(guī)的全基因組顯著性閾值5.0E-8之外，還有多種方法可以處理多重檢測。

1.13 p值計(jì)算

上一步給出了Bonferroni校正的p值，以及FDR和其他文件。本步是一個(gè)計(jì)算密集的步驟。本教程相對應(yīng)的文章中描述了這種方法的進(jìn)一步優(yōu)缺點(diǎn)，它可以用于關(guān)聯(lián)和處理多個(gè)測試(https://www.ncbi.nlm.nih.gov/pubmed/29484742)為了減少計(jì)算時(shí)間，我們只對22號染色體的一部分SNP進(jìn)行了測試。EMP2 列為多次測試提供了校正的p值。根據(jù)上述步驟生成的SNP子集篩選您的bfile。執(zhí)行1000000次測試。按從最低到最高的p值順序排列數(shù)據(jù)。檢查有序排列結(jié)果

1.14?曼哈頓圖和QQ圖的繪制

本步驟的話主要是通過可視化腳本繪制曼哈頓圖和QQ圖。

二?慣例小結(jié)

恭喜讀者，到了這里，我們共同學(xué)習(xí)了如何對GWAS數(shù)據(jù)進(jìn)行質(zhì)控以及如何運(yùn)用GWAS進(jìn)行群體分析。事實(shí)上，就如同開篇所提到的那樣，GWAS的分析遠(yuǎn)遠(yuǎn)不僅如此，他還可以與其他分析進(jìn)行聯(lián)合，進(jìn)一步挖掘信息。如近期發(fā)表在Nature medicine上的一篇文章就利用MR整合GWAS和eQTL發(fā)現(xiàn)復(fù)雜疾病的causal genes（具體內(nèi)容可見文獻(xiàn)Actionable druggable genome-wideMendelian randomization identifies repurposing opportunities for COVID-19）。因此，讀者可以在GWAS分析的基礎(chǔ)上結(jié)合其他組學(xué)或者實(shí)驗(yàn)，實(shí)現(xiàn)文章在質(zhì)上的飛躍。也希望各位讀者天天有大paper。

本公眾號開發(fā)的相關(guān)軟件，Multi-omics Hammer軟件和Multi-omics Visual軟件歡迎大家使用。文末是本公眾號在其他平臺的賬戶，也歡迎大家關(guān)注并多提意見。

簡書：WJ的生信小院

公眾號：生信小院

博客園：生信小院

最后，也歡迎各位大佬能夠在本平臺上：1傳播和講解自己發(fā)表的論文；2：發(fā)表對某一科研領(lǐng)域的看法；3：想要達(dá)成的合作或者相應(yīng)的招聘信息；4：展示自己以尋找博后工作或者博士就讀的機(jī)會;5:博導(dǎo)提供博后工作或者博士攻讀機(jī)會，都可以后臺給筆者留言。希望本平臺在進(jìn)行生信知識分享的同時(shí)，能夠成為生信分析者的交流平臺，能夠?qū)崿F(xiàn)相應(yīng)的利益互補(bǔ)和雙贏（不一定能實(shí)現(xiàn)，但是夢想總得是有的吧）。

另外，怎么說呢，投幣也可，不強(qiáng)求，但奢求。