抄読会でRNA-seqのデータ処理についての話を聞いて、その中で出てきたプロットが気になったのでpythonでやってみた。実際に可視化するときはRのパッケージを使ったほうが楽だと思う。
MAプロット
- 参考サイト
- MAプロット: MAプロット | マイクロアレイ解析やRNA-Seq解析で遺伝子発現量の偏りを視覚化,Rのサンプルコード
- マイクロアレイでの正規化: Normalization of microarray
- pandas: Pandas でデータフレームから特定の行・列を取得する – Python でデータサイエンス
- MAプロット: MAプロット | マイクロアレイ解析やRNA-Seq解析で遺伝子発現量の偏りを視覚化,Rのサンプルコード
参考サイトでは、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()
平均値が0のものを無理やりプロットしたものが以下(wikipediaより)
RAプロット
MAプロットと事実上同じなのだけれど、0になる値に$\epsilon$だけ足してやってからプロットするらしい。wikipediaより。
Bland-Altmanプロット
- 参考サイト
logを取らずに平均を横軸、差を縦軸にとってプロットしたもの。系統誤差を調べるときに用いる。