Pythonで作るVolcano plot

# Volcano plotとは
Volcano plotは、RNA-seqやマイクロアレイで二群の遺伝子発現量を比較する際に、遺伝子の発現比と統計的有意性（p値）でプロットした図です。

x軸を発現比、y軸を統計的有意性としたときのプロットが一般的な描き方です。

何千、何万という遺伝子の発現を視覚化して、重要な変化をしている遺伝子を特定することができます。

この図だと赤が統計的有意に発現量が増加した遺伝子で、青が統計的優位に発現量が減少した遺伝子を示しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/80f36454-5f0a-d806-52a7-fdde49fd4bd0.png)

エクセルでも作れますが、Pythonで**bioinfokit**ライブラリを使えば、めっちゃ簡単に、しかもかっこいい図が作れちゃいます。

今回はこのVolcano plotの作り方を記していきます！

# ライブラリのインストール
**pip**でも**conda**でもインストールできますが、今回は**conda**インストールをやっていきます。

下のコマンドを実行します。

```console
conda install bioinfokit
```

これで**bioinfokit**ライブラリはインストールできましたが、**bioinfokit**ライブラリは他に色々なライブラリに依存しているので、入ってなければ必要なほかのライブラリ（**adjusttext**とか**textwrap3**など）も入れていきます。

このあと``from bioinfokit import analys, visuz``を実行したときに、足りないものがあるとerrorになるので、errorとして出てきたものをインストールしていけば大丈夫です。

下のリンクの「Depends」を見れば必要なライブラリ一覧がわかります。

[bioinfokit](https://bioconda.github.io/recipes/bioinfokit/README.html)

# データの準備

必要なデータセットは、1) 遺伝子名、2) 比較したい二群（A、B）の遺伝子ごとのt検定の結果（p値）、3) 二群の発現比（A/B）です。つまり3カラム分です。

今回のデモデータは以下のURLのものを使います。

[データセット](https://github.com/vappiah/bioinfoscripts/blob/main/testvolcano.csv)

# Volcano plotを描いてみる
早速描いていきます！

まずはライブラリのインポートです。使うライブラリは**pandas**と**bioinfokit**の二種類です。

```python
import pandas as pd
from bioinfokit import analys, visuz
```

ここで、さっき言ったライブラリ（**adjusttext**とか**textwrap3**など）が入っていないというエラーがあれば、condaインストールで入れていきます。

つぎにデータを読み込んで、中身を確認します。

```python
df = pd.read_csv('testvolcano.csv')
print(df.head())
```

```python:出力
          GeneNames    log2FC       p-value
0  LOC_Os09g01000.1 -1.886539  1.250000e-55
1  LOC_Os12g42876.1  3.231611  1.050000e-55
2  LOC_Os12g42884.2  3.179004  2.590000e-54
3  LOC_Os03g16920.1  5.290677  4.690000e-54
4  LOC_Os05g47540.4  4.096862  2.190000e-54
```

遺伝子名（GeneNames）、遺伝子の発現比（log2FC）、p値（p-value）の必要な3カラムを確認しました。

これでひとまず、必要なものは揃いました！

めっちゃ簡単です！

さっそく描いていきましょう！

引数は1) データフレーム名（df）、2) 発現比（log2FC）、3) p値（p-value）を与えるだけです。

```python
visuz.GeneExpression.volcano(df=df, lfc='log2FC', pv='p-value')
```

そうすると、、、

![volcano.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/016fa75e-d48b-3d39-9698-008d69dc1381.png)

こんな簡単に描けました！しかも、10行以内のコードで！

図は作業ディレクトリにPNG形式で保存されています。

緑色が統計的有意（p < 0.05）に発現が2倍以上増加した遺伝子群、赤色が統計的に有意に発現が2倍以上減少した遺伝子群です。

# 整えます
せっかくpythonで作ったのにこのままだと味気ないので、色々いじってみます。

まずは色を変えてみます。発現が上がったらなんとなく赤色で、下がったらなんとなく青色な先入観があるので変えてみます。

```python:プロットの色を変える
visuz.GeneExpression.volcano(df=df, lfc='log2FC', pv='p-value',
                             color=('red', 'grey', 'blue')) #👈
```

![volcano.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/503beb81-2aa6-714a-8211-2c4ed42f9a0a.png)

しきい値がどこらへんなのか、ラインがあるとわかりやすいですよね。しきい値ライン、入れてみます。

```python:しきい値を入れる
# ボルケーノプロット3
visuz.GeneExpression.volcano(df=df, lfc='log2FC', pv='p-value',
                             color=('red', 'grey', 'blue'),
                             sign_line=True) #👈
```

![volcano.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/70ab230b-3f5d-a7c1-8be0-09e3e4102537.png)


しきい値自体を変えることもできます。

デフォルトだと、低2の対数をとった発現比（log2FC）が1以上（つまり発現比に戻すと2倍以上）、p値が0.05未満になっています。

発現減少側（青色）のしきい値は変更せずに、発現増加側（赤色）のしきい値をlog2FCが2以上、p値が0.01未満としてみましょう。

```python:しきい値を変える
visuz.GeneExpression.volcano(df=df, lfc='log2FC', pv='p-value',
                             color=('red', 'grey', 'blue'),
                             sign_line=True,
                             lfc_thr=(2, 1), #👈 発現比
                             pv_thr=(0.01, 0.05) #👈 p値
)
```

![volcano.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/02ef5548-0cad-54dd-bbb3-b1212ee15bc6.png)

しきい値が二段階になっているので、ラインが機能していません。こういう場合は、ラインは消したほうが良さそうですね。

プロットに遺伝子名を付与していきます。引数``geneid``に遺伝子名のカラムを指定して、入れたい遺伝子名を引数``genenames``で指定していきます。

```python:プロットに遺伝子名を付与する
visuz.GeneExpression.volcano(df=df, lfc='log2FC', pv='p-value', geneid='GeneNames', #👈
                             color=('red', 'grey', 'blue'),
                             sign_line=True,
                             genenames=('LOC_Os12g42876.1', 'LOC_Os09g01000.1', 'LOC_Os09g27030.2') #👈
)
```

![volcano.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/471633/1c4ce525-f2ac-fabe-ec12-5f9abf3042ee.png)

簡単にプロットと遺伝子名をリンクすることができました！

めっちゃ簡単です！

ちなみに、SpyderなどのIDEを使っていれば、引数に``show=True``と入れれば、プロットのところに図が出力されます。

他にも色々ビジュアルを変えられるので試してみてください！