講義で教える、『お手軽マンハッタンプロット』

ゲノムワイドの関連検定データが出たとする。DNA上の物理位置を横軸に、関連検定統計量を縦軸にプロットするのはルーチン。
ルーチンはだれでもできることが望ましい(できれば自分でやらずにだれかにやって欲しい…)
次のようなファイルがテキストファイルがあるとする。染色体番号とその塩基番号と統計量(この場合はp)

ch	nt	p
20	4752	0.147597304
17	6753	0.085212024
20	3723	0.993693715
19	1183	0.608450081
1	3327	0.058905169
19	6076	0.670302559
3	8968	0.521401245
2	5699	0.31427928
21	6782	0.189641514
16	1244	0.875761717
7	9385	0.487822156
11	5611	0.802256673
19	935	0.612314601
22	5513	0.895206958
19	1852	0.360684214
4	9926	0.473783418
20	8646	0.898933488
18	304	0.486757311
...

chのカラムとntのカラムの値を使って横軸の値を決め、それを使ってプロットするわけである。

まず、このファイル"data.txt" (タブ区切り)をヘッダーつきで読み込む。

data<-read.table(file="data.txt",header=T)

data$chに染色体番号、data$ntに染色体上の塩基番号、data$pにp値が入っている。
物理地図の位置をntの値を反映させてきちんとやることもできる(後述)が、「簡単マンハッタンプロット」では、マーカーの順序さえ守れればよいとすると次のように出来る。
1本の染色体の長さは、ゲノム全長3*10^9より短いことは確かめなくてもわかるから、染色体番号にゲノム全長をかけた値の下駄をはかせて、位置の順序を取る

order<-order(data$ch*3*10^9+data$nt)

この順序でプロットする

plot(data$p[order])

もっと、縮めて

genomelength<-3*10^9
plot(data$p[order(data$ch*genomelength+data$nt)])

物理位置を反映させるやり方。
今、染色体のnt長さの最大値をRに教えてやる必要がある。X染色体を23番扱いして、23個の数字で以下のように。

chlen<-c(279000000,251000000,221000000,197000000,198000000,176000000,163000000,148000000,140000000,143000000,148000000,142000000,118000000,107000000,100000000,104000000,88000000,86000000,72000000,66000000,45000000,48000000,163000000)

染色体番号が後になるに従って、累計していかないといけないのであるが、それをループを使うよりは、行列を使う方がRらしいので、そうする。

mat<-matrix(1,ncol=23,nrow=23)
mat[upper.tri(mat)]<-0
accumchlen<-mat%*%chlen

accumchlen[data$ch]が染色体番号に応じて加算されるべき値なので

horizontalvalue<-accumchlen[data$ch]+data$nt

で横軸の値が決まる

plot(horizontalvalue,data$p)