読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

MAプロット、RAプロット、Bland-Altmanプロット

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

MAプロット

参考サイトでは、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()

f:id:kuyata:20160726010430p:plain

平均値が0のものを無理やりプロットしたものが以下(wikipediaより)

https://upload.wikimedia.org/wikipedia/commons/thumb/c/c0/MaPlot-edgeR.smear-wikipedia.png/450px-MaPlot-edgeR.smear-wikipedia.png

RAプロット

MAプロットと事実上同じなのだけれど、0になる値に$\epsilon$だけ足してやってからプロットするらしい。wikipediaより。

https://upload.wikimedia.org/wikipedia/commons/thumb/b/b3/RaPlot-orange.uniques-wikipedia.png/450px-RaPlot-orange.uniques-wikipedia.png

Bland-Altmanプロット

logを取らずに平均を横軸、差を縦軸にとってプロットしたもの。系統誤差を調べるときに用いる。