講義で教える、『お手軽マンハッタンプロット』
ゲノムワイドの関連検定データが出たとする。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の値を反映させてきちんとやることもできる(後述)が、「簡単マンハッタンプロット」では、マーカーの順序さえ守れればよいとすると次のように出来る。
１本の染色体の長さは、ゲノム全長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)
