volcano plotとは
Volcano plotは、一般的にはx軸を発現比、y軸を統計的有意性とした散布図の一種で、主にRNA-seqやDNAマイクロアレイなどで得られた二群の遺伝子発現量を比較する際に用いられます。
数千、数万という遺伝子発現の変化を可視化することにより、重要な要素を特定する助けとなります。
例として、この図だと赤が統計的有意に発現量が増加した遺伝子で、青が統計的優位に発現量が減少した遺伝子を示しています。
今回はこのvolcano plotをもっとカッコ良くしていこうと思います。(元データ)
完成図
実践
必要なライブラリのインポート
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from adjustText import adjust_text
import random
#randomは必須じゃないです
次にデータを読み込んで、必要なデータのみを抽出します。
df = pd.read_csv('volcanoplot_test.csv',usecols = [1, 3, 6]).dropna()
#データの読み込みおよび欠損値NaNを含む行の削除
df = df.rename(columns={'P.Value' : 'qvalue'})
#列名の変更
df.head()
出力
SYMBOL logFC qvalue
0 COX1 -0.234558 0.996856
1 COX3 -0.132968 0.999574
2 CYTB -0.203188 0.997867
3 COX2 -0.273837 0.993166
4 ATP6 -0.156514 0.999473
仮遺伝子名（SYMBOL）、遺伝子の発現比（logFC）、q値（q-value）の必要な3カラムを確認。
次に各遺伝子の正規化された平均発現量(basemean)、遺伝子産物（biotype）をランダムに新たな要素として作成し(元データにはなかったので、あると仮定して作成)、y軸用の値(-log10q-value)も追加します。
df['basemean'] = [random.uniform(1,100) for i in range(len(df))]
#各遺伝子にランダムなbasemeanを追加する
df['nlog10'] = -np.log10(df.qvalue)
#y軸用の値
l = ['protein-coding', 'miRNA', 'siRNA','piRNA']
df['biotype'] = [random.choice(l) for i in range(len(df))]
#各遺伝子にランダムなbiotypeを追加する
df.head()
出力
SYMBOL logFC qvalue basemean nlog10 biotype
0 COX1 -0.234558 0.996856 19.200694 0.001368 piRNA
1 COX3 -0.132968 0.999574 21.929864 0.000185 siRNA
2 CYTB -0.203188 0.997867 6.621319 0.000927 protein-coding
3 COX2 -0.273837 0.993166 99.418353 0.002978 miRNA
4 ATP6 -0.156514 0.999473 99.971236 0.000229 piRNA
ここまでで元データに必要なデータの追加が終了。
ここからはvolcano plotを書くためのデータ整形・加工になります。
volcano plotの色分けのために新たなカラムを作成し、遺伝子発現量の差が統計的に有意であるものとそうでないものを分ける。
def dot_color(a):
logFC, nlog10 = a
if abs(logFC) <= 1 or nlog10 <= 2:
return 'Not significant' #有意でないもの
if logFC > 1 and nlog10 >2:
return 'Up' #有意に発現量が上がっているもの
if logFC < -1 and nlog10 >2:
return 'Down' #有意に発現量が減少しているもの
df['color'] = df[['logFC', 'nlog10']].apply(dot_color, axis = 1)
#新たなカラムcolorを作成し、各要素を追加する。
df.head()
出力
SYMBOL logFC qvalue basemean nlog10 biotype color
0 COX1 -0.234558 0.996856 19.200694 0.001368 piRNA Not significant
1 COX3 -0.132968 0.999574 21.929864 0.000185 siRNA Not significant
2 CYTB -0.203188 0.997867 6.621319 0.000927 protein-coding Not significant
3 COX2 -0.273837 0.993166 99.418353 0.002978 miRNA Not significant
4 ATP6 -0.156514 0.999473 99.971236 0.000229 piRNA Not significant
これでデータフレームの完成。
実際に使用してvolcano plotを描画してみる。
plt.figure(figsize = (10,10)) #サイズの指定
#各引数に作成したdataframeのカラム名を追加していきます。
#色は遺伝子の発現量で、dotの形はbiotypeで、dotのサイズはbasemeanで指定
ax = sns.scatterplot(data = df, x = 'logFC', y = 'nlog10',
hue = 'color', hue_order = ['Not significant', 'Up', 'Down'],
palette = ['dimgrey', 'firebrick', 'darkslateblue'],
style = 'biotype', style_order = ['protein-coding', 'miRNA', 'siRNA','piRNA'],
markers = ['o', 's', '^', 'X'],
size = 'basemean', sizes = (10, 200))
#統計的に有意な範囲を点線で区切る
ax.axhline(2, zorder = 0, c = 'k', lw = 2, ls = '--')
ax.axvline(1, zorder = 0, c = 'k', lw = 2, ls = '--')
ax.axvline(-1, zorder = 0, c = 'k', lw = 2, ls = '--')
#一定以上の統計的優位性を持つもののにannotationとして遺伝子名を表示させる
texts = []
for i in range(len(df)):
if df.iloc[i].nlog10 > 2 and abs(df.iloc[i].logFC) > 6:
texts.append(plt.text(x = df.iloc[i].logFC, y = df.iloc[i].nlog10, s = df.iloc[i].SYMBOL,
fontsize = 12, weight = 'bold'))
adjust_text(texts, arrowprops = dict(arrowstyle = '-', color = 'k'))
plt.legend(loc = 1, bbox_to_anchor = (1.2,1), frameon = False, prop = {'weight':'bold'})
for axis in ['bottom', 'left']:
ax.spines[axis].set_linewidth(2)
#上と右の囲いを消す
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.tick_params(width = 2)
plt.xticks(size = 12, weight = 'bold')
plt.yticks(size = 12, weight = 'bold')
plt.xlabel("$log_{2}$ fold change", size = 15)
plt.ylabel("-$log_{10}$ q-value", size = 15)
#グラフの保存
plt.savefig('Qiita_test.png', dpi = 300, bbox_inches = 'tight')
plt.show()
出力
出来上がり！
どうでしょう...
少しごちゃごちゃしているような気もしますが、dotの形とサイズを新たな変数として設けたことで、一枚の図から個々の遺伝子発現量の規模や産物の種類までがわかるようになりました！