vulgarから描く¶
vulgarはexonerateが出力するマッチングの1行表示です。以下のような感じです。
In [1]:
scan("AT5G22650.vulgar", sep = "\n", what = 'c')
vulgarの項目はスペースで区切られています。 先頭の4個と続く4個が比較した二つの配列の情報を示していて、9番目はマッチングのスコアを示しています。 先頭の4項目、続く4項目は「配列の名前」、「どこから」、「どこまで」、「マッチの方向」から成ります。
10番目以降がマッチした領域の情報で、3個一組になっています。 各組の1番目がMはマッチ、Iはイントロン、Gはギャップなどになっていて、2番目と3番目は二つの塩基配列のマッチした長さになっています。
例えば100bpと200bpのエキソンが150bpのイントロンで分断されていると、M 100 100 5 0 2 I 0 146 3 0 2 M 200 200 となります。M 100 100はcDNAとゲノムが100塩基ずつマッチしていて、5 0 2はイントロンの5'側がcDNAの0塩基とゲノムの2塩基(GT)が一致していることを示しています。
構造を図示するのに必要な情報は10番目以降で、3個一組を配列(3行)に並び直します。
In [2]:
vulgar <- scan("AT5G22650.vulgar", sep = " ", what = 'c')
vulgar <- matrix(tail(vulgar, - 9), nrow = 3)
vulgar
数値の型を修正し、累積します。Rで図(四角)を描くには起点・長さではなく起点(f)・終点(t)のほうがやりやすいためです。
In [3]:
t <- cumsum(as.numeric(vulgar[3,]))
f <- t - as.numeric(vulgar[3,])
head(cbind(vulgar[1,],f, t))
エキソンの情報に絞って、図を描きます。
In [4]:
exons <- data.frame(f, t)[vulgar[1,] == 'M',]
options(repr.plot.width=10, repr.plot.height=2.5)
plot(c(0, max(t)), c(0.5, 0.5), type = 'l', lwd = 10,
xlim = c(0, max(exons$t)), ylim = c(0, 1),
axes = F, xlab = '', ylab = ''
)
rect(exons$f, 0, exons$t, 1, col = 'gray')
エキソンをつないだようにするには以下のようにしてみます。
In [5]:
plot(c(0, max(t)), c(0.5, 0.5), type = 'n',
xlim = c(0, max(exons[,'t'])), ylim = c(-0.5, 1),
axes = F, xlab = '', ylab = ''
)
rect(exons$f, 0, exons$t, 1, col = 'gray')
introns <- data.frame(f = head(exons$t, -1), t = tail(exons$f, -1))
apply(introns, 1, function (intron){
points(c(intron[1], mean(intron), intron[2]), c(0, -0.5, 0), type = 'l')
})
