抄読会でRNA-seqのデータ処理についての話を聞いて、その中で出てきたプロットが気になったのでpythonでやってみた。実際に可視化するときはRのパッケージを使ったほうが楽だと思う。

MAプロット

参考サイト
- MAプロット: MAプロット | マイクロアレイ解析やRNA-Seq解析で遺伝子発現量の偏りを視覚化,Rのサンプルコード
  - マイクロアレイでの正規化: Normalization of microarray
- pandas: Pandas でデータフレームから特定の行・列を取得する – Python でデータサイエンス

参考サイトでは、A、Bの2群について遺伝子ごとに値の平均を取った後、平均の「積のlogを表す横軸A」と「比のlogを表す縦軸M」を計算、プロットしている。2群は基本的に発現量は同じだと考えられるので、データを正規化してから利用する。正規化の方法としてはマイクロアレイデータだとLOWESS正規化が使われるらしい。今回は人ごとの発現量の違いだけ調整している。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# データはリンク先のものを使わせてもらった
data = pd.read_csv("http://bi.biopapyrus.net/data/count2g.txt", sep='\t')
"""
In [45]: data.head()
Out[45]: 
   GeneID   A1   A2   A3   B1  B2  B3
0  gene_1    3    2   40    4   2   0
1  gene_2  680  807  594  165  91  73
2  gene_3   35   38   31    7   9  11
3  gene_4  103  165  110   38  31  28
4  gene_5    2    9    1    3   1   0
"""

groupA = data.iloc[:,1:4].mean(axis=1)
groupB = data.iloc[:,4:7].mean(axis=1)

# 正規化
groupA = groupA * 1000000 / groupA.sum(axis=0)
groupB = groupB * 1000000 / groupB.sum(axis=0)

M = np.log2(groupA) - np.log2(groupB)
A = 1/2 * np.log2(groupA) + np.log2(groupB)

plt.scatter(A, M)
plt.show()