連続する要素の数を数える

ベクトル中の連続する要素の連続数は

rle

で計算できる。
　
いま、ある塩基配列

s1 <- c("ctgactgagagactttacgtaggagggcccatcgcgccagcggggtgctaaccactagcatgggtaacaccatcacatgctatgtgaaagccctagcggcctgcaaggctgcggggatagttgcgcccacgaaagccaggggactgagaggacgagcggaacctgagggaggccatgagtactctgcccctcctggtgat",
"ccccccagaccggaatatgcaaatgtgtccggggccgccgcagatactacctgaccagagacccaaccactccaccccgggctgcctgggaaacagttagacactcccctatcaattcatggctgggatatgggttcgcatggtcctaatgacacgtccaagacaccctggaccagaacctcaactttgagatgtatgga",
"tcagtatactccgtgaatcctttggaccttccagccataattgagaggttacacgggcttgacgccttttctatgcacacatactctcaccacgaactgacgcgggtggcttcagccctcagaaaacttggggcgccacccctcagggtgtggaagagtcgggctcgcgcagtcagggcgtccctcatctcccgtggagg",
"gaaagcggccgtttgcggccgatatctcttcaattgggcggtgaagaccaagctcaaactcactccattgccggaggcgcgcctactggacttatccagttggttcaccgtcggcgccggcgggggcgacatttttcacagcgtgtcgcgcgcccgaccccgctcattactcttcggcctactcctacttttcgtagggg",
"gtaggcctcttcctactccccgctcggtagagcggcacacactaggtacactccatagctaactgttcctttttttttttttttttttttttttttttttttttttttttttcttttttttttttttccctctttcttcccttctcatcttattctactttcttttggtggctccatcttagccctagtcacggctctaa")
s2 <- paste(s1, sep="", collapse="")

のAGCT各塩基が連続して出現する位置を確認して図示したいとする。

library(seqinr)
s2 <- toupper(s2) # 大文字にしただけ
s3 <- s2c(s2) # 一文字ずつバラバラにする
AGCT <- c("A", "G", "C", "T")

# rle関数で連続した塩基を計算
cont_max <- 12 # 連続しすぎると、2〜3個の連続する塩基が潰れてしまうので強制変換する上限
cols <- grey(seq(1, 0, length=10))
res0 <- matrix(0, nr=length(AGCT), nc=length(s3), dimnames=list(AGCT, seq(s3)))
for(i in seq(AGCT)){
	r0 <- rle(s3)
	r1 <- unlist(mapply(rep, r0$lengths, r0$length))
	res0[i, ] <- replace(r1, s3 != AGCT[i], 0)
}

res0[res0 >= cont_max] <- cont_max
layout(t(as.matrix(c(rep(1, 40), 2))))
par(mar = c(5, 2, 1, 0.5))
image(seq(s3), seq(AGCT), t(res0), col=cols, axes=FALSE, xlab="base position", ylab="")
axis(1, tick=FALSE, las=2, line=-0.5)
axis(2, at=seq(AGCT), label=AGCT, tick=FALSE, line=-1)

par(mar= c(5, 0, 1, 2))
image(y=0:max(res0), z=t(as.matrix(0:max(res0))), col=cols, axes=FALSE, ylab="")
axis(4, tick=FALSE, las=2, line=-0.8)

Tは配列の5'側にはほとんどなくて、3'側に密集していることがわかる。