めちゃくちゃかっこいいvolcano plotを作りたい！

# volcano plotとは
Volcano plotは、一般的にはx軸を発現比、y軸を統計的有意性とした散布図の一種で、主にRNA-seqやDNAマイクロアレイなどで得られた二群の遺伝子発現量を比較する際に用いられます。

数千、数万という遺伝子発現の変化を可視化することにより、重要な要素を特定する助けとなります。
例として、この図だと赤が統計的有意に発現量が増加した遺伝子で、青が統計的優位に発現量が減少した遺伝子を示しています。
![volcanoplot_example](https://camo.qiitausercontent.com/176ca187fe2e3d2b9a58891ff4891081f27ea7a6/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e61702d6e6f727468656173742d312e616d617a6f6e6177732e636f6d2f302f3437313633332f38306633363435342d356630612d643830362d353261372d6664646534396664346264302e706e67)

今回はこのvolcano plotをもっとカッコ良くしていこうと思います。([元データ](https://training.galaxyproject.org/archive/2022-04-01/topics/transcriptomics/tutorials/rna-seq-viz-with-volcanoplot-r/tutorial.html#import-data))

完成図
![Qiita_test.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2650928/bcce6e75-881a-f348-f7d5-62a399da9291.png)


# 実践
必要なライブラリのインポート
```python
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は必須じゃないです
```
次にデータを読み込んで、必要なデータのみを抽出します。

```python
df = pd.read_csv('volcanoplot_test.csv',usecols = [1, 3, 6]).dropna()
#データの読み込みおよび欠損値NaNを含む行の削除
df = df.rename(columns={'P.Value' : 'qvalue'})
#列名の変更
df.head()
```
```python:出力
   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)も追加します。
```python
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()
```
```python:出力
	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の色分けのために新たなカラムを作成し、遺伝子発現量の差が統計的に有意であるものとそうでないものを分ける。
```python
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()
```
```python:出力
    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を描画してみる。
```python
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()
```
出力
![Qiita_test.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2650928/bcce6e75-881a-f348-f7d5-62a399da9291.png)

出来上がり！

どうでしょう...
少しごちゃごちゃしているような気もしますが、dotの形とサイズを新たな変数として設けたことで、一枚の図から個々の遺伝子発現量の規模や産物の種類までがわかるようになりました！