怎樣使用Python繪制GWAS分析中的曼哈頓圖和QQ圖

這篇文章給大家介紹怎樣使用Python繪制GWAS分析中的曼哈頓圖和QQ圖，內(nèi)容非常詳細(xì)，感興趣的小伙伴們可以參考借鑒，希望對(duì)大家能有所幫助。

曼哈頓圖和QQ圖是兩個(gè)在全基因組關(guān)聯(lián)（GWAS）分析里面最常出現(xiàn)的圖形，基本上已經(jīng)是GWAS的標(biāo)配，幾乎在每篇GWAS的文章都會(huì)見到，它們的作用和所要傳達(dá)出來的信息我也在上一篇關(guān)于GWAS的博文中做了些說明，在這里我們就只集中在如何用Python和geneview將其有效地展現(xiàn)出來。

首先，準(zhǔn)備一些數(shù)據(jù)來作為例子。

我這里用來展現(xiàn)的數(shù)據(jù)是2011年丹麥人所做過的一個(gè)關(guān)于年輕人過度肥胖的GWAS研究——GOYA，數(shù)據(jù)也是從他們所發(fā)表的結(jié)果中獲得，總共有5,373個(gè)樣本，其中超重的個(gè)體（case）有2,633個(gè)，正常的個(gè)體（control）是2,740個(gè)，從樣本量上看還算可以。為了方便使用，我對(duì)其做了相關(guān)的處理，包括從PED和MAP文件到GEN文件的生成，并重復(fù)了一次case-control的關(guān)聯(lián)性分析，計(jì)算出了芯片上所研究的各個(gè)SNP位點(diǎn)與肥胖相關(guān)的顯著性程度（即p-value），最后又將結(jié)果數(shù)據(jù)抽取出來做成數(shù)據(jù)集——放在這里供下載（15.6Mb，csv格式）。

【注】以上內(nèi)容雖提及到了一些領(lǐng)域內(nèi)術(shù)語和相關(guān)文件格式，但若不懂也請(qǐng)不必糾結(jié)，因?yàn)楹罄m(xù)處理都是基于這個(gè)最終的數(shù)據(jù)集來完成的。

接著，需要將geneview軟件包加入到你的Python中，有多種不同的方式，但推薦直接使用pip，以下是安裝比較穩(wěn)定的發(fā)布版，直接在終端命令行下(Linux or Mac)輸入：

pip install geneview

或者，也可以直接從github上安裝正在開發(fā)的版本：

pip install git+git://github.com/ShujiaHuang/geneview.git#egg=geneview

第三種辦法就是直接下載源碼，然后自行編譯，雖然不推薦這種做法（因?yàn)檫€有依賴包必須自行下載安裝，過程會(huì)比較麻煩低效），但對(duì)于某些不能連接外網(wǎng)的集群也只能如此，這三種方式都是可行的。

曼哈頓圖

將示例數(shù)據(jù)下載下來：

wget https://raw.githubusercontent.com/ShujiaHuang/geneview-data/master/GOYA.csv 

先簡單地查看一下數(shù)據(jù)的格式:

chrID,rsID,position,pvalue

1,rs3094315,742429,0.144586

1,rs3115860,743268,0.230022

1,rs12562034,758311,0.644366

1,rs12124819,766409,0.146269

1,rs4475691,836671,0.458197

1,rs28705211,890368,0.362731

1,rs13303118,908247,0.22912

1,rs9777703,918699,0.37948

1,rs3121567,933331,0.440824

一共是4列（逗號(hào)分隔），分別為：[1]染色體編號(hào)，[2]SNP rs 編號(hào)，[3] 位點(diǎn)在染色體上的位置，[4]顯著性差異程度（pvalue）。在本例曼哈頓圖中我們只需要使用第1,3和4列；而QQ圖則只需要第4列——pvalue。

下面先從繪制曼哈頓圖開始。我們先將需要的數(shù)據(jù)讀取到一個(gè)列表中，可以這樣做：

import csv

data = []

with open("GOYA.csv") as f:

f_csv = csv.reader(f)

headers = next(f_csv)

data = [[row[0], int(row[2]), float(row[3])] for row in f_csv]

現(xiàn)在GOYA.csv中的數(shù)據(jù)就都存放在data列表中了，由于Python在讀取文件中數(shù)據(jù)時(shí)，都是以string類型存放，因此對(duì)于第3和第4列的數(shù)據(jù)有必要事先把做點(diǎn)類型轉(zhuǎn)換。

接下來，調(diào)用geneview中的曼哈頓圖函數(shù)。

import matplotlib.pyplot as plt

from geneview.gwas import manhattanplot

ax = manhattanplot(data, xlabel="Chromosome", ylabel="-Log10(P-value)")  # 這就是Manhattan plot的函數(shù)

plt.show()

只需這樣的一句代碼就能創(chuàng)建一個(gè)漂亮的曼哈頓圖，有必要再次指出的是，geneview是以matplotlib為基礎(chǔ)開發(fā)出來的，所創(chuàng)建的圖形對(duì)象實(shí)際上仍屬于matplotlib，geneview內(nèi)部自定義了很多圖形風(fēng)格，同時(shí)封裝了大量只屬于基因組數(shù)據(jù)的圖表類型，但圖形的輸出格式以及界面顯示都仍和matplotlib一樣，因此在這里我們使用matplotlib.pyplot的show()函數(shù)(上例中：plt.show())將所繪制出來的曼哈頓圖顯示出來。如果要將圖形保存下來，則只需執(zhí)行`plt.savefig("man.png")`，這樣就會(huì)在該目錄下生成一個(gè)名為『man.png』png格式的曼哈頓圖，若是要存為pdf格式，則只需將所要保存的文件名后綴改成『.pdf』（plt.savefig("man.pdf")）就可以了。下面這些格式：emf, eps, pdf, png, jpg, ps, raw, rgba, svg, svgz等都是支持的，至于最新的還有多少種，還請(qǐng)參照matplotlib文檔中說明。

此外，geneview中的每個(gè)畫圖函數(shù)都有著足夠的靈活性，我們也可以根據(jù)自己的需要做一些調(diào)整，比如：

xtick = ['1', '2','3','4','5','6','7','8','9','10','11','12','13','14','16','18', '20','22']

manhattanplot(data,

              xlabel="Chromosome", # 設(shè)置x軸名字

              ylabel="-Log10(P-value)", # 設(shè)置y軸名字

              xtick_label_set = set(xtick), # 限定橫坐標(biāo)軸上的刻度顯示

              s=40, # 設(shè)置圖中散點(diǎn)的大小

              alpha=0.5, # 調(diào)整散點(diǎn)透明度

              color="#f28b1e,#9a0dea,#ea0dcc,#63b8ff", # 設(shè)置新的顏色組合

              ) 

實(shí)現(xiàn)新的顏色組合、限定x軸上的刻度顯示和散點(diǎn)大小的調(diào)節(jié)。甚至還可以將散點(diǎn)改為線：

manhattanplot(data，

              xlabel="Chromosome", # 設(shè)置x軸名字

              ylabel="-Log10(P-value)", # 設(shè)置y軸名字

              xtick_label_set = set(xtick), # 限定橫坐標(biāo)軸上的刻度顯示

              alpha=0.5, # 調(diào)整散點(diǎn)透明度

              color="#f28b1e,#9a0dea,#ea0dcc,#63b8ff", # 設(shè)置新的顏色組合

              kind="line"

              ) 

其它方面的調(diào)整請(qǐng)查看geneview文檔中的相關(guān)說明。

Q-Q圖

qq圖只需用到上例中的pvalue那一列：

import csv

import matplotlib.pyplot as plt

from geneview.gwas import qqplot

pvalue=[]

with open("GOYA.csv") as f:

    f_csv = csv.reader(f)

    headers = next(f_csv)

    pvalue = [float(row[3]) for row in f_csv]

ax = qqplot(pvalue, color="#00bb33", xlabel="Expected p-value(-log10)", ylabel="Observed p-value(-log10)") # Q-Q 圖

plt.show()

同樣，也可以根據(jù)自己的需要對(duì)改圖進(jìn)行相關(guān)的調(diào)整。

以上，便是如何使用Python來制作Manhattan圖和QQ圖的方法，geneview的集成函數(shù)簡化了這樣的一個(gè)過程。

另外，如果你也看過丹麥人的這個(gè)GOYA研究，就會(huì)發(fā)現(xiàn)實(shí)際以上的兩個(gè)圖和其文章中的基本是一致的，當(dāng)然我自己做了些數(shù)據(jù)清洗的操作，結(jié)果上仍然會(huì)有些許的不同。雖然此刻下結(jié)論還有點(diǎn)為時(shí)尚早，但總的來講，我應(yīng)該也可以通過這個(gè)數(shù)據(jù)集比較順利的將其結(jié)果重復(fù)出來了。

最后，附上利用geneview畫曼哈頓圖和QQ圖的代碼：

（1）曼哈頓圖：

（2）QQ圖：

關(guān)于怎樣使用Python繪制GWAS分析中的曼哈頓圖和QQ圖就分享到這里了，希望以上內(nèi)容可以對(duì)大家有一定的幫助，可以學(xué)到更多知識(shí)。如果覺得文章不錯(cuò)，可以把它分享出去讓更多的人看到。