バイオインフォマティクス、勉強中

エンコーディングの問題

2017-07-21T00:00:00.000+09:00

たまたま入手したファイルが開けず苦労した。
mi で調べると、エンコーディングが Windows-1252 だった。

これを R で解決する方法。
readLines() で encoding を指定して読み込み、write() で書き出す。

gff <- readLines("XXX.gff3", encoding = "windows-1252")
write(gff, file = "XXX_UTF8.gff3")

これでとりあえず問題回避できる。

homebrew: tbl2asn installation failed

2017-07-06T00:00:00.000+09:00

homebrew で tbl2asn をインストールできないときの対処法。
When installation of tbl2asn through homobrew failed.

Error: SHA256 mismatch
Expected: 429d63ee3c36d1f2f6322c62c6089d5ee8a8b089e5cc9373e298e017bcbbb9ec
Actual: bfaff2308fd58b94730639b9c6b520ca1884ede9979389cfc9bc4902ae75702c
Archive: /Users/XXX/Library/Caches/Homebrew/tbl2asn-25.3.gz
To retry an incomplete download, remove the file above.

下のファイルをテキストで開いて、変更を加える。
Open the .rb file as a text and change the sha256 as bellow.

File location:
/usr/local/Homebrew/Library/Taps/homebrew/homebrew-science/tbl2asn.rb

### original
# if OS.mac?
# url "ftp://ftp.ncbi.nih.gov/toolbox/ncbi_tools/converters/by_program/tbl2asn/mac.tbl2asn.gz"
# sha256 "429d63ee3c36d1f2f6322c62c6089d5ee8a8b089e5cc9373e298e017bcbbb9ec"

### fixed
# if OS.mac?
# url "ftp://ftp.ncbi.nih.gov/toolbox/ncbi_tools/converters/by_program/tbl2asn/mac.tbl2asn.gz"
# sha256 "bfaff2308fd58b94730639b9c6b520ca1884ede9979389cfc9bc4902ae75702c"

R: pdflatex is not available

2016-08-03T00:00:00.000+09:00

Rで「pdflatex is not available」というエラーに遭遇したときの対処方法。
こちらのサイトを参考にした。

まず、ターミナルで「pdflatex」の場所を確かめる。whichを使う。

$ which pdflatex
/usr/texbin/pdflatex

ちなみに「R」の場所はこちら。

$ which r
/usr/bin/r

「pdflatex」が入っていることを確認した。mac側ではなくR側の問題ということになる。

Rがソフトウェアを認識するかどうかは、Sys.which()で知ることができる。
たとえば「R」を探すなら、

> Sys.which("r")
r
"/usr/bin/r"

という具合。ターミナルでみたときと同じ。

上のエラーが出ている状態で「pdflatex」を探してみると、

> Sys.which("pdflatex")
pdflatex
""

どうやら認識していないことがわかる。

これはPATHの問題。RのPATHを調べるには、Sys.getenv()を使う。

> Sys.getenv("PATH")
[1] "/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin"

いくつかPATHがあるが、「/usr/texbin」がない。
これがエラーの元凶。

PATHを通す作業を行う。そのためには、Sys.setenv()を使う。
既存のPATHと「/usr/texbin」をpaste()で繋げるところがRっぽい。

Sys.setenv(PATH = paste(Sys.getenv("PATH"), "/usr/texbin", sep=":"))

これで、PATHが追加された。

> Sys.getenv("PATH")
[1] "/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/usr/texbin"

Rが「pdflatex」の場所を知ることができるようになった。

> Sys.which("pdflatex")
pdflatex
"/usr/texbin/pdflatex"

これで、「pdflatex is not available」というエラーが出なくなる。

R: オブジェクトのサイズを考える

2016-07-22T00:00:00.000+09:00

関数object.size()を使えば、オブジェクトのサイズを知ることができる。大きなデータを扱う場合、サイズは意識しておいた方が良いと思う。

とりあえず、character vectorをつくってサイズを眺めてみる。

# > x <- ""
# > nchar(x)
# [1] 0
# > object.size(x)
# 96 bytes

# > object.size("1")
# 96 bytes

# > object.size("1234567")
# 96 bytes

# > object.size("12345678")
# 104 bytes

# > object.size("123456789012345")
# 104 bytes

# > object.size("1234567890123456")
# 120 bytes

# > object.size("1234567890123456789012345678901")
# 120 bytes

# > object.size("12345678901234567890123456789012")
# 136 bytes

文字が長くなると表示しづらいので、for()で伸ばしてnchar()で文字数を確認し、object.size()で調査する。

# > x <- ""
# > for(i in 1:47) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 47
# > object.size(x)
# 136 bytes

# > x <- ""
# > for(i in 1:48) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 48
# > object.size(x)
# 152 bytes

# > x <- ""
# > for(i in 1:63) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 63
# > object.size(x)
# 152 bytes

# > x <- ""
# > for(i in 1:64) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 64
# > object.size(x)
# 216 bytes

文字の数とサイズは、どうやら以下の関係になっている。
0~7: 96 bytes
8~15: 104 bytes (+8)
16~31: 120 bytes (+16)
32~47: 136 bytes (+16)
48~63: 152 bytes (+16)
64~?: 216 bytes (+64)

一見すると秩序がありそう。
ということで、以下のコマンドでプロットしてみる。

main <- "Relationship between nchar and object.size"
sub <- c("64(blue), 128(green), 136(red), 144(orange)")
y <- ""
plot(x = nchar(y), y = object.size(y),
sub = sub, main = main,
xlim = c(0,550), ylim = c(0,700), pch = 20)
for(i in 1:550){
y <- paste(y, "y", sep = "")
points(x = nchar(y), y = object.size(y), pch = 20)
}
abline(v = c(64, 128, 136, 144),
col = c("blue", "green", "red", "orange"))

これはもう、本当に変なプロットになった。
どうやら、127文字を越えると取扱い方が変わるらしい。

# > x <- ""
# > for(i in 1:127) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 127
# > object.size(x)
# 216 bytes

# > x <- ""
# > for(i in 1:128) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 128
# > object.size(x)
# 224 bytes

# > x <- ""
# > for(i in 1:135) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 135
# > object.size(x)
# 224 bytes

# > x <- ""
# > for(i in 1:136) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 136
# > object.size(x)
# 232 bytes

# > x <- ""
# > for(i in 1:143) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 143
# > object.size(x)
# 232 bytes

# > x <- ""
# > for(i in 1:144) x <- paste(x, "a", sep = "")
# > nchar(x)
# [1] 144
# > object.size(x)
# 240 bytes

まとめると、文字の数とサイズは以下の関係になっている。
0~7: 96 bytes
8~15: 104 bytes (+8)
16~31: 120 bytes (+16)
32~47: 136 bytes (+16)
48~63: 152 bytes (+16)
64~127: 216 bytes (+64)
128~135: 224 bytes (+8)
136~143: 232 bytes (+8)
144~152: 240 bytes (+8)
153~xxx: 8文字毎に8バイト増える。

基本的に１バイトに１文字が格納される計算になるのだが、そのオブジェクトの占有メモリは１文字毎に増えるのではない。あらかじめ決められた規則に則って、時には16文字分、あるいは64文字分のメモリが割り当てられる。しかも、どうやら128文字目以降は8文字単位でメモリの割り当てが増やされるらしい。

もうひとつの奇妙な現象を見つけた。
以下のコマンドを実行してもらうとわかるのだが、ひとつのベクターが２つの文字列をもち、しかもそれらの文字列が同一であれば、オブジェクトサイズは期待値よりも小さい。もしかしてfactorとして管理してるのかなぁ。

x <- ""
for(i in 1:10) x <- paste(x, "a", sep = "")
y <- c(x, x)
z <- c(x, paste(x, "b", sep = ""))
w <- c(paste(x, "b", sep = ""), paste(x, "b", sep = ""))
x
y
z
w
nchar(x)
nchar(y)
nchar(z)
nchar(w)
object.size(x)
object.size(y)
object.size(z)
object.size(w)

R: latticeのlevelplot()で枠を消す

2016-07-20T00:00:00.000+09:00

表に入っているデータを俯瞰したいとき、レベルプロットが役に立つ。レベルプロットはlatticeパッケージのlevelplot()で実現できる。

しかし通常のやり方だと関数内で縦と横がスケーリングされてしまうので、ピクセル単位で期待通りのサイズの絵を出力するのが難しい。

マージンをゼロにしてプロット領域のみにし、プロットの外枠を消し、カラーキーを消すと、１ピクセルに１つの値を配置することができる。

R: lattice levelplot(), Remove the border line.

When you would like to see the information in a table or matrix at a glance, draw a level plot. It can be done with the function levelplot() in the lattice package.

It seems that the function automatically scales the plot, making it difficult to get a plot with defined scales.

By removing the margins, border lines and colorkey, you can get a plot in which a square of 1 x 1 pixel is filled with a value-related color.

# テスト用のマトリックスを作る。乱数の絶対値。
# Make a matrix with randomly produced values for a test.
x <- rnorm(20, mean = 0.5, sd = 0.3)
x <- abs(x)
m <- matrix(x, nrow = 4)

# パッケージをロードする。
# Load the package.
library(lattice)

# マージンゼロのテーマを定義する。参考サイトはこちら。
# Define the theme.
theme.x <- list(
layout.heights = list( # マージンをゼロに
top.padding = 0,
main.key.padding = 0,
key.axis.padding = 0,
axis.xlab.padding = 0,
xlab.key.padding = 0,
key.sub.padding = 0,
bottom.padding = 0),
layout.widths = list(
left.padding = 0,
key.ylab.padding = 0,
ylab.axis.padding = 0,
axis.key.padding = 0,
right.padding = 0
),
axis.line = list(lwd = 0) # 外枠を消す
)

# 保存先に移動する。
# Change the working directory.
setwd("test")

# レベルプロットをPNGとして保存する。
# Draw a level plot and save it as a PNG.
width <- nrow(m)
height <- ncol(m)
filename <- "m.png"
png(filename = filename, width = width, height = height)
levelplot(m, par.settings = theme.x, xlab = NULL, ylab = NULL,
scales = list(draw = FALSE), colorkey = FALSE)
dev.off()

こうして得たPNGファイルは、X軸とY軸の表示もなく、色と値の関係も分からないが、１ピクセルに１つの値が入っているので取扱いやすい。
# The resulting PNG file does not have the XY-axises or colorkey. However, as it is completely in a pixel-based format, it may be good for some purposes.

と、思いきや、Keynote ver. 6 以降に張り付けようとすると、PNGの画質が劣化する問題が生じる。この問題を回避するには、png()ではなくpdf()を使用するのが良さそう。pdf()は細い外枠が加えられるが、おそらく支障のないレベル。

# レベルプロットをPDFとして保存する。
# pdf()のサイズはピクセルではなくインチなので注意。
width <- nrow(m) / 72
height <- ncol(m) / 72
filename <- "m.png"
pdf(filename = filename, width = width, height = height)
levelplot(m, par.settings = theme.x, xlab = NULL, ylab = NULL,
scales = list(draw = FALSE), colorkey = FALSE, useRaster = FALSE)
dev.off()

pdf()を使うとき、levelplot()にuseRaster=TRUEを与えると、keynoteに載せたときに絵が変になる。

R: strsplit() で、ピリオドを区切りとして文字列を分ける

2016-07-12T00:00:00.001+09:00

strsplit()は、指定した文字を区切りとして文字列を分割する関数。
splitというArgumentで区切り文字を指定する。

ピリオド「.」はワイルドカードとして認識されてしまう、というのが今回の問題。手元にピリオドを含む文字列があり、ピリオドを区切りとして２つの文字列に分割したいとき、エスケープする必要がある。

# たとえば「ABC.01」という文字列を考える。
x <- "ABC.01"

# ピリオドをそのまま与えると、ワイルドカードになる。
strsplit(x, split = ".")[[1]]

# 「\\」でエスケープすると、ピリオドが区切りとみなされる。
strsplit(x, split = "\\.")[[1]]

R: Map() でマップする

2016-07-12T00:00:00.000+09:00

Map the values with Map().

Map()はapply系の関数。mapply()のラッパー。

Map()は沢山のArgumentsを受け付けてくれるのでとても便利。
マルチコア版のmcMap()もある。

Map() is a wrapper function to mapply. Map() is useful as it takes lots of arguments. You can also use mcMap(), the multicore version of it.

#######################################################
たとえば、以下のようなデータフレームがあるとする。
Assume that there is a data frame like this.

x <- data.frame(strA = c("A", "B", "C", "D", "E", "F"),
strB = c("a", "b", "c", "d", "e", "f"),
strC = c("1", "2", "3", "4", "5", "6"))

> x
strA strB strC
1 A a 1
2 B b 2
3 C c 3
4 D d 4
5 E e 5
6 F f 6

各row毎に、strAとstrBとstrCを繋げた文字列を得たいとする。
ついでに任意の文字（ここでは「-」）を加える。
And, you would like to concatenate character strings in the columns strA, strB and strC. Here, you may want the characters to be linked with '-'.

当然ながら、for()ループで実現できる。
You can do it with the for().

x$strABC <- ""
for(i in 1:nrow(x)){
strABC <- paste(x$strA[i], x$strB[i], x$strC[i], sep = "-")
x$strABC[i] <- strABC
}

for()ループの結果。
The result of the for loop.

> x
strA strB strC strABC
1 A a 1 A-a-1
2 B b 2 B-b-2
3 C c 3 C-c-3
4 D d 4 D-d-4
5 E e 5 E-e-5
6 F f 6 F-f-6

しかし、rowの数が増えるほど計算が遅くなってしまう。Map()だと速い。
However, the more row numbers you have, the more time it takes. Map() does this more rapidly.

Map()を使うには、まず関数をつくる必要がある。
In order to use Map(), you need to make a function to be used by the function Map().

get.strABC <- function(strA, strB, strC){
return(paste(strA, strB, strC, sep = "-"))
}

この関数を使ってみるとこんな結果が帰ってくる。
This function acts like this.

> get.strABC("A", "a", "1")
[1] "A-a-1"

上のfor()ループの中でこの関数を使うと同じ結果を得ることができる。
Of course, you can use this function in a for loop.

x$strABC <- ""
for(i in 1:nrow(x)){
x$strABC[i] <- get.strABC(x$strA[i], x$strB[i], x$strC[i])
}

Map()でget.strABC()を実行する。Map()はリストで結果を返すので、unlist()で結果をシンプルにする必要がある。
Let's use Map(). You need to unlist the output of Map() to map the values in a preexisting container.

x$strABC <- unlist(Map(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC))

> x
strA strB strC strABC
1 A a 1 A-a-1
2 B b 2 B-b-2
3 C c 3 C-c-3
4 D d 4 D-d-4
5 E e 5 E-e-5
6 F f 6 F-f-6

Map()のマルチコア版であるmcMap()も使える。この場合、parallelパッケージをロードする必要がある。
You can also use mcMap(), which is the multicore version of Map().
You need load the parallel package to do it.

library(parallel)
x$strABC <- unlist(mcMap(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC))

mcMap()のArgumentであるmc.coresで、使うコアの数を制御できる。
The number of cores can be regulated with the argument mc.cores.

#######################################################
rowの数を増やしてfor()とMap()を比べる。

x <- data.frame(strA = c("A", "B", "C", "D", "E", "F"),
strB = c("a", "b", "c", "d", "e", "f"),
strC = c("1", "2", "3", "4", "5", "6"))

for(i in 1:12) x <- rbind(x, x)

> nrow(x)
[1] 24576

for()ループでやると、4.286秒かかる。

x$strABC <- ""
system.time(
for(i in 1:nrow(x)){
x$strABC[i] <- get.strABC(x$strA[i], x$strB[i], x$strC[i])
}
)

ユーザシステム経過
4.169 0.137 4.286

Map()でやると、0.940秒。for()より速い。

system.time(
x$strABC <- unlist(Map(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC))
)

ユーザシステム経過
0.892 0.009 0.940

mcMap()でやると、0.711秒。Map()との差は小さい。

system.time(
x$strABC <- unlist(mcMap(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC, mc.cores = 4))
)

ユーザシステム経過
1.697 0.155 0.711

#######################################################
rowの数を増やしてMap()とmcMap()を比べる。

x <- data.frame(strA = c("A", "B", "C", "D", "E", "F"),
strB = c("a", "b", "c", "d", "e", "f"),
strC = c("1", "2", "3", "4", "5", "6"))

for(i in 1:16) x <- rbind(x, x)

> nrow(x)
[1] 393216

Map()でやると、16.926秒かかる。for()は怖くて試せない。

system.time(
x$strABC <- unlist(Map(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC))
)

ユーザシステム経過
16.596 0.107 16.926

mcMap(mc.cores = 4)でやると、8.939秒。Map()の半分の時間。

system.time(
x$strABC <- unlist(mcMap(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC, mc.cores = 4))
)

ユーザシステム経過
1.104 0.172 8.939

mcMap(mc.cores = 20)でやると、3.412秒。かなり速い。

system.time(
x$strABC <- unlist(mcMap(f = get.strABC,
strA = x$strA,
strB = x$strB,
strC = x$strC, mc.cores = 20))
)

ユーザシステム経過
80.537 3.933 3.412

R: read.table() EOFの警告

2016-07-07T00:00:00.000+09:00

read.table()で以下のような警告が出る場合、read.table()に「quote = ""」を与えると、とりあえず警告が消える。
When you are using the function read.table() and encounter the warnings like below, just give 'quote = ""' to the read.talbe() to avoid the problem.

【日本語の警告】
警告メッセージ:
scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, で:
クオートで囲まれた文字列中にEOFがあります

【Warning in English】
Warning message:In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, :
EOF within quoted string

上は、「"」や「'」が悪さをしているときの解決法になります。

もともと "\"'" がデフォルトなので、取り込み対象の書類に２カ所のダブルクォートあるいはシングルクォートがあると、警告なしに、一見問題なく取り込まれます。

異なる２つの行に個別に奇数のクォートがあると、その間の行はトリミングされます。しかも恐ろしいことに、警告がでません。

EOF のメッセージが出るのは、たまたま奇数のクォートが奇数行に存在するときだけらしいです。

R：文字列をつなぐ

2016-05-27T00:00:00.000+09:00

Rで文字を繋ぐためにはpaste()を使う。
You can concatenate character strings with the function paste().

ひとつの文字列ベクターに含まれる各文字列を繋いで１つにしたいとき、collapseを使うと良い。
If you want to concatenate the strings in a character vector to obtain one string in which all the strings in the vector are linked together, use the argument collapse.

x <- "ABC"
y <- "DEF"
z <- c(x, y)

# zには２つの文字列が入っている。
# There are two strings in the character vector z.
# > z
# [1] "ABC" "DEF"

# collapseを使わないとき、文字列は２つのまま、文字列が繋がる。
# When you use paste() without the collapse argument, the resultant object will contain two strings.
# > paste(z, sep = "")
# [1] "ABC" "DEF"

# collapseを使うと、文字列が１つにまとまる。
# If you use the collapse, the strings will be linked together.
# > paste(z, sep = "", collapse = "")
# [1] "ABCDEF"

# つなぎを指定することもできる。
# You can add your favorite characters in between the strings.
# > paste(z, sep = "", collapse = "*")
# [1] "ABC*DEF"

Fortranにやってもらう

2016-04-21T00:00:00.000+09:00

Rをある程度習熟すれば、考えたことをなんでも実行できる。塩基配列の取扱いはBiostringsパッケージが強力だし、parallelパッケージの並列処理で多くの作業を高速化できる。

だけれども、やりたい計算があまりにもマニアックでapply系が使えないとき、とても困る。for()を多用せざるをえないような計算だ。

一応、for()も高速化はできる。先に適切なオブジェクトをつくっておき、その中のある座標だけを入れ替えるようにfor()ループを組み立てておけば、それなりに速く計算が終わる。宣言する感じ。そのつどオブジェクトを作り直すように組み立ててしまうと極端に遅い。それに、巨大なデータを用いて計算すると、速いといっても大したことない。

そこでFortranの出番となる。

RとFortranをつなぐのは.Fortran()という関数。Fortranのコードはプログラムではなくサブルーチンとして用意する。RからFortranに変数を渡せるし、FortranからRに変数を返してもらえる。ただし、サブルーチンの記述通りの順番で変数を渡す必要があるし、渡した変数しか返してもらえない。

サブルーチンの最後まで行った時点で保持していた変数のうち、Rから受け取ったものと同じもののみが、Rに帰ってくる。リストとして束ねられている。Fortranサイドで得た計算結果を受け取るためには、あらかじめRから初期値を入れた変数を与えて、それをFortranサイドで上書きしてしまえば良い。

RからFortranに変数が渡されると、Fortranが計算する。Fortranから変数が戻ってくるとき、Rのオブジェクトの中身だけを入れ替える形で戻ってくるらしい。だから、Fortranサイドでの宣言の仕方によってはおかしなことになる。

たとえば character(3) x みたいな感じで３個しか文字を格納できないようにしていると、Rから「x = "ABCDE"」というふうに５文字を受け取ったとしても Fortranで使うのは最初の３文字だけとなる（"ABC"）。この３文字を仮に別の文字に入れ替えてみるとする。たとえば、「x(3:3) = 'a'」とすると"ABa"になるのだが、このxをRに返してやると、上書きで返すのでRが受け取る変数は最初の３文字だけが改変された５文字、つまり"ABaDE"ということになる。

数字の並びを部分的に置換したいとき

2016-04-01T00:00:00.000+09:00

数字の並びを部分的に置換することを考える。

たとえば、あるinteger vectorの中に存在する「3, 4, 5」を「5, 4, 3」に置き換えたいとき。
ここでは、文字列用の関数であるpaste()とgsub()を使ったトリックを紹介する。

適当なinteger vectorを作る。
int <- 0:9
int.vec <- c(int, int, int)

> int.vec
[1] 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9

【for()ループでやる場合】

int.vec.new <- int.vec
for(i in 1:length(int.vec)){
if(int.vec[i] == 3 & int.vec[i+1] == 4 & int.vec[i+2] == 5){
int.vec.new[i:(i+2)] <- c(5, 4, 3)
i <- i + 2
}

> int.vec.new
[1] 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9

int.vecの中の「3, 4, 5」を「5, 4, 3」に置き換える仕事をfor()ループでやろうとすると、
このように一応はできるのだけれども、for()ループは遅い。
とくにvectorが長いほど遅さを実感する。

そこで、paste()とgsub()を使ったトリックの登場。

【paste()とgsub()でやる場合】

chr.vec <- paste(int.vec, sep = "", collapse = "")
int.vec.new <- as.integer(strsplit(gsub("345", "543", chr.vec), split = "")[[1]])

> int.vec.new
[1] 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9

まず、int.vecの値を文字列として繋げてしまう。
paste() はintegerを与えられたとしても、それをcharacterとして扱う。
collapseをうまく使えばすべての文字が繋がって、値がひとつのcharacter vectorになる。

chr.vec <- paste(int.vec, sep = "", collapse = "")

> chr.vec
[1] "012345678901234567890123456789"

こうしてひとつなぎの文字列にしてしまえば、gsub()が使える。
sub()は１つ目のヒットだけを置換するが、gsub()はすべてのヒットを置換する。

gsub("345", "543", chr.vec)

> gsub("345", "543", chr.vec)
[1] "012543678901254367890125436789"

あとは、strsplit()で１文字ずつ切り出し、integerに戻す。
strsplit() はリストで返すので、その１つ目を採用する。

> strsplit(gsub("345", "543", chr.vec), split = "")
[[1]]
[1] "0" "1" "2" "5" "4" "3" "6" "7" "8" "9" "0" "1" "2" "5" "4" "3" "6" "7" "8"
[20] "9" "0" "1" "2" "5" "4" "3" "6" "7" "8" "9"

回りくどい仕事のように見えるが、対象が大きくなるほど早さを実感する。

あるドメインをもつタンパク質をコードする遺伝子の数

2016-03-22T00:00:00.000+09:00

興味のあるドメインがあるとして、そのドメインをもつタンパク質がヒトゲノムにどれだけコードされているかを調査する方法のメモ。

例えば、HMGというドメイン。このドメインの「SMART ID」は「SM00398」。
http://smart.embl.de/smart/do_annotation.pl?BLAST=DUMMY&DOMAIN=SM00398

ポイントは「SMART MODE」を「GENOMIC」にしておくこと。
SMARTのサイトのEvolutionという項目をたどればヒトにいくつのタンパク質があるかわかる。
「SMART MODE」が「NORMAL」だと147。
「SMART MODE」が「GENOMIC」だと37。

この37の数字をクリックして「Expand all nodes」をクリックするとすべての候補の名前を確認することができる。
でも、「ENSP...」というEnsemblのIDが付与されたものや、「..._HUMAN」という感じで遺伝子名が付与されたものや、UniprotのIDが付与されたものが混在している。だから、結局のところ現行ゲノムにどれだけ遺伝子が存在するのかよくわからない。そもそも「GENOMIC」での絞り込みが妥当かどうかがよくわからない。

この問題を回避するため、EnsemblのBiomartを使う。
http://www.ensembl.org/biomart/martview

セッティングは以下のとおり。
Dataset:
CHOOSE DATABASE: Ensembl Genes 84
CHOOSE DATASET: Homo sapience genes (GRCh38.p5)
Filters:
PROTEIN DOMAINS AND FAMILIES
Limit to genes with these family or domain IDs
SMART ID(s): SM00398
Attributes:
Features
GENE: Associated Gene Name, Description
PROTEIN DOMAINS AND FAMILIES: Domains: SMART ID

ポイントは、Ensemble Gene IDなどを呼び出さないこと。これをやっちゃうと変なものまで出てくる。
Countすると、52件のヒット。
Resultsの「Unique results only」で限定してダウンロードすると、52遺伝子の情報を得ることができる。

こうすれば、特定のSMART IDをもつ遺伝子が現行ゲノムにいくつあるかを知ることができる。
Transcript variantがやたらとあるので遺伝子単位で数を数えるのが妥当だと思う。

R: new.env() で環境をつくる

2013-08-28T01:04:00.000+09:00

CRAN のポリシーのひとつに、大局的環境（.GlobalEnv）を操作してはならない、というのがある。自作パッケージの関数の中でオブジェクトを .GlobalEnv に assgin() したり、そこから get() していたので、このポリシーの存在には正直へこんだ。少なくとも、最終段階で気づくべきことではない。

google 先生に色々訊ねた結果、どうやらこの道の常識では .GlobalEnv は綺麗にしておくものらしい。

ポリシーに従ったパッケージにするためには、パッケージの環境に new.env() で小環境を作れば良い。たとえば、DATA という環境をつくる。そして、パッケージの中の assign(pos = .GlobalEnv) や get (pos = .GlobalEnv) を、assign(pos = DATA) や get(pos = DATA) に置換すればうまくいくはず。

以下のコードを R コンソールにコピペして、各関数を実行してみたら理解しやすいと思う。ls() で調べれば、.GlobalEnv の中に小環境 DATA が作られるのがわかる。ls(DATA) で小環境 DATA の中のオブジェクトを確認できる。

# create an environment named DATA.

DATA <- new.env(parent = emptyenv())

# assign an object to DATA.

assignToDATA <- function (name = "unknown", value = "unknown"){

 assign(name, value, pos = DATA)

}

# get an object from DATA.

getFromDATA <- function (name = "unknown"){

 get(name, pos = DATA)

}

# remove an object from DATA.

rmFromDATA <- function (list = "unknown"){

 rm(list = list, pos = DATA)

}

パッケージの「R」ディレクトリに入れる「.r」ファイルの１つに「DATA <- 」の行を書き込んでおき、NAMESPACE に export("DATA") を加えておけば、パッケージをインストールして library() でロードしたときに、パッケージの環境中に小環境として DATA が作られる。これなら .GlobalEnv を操作しないで済む。DATA の中身は ls(DATA) で閲覧できる。

これはわかったけど、パッケージ中の至る所に get や assign があるので、これらを全部修正するのは非常に面倒。最初に知っておきたかった。

メモ：CRAN にパッケージを投稿する際の問題

2013-08-26T00:00:00.000+09:00

メモです。

パッケージのチェックには「--as-cran」を使う。
R CMD check --as-cran test.package

Mac の Finder が自動生成する隠しファイル（.DS_Store）の存在が邪魔。これは rm で取り除くことができるから OK。

* checking for hidden files and directories ... NOTE
Found the following hidden files and directories:
.DS_Store
These were most likely included in error. See section ‘Package
structure’ in the ‘Writing R Extensions’ manual.

assign で .GlobalEnv にオブジェクトを作ったらダメ。この解決法を知りたい。別の環境をつくってそこに放り込む？

* checking R code for possible problems ... NOTE

Found the following assignments to the global environment:

File ‘test.package/R/assign.obj.R’:

assign(name, value, pos = .GlobalEnv)

R: やっぱりapply系を使い慣れておかないと

2012-12-23T01:49:00.000+09:00

Rはベクトル計算が強い。ていうかforループの多用は格好悪い。
いっぱい関数つくったりする前に、やっぱりapply系を使いこなしておきたい。

ということで、apply()の練習。

apply()は行列（matrix）の行（1）あるいは列（2）を対象として関数を適用する。

apply(X,MARGIN,FUN,OA)

　　X: 対象とする行列
　　MARGIN: 1 あるいは 2
　　FUN: 関数
　　OA: FUNに与える X 以外の引数（必要時のみ）

まず、適当な行列をつくる

set.seed(0) # 乱数のseedを規定
numbers <- rnorm(12,0,1)
matrix1 <- matrix(numbers,nrow=4)

> matrix1
[,1] [,2] [,3]
[1,] 1.2629543 0.4146414 -0.005767173
[2,] -0.3262334 -1.5399500 2.404653389
[3,] 1.3297993 -0.9285670 0.763593461
[4,] 1.2724293 -0.2947204 -0.799009249

行単位でmax()。

> apply(matrix1,1,max)

[1] 1.262954 2.404653 1.329799 1.272429

列単位でmax()。

> apply(matrix1,2,max)
[1] 1.3297993 0.4146414 2.4046534

自作関数をつくってみる。関数名はf1220.1()。
max()で求めたxの最大値にyを足す。

f1220.1 <- function(x,y){
z <- max(x) + y
return(z)
}

関数の動作確認。1, 2, 3の最大値に1を足すと4。

num123 <- c(1,2,3)

> f1220.1(num123,1)
[1] 4

apply()にはFUNのoptional argumentsを与えることができる。
行あるいは列単位でf1220.1()を実行。引数は1。

> apply(matrix1,1,f1220.1,1)
[1] 2.262954 3.404653 2.329799 2.272429

> apply(matrix1,2,f1220.1,1)
[1] 2.329799 1.414641 3.404653

引数は複数でもOK。

f1220.2 <- function(x,y,a,b){
z <- max(x) + y + a + b
return(z)
}

> apply(matrix1,1,f1220.2,1,1,1)
[1] 4.262954 5.404653 4.329799 4.272429

> apply(matrix1,2,f1220.2,1,1,1)
[1] 4.329799 3.414641 5.404653

apply()の対象となる行列に、行と列の名前がある場合。

colnames(matrix1) <- c("A","B","C")
rownames(matrix1) <- c("first.row","second.row","third.row","fourth.row")

> matrix1
A B C
first.row 1.2629543 0.4146414 -0.005767173
second.row -0.3262334 -1.5399500 2.404653389
third.row 1.3297993 -0.9285670 0.763593461
fourth.row 1.2724293 -0.2947204 -0.799009249

結果はこんな風になる。

> apply(matrix1,1,max)
first.row second.row third.row fourth.row
1.262954 2.404653 1.329799 1.272429

> apply(matrix1,2,max)
A B C
1.3297993 0.4146414 2.4046534

次に、lapplyとsapplyの練習。
まず、適当なリストをつくる。

set.seed(0)
a <- rnorm(10,0,1)
set.seed(1)
b <- rnorm(10,0,1)
list1 <- list(a,b)

> list1
[[1]]
[1] 1.262954285 -0.326233361 1.329799263 1.272429321
[5] 0.414641434 -1.539950042 -0.928567035 -0.294720447
[9] -0.005767173 2.404653389

[[2]]
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
[6] -0.8204684 0.4874291 0.7383247 0.5757814 -0.3053884

lapplyはリストを対象とする。リストの各要素に対して関数を適用する。

> lapply(list1,mean)
[[1]]
[1] 0.358924

[[2]]
[1] 0.1322028

sapplyはlapplyの結果をsimplifyして返す。

> sapply(list1,mean)
[1] 0.3589240 0.1322028

自作関数も動く。

f1220.3 <- function(x){

y <- x + 1

return(y)

}

> lapply(list1,f1220.3)

[[1]]

[1] 2.26295428 0.67376664 2.32979926 2.27242932 1.41464143

[6] -0.53995004 0.07143297 0.70527955 0.99423283 3.40465339

[[2]]

[1] 0.3735462 1.1836433 0.1643714 2.5952808 1.3295078 0.1795316

[7] 1.4874291 1.7383247 1.5757814 0.6946116

> sapply(list1,f1220.3)

[,1] [,2]

[1,] 2.26295428 0.3735462

[2,] 0.67376664 1.1836433

[3,] 2.32979926 0.1643714

[4,] 2.27242932 2.5952808

[5,] 1.41464143 1.3295078

[6,] -0.53995004 0.1795316

[7,] 0.07143297 1.4874291

[8,] 0.70527955 1.7383247

[9,] 0.99423283 1.5757814

[10,] 3.40465339 0.6946116

引数を必要とする関数も動かせる。

f1220.4 <- function(x,b){

y <- x + b

return(y)

}

> lapply(list1,f1220.4,2)

[[1]]

[1] 3.262954 1.673767 3.329799 3.272429 2.414641 0.460050 1.071433

[8] 1.705280 1.994233 4.404653

[[2]]

[1] 1.373546 2.183643 1.164371 3.595281 2.329508 1.179532 2.487429

[8] 2.738325 2.575781 1.694612

> sapply(list1,f1220.4,2)

[,1] [,2]

[1,] 3.262954 1.373546

[2,] 1.673767 2.183643

[3,] 3.329799 1.164371

[4,] 3.272429 3.595281

[5,] 2.414641 2.329508

[6,] 0.460050 1.179532

[7,] 1.071433 2.487429

[8,] 1.705280 2.738325

[9,] 1.994233 2.575781

[10,] 4.404653 1.694612

リストの中のリストも対象にできる。

まず、入れ子状のリストをつくる。

names <- c("a","b","c","d","e")

set.seed(0)

value.a <- rnorm(10,0,1)

set.seed(1)

value.b <- rnorm(10,0,1)

set.seed(2)

value.c <- rnorm(10,0,1)

set.seed(3)

value.d <- rnorm(10,0,1)

set.seed(4)

value.e <- rnorm(10,0,1)

list2 <- list()

list2$names <- names

values <- list(a=value.a, b=value.b, c=value.c, d=value.d, e=value.e)

list2$values <- values

> list2

$names

[1] "a" "b" "c" "d" "e"

$values

$values$a

[1] 1.262954285 -0.326233361 1.329799263 1.272429321

[5] 0.414641434 -1.539950042 -0.928567035 -0.294720447

[9] -0.005767173 2.404653389

$values$b

[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078

[6] -0.8204684 0.4874291 0.7383247 0.5757814 -0.3053884

$values$c

[1] -0.89691455 0.18484918 1.58784533 -1.13037567 -0.08025176

[6] 0.13242028 0.70795473 -0.23969802 1.98447394 -0.13878701

$values$d

[1] -0.96193342 -0.29252572 0.25878822 -1.15213189 0.19578283

[6] 0.03012394 0.08541773 1.11661021 -1.21885742 1.26736872

$values$e

[1] 0.2167549 -0.5424926 0.8911446 0.5959806 1.6356180

[6] 0.6892754 -1.2812466 -0.2131445 1.8965399 1.7768632

クラスを確認。

> class(list2)

[1] "list"

> class(list2$values)

[1] "list"

入れ子のリストにlapply()を適用。

> lapply(list2$values,mean)

[1] 0.358924

[1] 0.1322028

[1] 0.2111516

[1] -0.06713568

[1] 0.5665293

sapply()も動く。

> sapply(list2$values,mean)
a b c d e
0.35892396 0.13220278 0.21115165 -0.06713568 0.56652929

リストの中の行列にはapply()が使える。
まず、行列を含むリストをつくる。

list3 <- list()
list3$names <- names
values <- cbind(a=value.a, b=value.b, c=value.c, d=value.d, e=value.e)
rownames(values) <- seq(1,10,1)
list3$values <- values

> list3
$names
[1] "a" "b" "c" "d" "e"

$values
a b c d e
1 1.262954285 -0.6264538 -0.89691455 -0.96193342 0.2167549
2 -0.326233361 0.1836433 0.18484918 -0.29252572 -0.5424926
3 1.329799263 -0.8356286 1.58784533 0.25878822 0.8911446
4 1.272429321 1.5952808 -1.13037567 -1.15213189 0.5959806
5 0.414641434 0.3295078 -0.08025176 0.19578283 1.6356180
6 -1.539950042 -0.8204684 0.13242028 0.03012394 0.6892754
7 -0.928567035 0.4874291 0.70795473 0.08541773 -1.2812466
8 -0.294720447 0.7383247 -0.23969802 1.11661021 -0.2131445
9 -0.005767173 0.5757814 1.98447394 -1.21885742 1.8965399
10 2.404653389 -0.3053884 -0.13878701 1.26736872 1.7768632

クラスを確認。

> class(list3)
[1] "list"

> class(list3$values)
[1] "matrix"

apply()を動かす。

> apply(list3$values,1,max)
1 2 3 4 5 6
1.2629543 0.1848492 1.5878453 1.5952808 1.6356180 0.6892754
7 8 9 10
0.7079547 1.1166102 1.9844739 2.4046534

> apply(list3$values,2,max)
a b c d e
2.404653 1.595281 1.984474 1.267369 1.896540

こういうのをちゃんと理解してから複雑な仕事をした方が良い。

R: グラフの誤差範囲（標準偏差）をグレー表示

2012-10-14T22:50:00.001+09:00

Rでグラフを作る際、エラーバーの代わりにグレーで視覚化する。
例えば、５サンプルのグラフ。

標準偏差をグレー表示

# インデックスは10。

len <- 10

# 5サンプルのデータを適当につくる。平均値は10。

a <- rnorm(len,10,1)

b <- rnorm(len,10,1)

c <- rnorm(len,10,1)

d <- rnorm(len,10,1)

e <- rnorm(len,10,1)

# 行列にする

main.table <- rbind(a,b,c,d,e)

# 平均と標準偏差を計算する

main.mean <- apply(main.table,2,mean)

main.sd <- apply(main.table,2,sd)

# ポリゴンのx座標とy座標をつくる

x <- 1:len

poly.x <- c(x,rev(x))

poly.y <- c(main.mean+main.sd,rev(main.mean-main.sd))

# プロットする。最初はデータを入れずy軸はエラー分の余裕をつくる。

# ポリゴンを入れる。グレーで埋めて、輪郭線は消す。

# 線をプロットする。

plot(main.mean,type="n",ylim=c(min(main.mean)-max(main.sd),max(main.mean)+max(main.sd)),)

polygon(poly.x,poly.y,col="gray",lty=0)

lines(main.mean,lwd=2)

# 別ウィンドウで個別データをプロット

quartz("each")

plot(main.table[a],type="n",ylim=c(min(main.mean)-max(main.sd),max(main.mean)+max(main.sd)),)

lines(main.table[a])

lines(main.table[b])

lines(main.table[c])

lines(main.table[d])

lines(main.table[e])

R: wigもつくる

2012-10-14T02:09:00.001+09:00

せっかく自動でgffをつくったので、wigもつくってみる。
まず、前回の「gff.gen()」に若干の修正を加えてgffを得る。

gffをつくる関数

gff.gen <- function(n=1000,chrnum=3,type="protein_coding"){

chrs <- paste("chr",1:chrnum,sep="")

source <- "R"

type <- type

dots <- rep(".",n)

chr <- sample(chrs,n,replace=T)

start <- sample(1:9990000,n)

end <- start+1500+(round(rnorm(n)/2*500))

strand <- sample(c("+","-"),n,replace=T,prob=c(0.5,0.5))

id <- paste("gene",1:n,sep="")

gff <- data.frame(chr,source,type,start,end,dots,strand,dots,id)

return(gff)

}

gffをつくってソート

gff <- gff.gen(1000,3,"protein_coding")

gff <- sort.left(gff)

このgffの中の一部の遺伝子の5'末端にふくらみのあるwigをつくる。

# wigをつくる関数（上のgffがあることが前提）

wig.gen <- function(chrnum=3){

# chrの名前をつくる

chrs <- paste("chr",1:chrnum,sep="")

# chr毎に作業する

for (i in 1:length(chrs)){

chr <- chrs[i]

# headerをつくる

wig.head1 <- paste("track type=wiggle_0 name=\"",chr,"\" description=\"generated by R\"",sep="")

wig.head2 <- paste("variableStep chrom=",chr," span=10",sep="")

wig.head <- c(wig.head1,wig.head2)

wig.head.file <- paste("head_",chr,".txt",sep="")

write(wig.head,file=wig.head.file)

# bodyをつくる

n <- 10000000

step <- 10

pos <- seq(1,n,step)

depth <- round(rnorm(1:length(pos),100,sd=2))

wig.body <- cbind(pos,depth)

# chr毎のgffをつくる

assign(paste("gff.",chr,".fw",sep=""),gff[gff$chr==chr&gff$strand=="+",])

assign(paste("gff.",chr,".rv",sep=""),gff[gff$chr==chr&gff$strand=="-",])

# bodyに変化を加える

y.rep <- nrow(get(paste("gff.",chr,".fw",sep="")))%/%5

for(y in 1:y.rep) {

tss <- sample(get(paste("gff.",chr,".fw",sep=""))[,4],1)%/%10

upst <- tss-19

dnst <- tss+20

hight <- sample(1:2000,1)

add.tag.upst <- round((cos(seq(1*pi,0,length=20))+1)*hight)

add.tag.dnst <- round((cos(seq(0,1*pi,length=20))+1)*hight)

wig.body[upst:tss,2] <- wig.body[upst:tss,2]+add.tag.upst

wig.body[(tss+1):dnst,2] <- wig.body[(tss+1):dnst,2]+add.tag.dnst

tss <- sample(get(paste("gff.",chr,".rv",sep=""))[,5],1)%/%10

upst <- tss-19

dnst <- tss+20

hight <- sample(1:2000,1)

add.tag.upst <- round((cos(seq(1*pi,0,length=20))+1)*hight)

add.tag.dnst <- round((cos(seq(0,1*pi,length=20))+1)*hight)

wig.body[upst:tss,2] <- wig.body[upst:tss,2]+add.tag.upst

wig.body[(tss+1):dnst,2] <- wig.body[(tss+1):dnst,2]+add.tag.dnst

}

# bodyを保存する

wig.body.file <- paste("body_",chr,".txt",sep="")

write.table(wig.body,file=wig.body.file,col.names=F,row.names=F,sep=" ")

# headerとbodyをつなげる

wig.file <- paste(chr,".wig",sep="")

sys.command <- paste("cat",wig.head.file,wig.body.file,">",wig.file,sep=" ")

system(sys.command)

# headerとbodyを消す

system(paste("rm",wig.head.file,sep=" "))

system(paste("rm",wig.body.file,sep=" "))

} # chr毎のループはここまで

} # この関数はここまで

上の関数を使ってwigをつくる

wig.gen(3)

R: gffをつくる

2012-10-12T22:31:00.001+09:00

いくつかの解析系のことを考えていると、gffフォーマットのファイルが欲しくなった。Rで自動生成してみる。

gffファイルをつくる関数

行数と染色体の数、それとタイプを指定してgffの中身をランダムに生成する。

# n: 行数

# chrnum: 染色体の数

# type: "protein_coding", etc.

gff.gen <- function(n,chrnum,type){

chrs <- paste("chr",1:chrnum,sep="")

source <- "R"

type <- type

dots <- rep(".",n)

chr <- sample(chrs,n,replace=T)

start <- abs(round(rnorm(n)*1000000))

end <- start+1500+(round(rnorm(n)/2*500))

strand <- sample(c("+","-"),n,replace=T,prob=c(0.5,0.5))

id <- paste("gene",1:n,sep="")

gff <- data.frame(chr,source,type,start,end,dots,strand,dots,id)

return(gff)

}

# usage

gff <- gff.gen(1000,3,"ncRNA")

# check

head(gff)

# min length

min(gff[,5]-gff[,4])

# distribution

hist(gff[,5]-gff[,4],main="gene length")

# sort.left

gff.s <- sort.left(gff)

# check

head(gff.s)

# number of each strand

nrow(gff.s[gff.s$strand=="+",])

nrow(gff.s[gff.s$strand=="-",])

関数「sort.left()」はこちらを参照。

上の「# usage」以下の実行結果。

> # usage

> gff <- gff.gen(1000,3,"ncRNA")

> # check

> head(gff)

chr source type start end dots strand dots.1 id

1 chr2 R ncRNA 272400 273456 . - . gene1

2 chr1 R ncRNA 1173705 1175235 . + . gene2

3 chr1 R ncRNA 726009 727790 . + . gene3

4 chr3 R ncRNA 2041368 2042869 . - . gene4

5 chr3 R ncRNA 607172 608515 . + . gene5

6 chr3 R ncRNA 556786 558183 . + . gene6

> # min length

> min(gff[,5]-gff[,4])

[1] 767

> # distribution

> hist(gff[,5]-gff[,4],main="gene length")

histの結果

> # sort.left

> gff.s <- sort.left(gff)

> # check

> head(gff.s)

chr source type start end dots strand dots.1 id

150 chr1 R ncRNA 1017 2463 . - . gene150

378 chr1 R ncRNA 9260 11146 . - . gene378

164 chr1 R ncRNA 9922 10844 . + . gene164

93 chr1 R ncRNA 14776 15981 . + . gene93

557 chr1 R ncRNA 17269 18323 . + . gene557

863 chr1 R ncRNA 21640 22834 . - . gene863

> # number of each strand

> nrow(gff.s[gff.s$strand=="+",])

[1] 495

> nrow(gff.s[gff.s$strand=="-",])

[1] 505

R: 左側のカラムを優先して全カラムをソートする

2012-10-12T12:22:00.003+09:00

PeakAnnotatorを使ってみたら、順番に並んでないから無理と言われた。

こんな結果。

$ PeakAnnotator TSS peaks.bed genes.bed out/tss.txt

Starting...

Reading BED File

Files must be sorted by chromosomes,start position, end position

ということで、左のカラムを優先して全カラムをソートする関数を作った。
行列xに適用すれば並べ替えてくれる。普遍的に使える。列数の限界なし。

sort.left <- function(x,decreasing=FALSE){

n <- ncol(x)

colnum <- n

for (i in 1:n){

x <- x[order(x[,colnum],decreasing=decreasing),]

colnum <- colnum-1

}

return(x)

}

以下、挙動の確認

４列のマトリックスに１から８をランダムに入れた。

> x <- matrix(sample(1:8,4*16,replace=T),ncol=4)

xの中身を確認。

> x

[,1] [,2] [,3] [,4]

[1,] 2 5 5 2

[2,] 3 4 8 2

[3,] 5 6 8 1

[4,] 2 6 2 4

[5,] 2 2 2 4

[6,] 5 6 4 7

[7,] 7 1 6 1

[8,] 7 1 8 7

[9,] 4 2 8 8

[10,] 8 4 2 6

[11,] 5 2 4 4

[12,] 5 1 4 5

[13,] 8 7 4 7

[14,] 4 5 7 5

[15,] 6 5 3 7

[16,] 6 5 5 6

関数「sort.left()」を適用。

１列目が同じ場合には２列目（１行目と２行目を参照）が、２列目まで同じ場合には３列目（９行目から１４行目参照）が、ちゃんと昇順になっているのがわかる。

> sort.left(x)

[,1] [,2] [,3] [,4]

[1,] 2 2 2 4

[2,] 2 5 5 2

[3,] 2 6 2 4

[4,] 3 4 8 2

[5,] 4 2 8 8

[6,] 4 5 7 5

[7,] 5 1 4 5

[8,] 5 2 4 4

[9,] 5 6 4 7

[10,] 5 6 8 1

[11,] 6 5 3 7

[12,] 6 5 5 6

[13,] 7 1 6 1

[14,] 7 1 8 7

[15,] 8 4 2 6

[16,] 8 7 4 7

引数の「decreasing」はデフォルトでFALSEだが、TRUEを与えると降順になる。

> sort.left(x,T)

[,1] [,2] [,3] [,4]

[1,] 8 7 4 7

[2,] 8 4 2 6

[3,] 7 1 8 7

[4,] 7 1 6 1

[5,] 6 5 5 6

[6,] 6 5 3 7

[7,] 5 6 8 1

[8,] 5 6 4 7

[9,] 5 2 4 4

[10,] 5 1 4 5

[11,] 4 5 7 5

[12,] 4 2 8 8

[13,] 3 4 8 2

[14,] 2 6 2 4

[15,] 2 5 5 2

[16,] 2 2 2 4

せっかくなので、右側優先の関数もつくった。

sort.right <- function(x,decreasing=FALSE){

n <- ncol(x)

colnum <- n

for (i in 1:n){

x <- x[order(x[,n-colnum+1],decreasing=decreasing),]

colnum <- colnum-1

}

return(x)

}

xに関数「sort.right()」を適用。

最右列から優先的にソートされる。

> sort.right(x)

[,1] [,2] [,3] [,4]

[1,] 7 1 6 1

[2,] 5 6 8 1

[3,] 2 5 5 2

[4,] 3 4 8 2

[5,] 2 2 2 4

[6,] 2 6 2 4

[7,] 5 2 4 4

[8,] 5 1 4 5

[9,] 4 5 7 5

[10,] 8 4 2 6

[11,] 6 5 5 6

[12,] 6 5 3 7

[13,] 5 6 4 7

[14,] 8 7 4 7

[15,] 7 1 8 7

[16,] 4 2 8 8

もちろん降順もできる。

> sort.right(x,T)

[,1] [,2] [,3] [,4]

[1,] 4 2 8 8

[2,] 7 1 8 7

[3,] 8 7 4 7

[4,] 5 6 4 7

[5,] 6 5 3 7

[6,] 6 5 5 6

[7,] 8 4 2 6

[8,] 4 5 7 5

[9,] 5 1 4 5

[10,] 5 2 4 4

[11,] 2 6 2 4

[12,] 2 2 2 4

[13,] 3 4 8 2

[14,] 2 5 5 2

[15,] 5 6 8 1

[16,] 7 1 6 1

列が増えても大丈夫。

> y <- matrix(sample(1:4,6*16,replace=T),ncol=6)

> sort.left(y)

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 1 1 3 4 3

[2,] 1 2 3 3 1 4

[3,] 1 3 1 1 1 2

[4,] 1 3 3 3 2 3

[5,] 1 4 1 3 4 1

[6,] 1 4 3 4 1 3

[7,] 2 1 3 4 2 3

[8,] 2 1 4 1 2 4

[9,] 2 1 4 2 4 4

[10,] 2 2 1 3 4 4

[11,] 2 3 1 4 3 3

[12,] 3 1 4 4 3 1

[13,] 4 3 1 1 4 3

[14,] 4 3 1 3 1 2

[15,] 4 4 3 2 1 3

[16,] 4 4 3 2 2 1

> sort.right(y,T)

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 2 2 1 3 4 4

[2,] 2 1 4 2 4 4

[3,] 2 1 4 1 2 4

[4,] 1 2 3 3 1 4

[5,] 1 1 1 3 4 3

[6,] 4 3 1 1 4 3

[7,] 2 3 1 4 3 3

[8,] 2 1 3 4 2 3

[9,] 1 3 3 3 2 3

[10,] 1 4 3 4 1 3

[11,] 4 4 3 2 1 3

[12,] 4 3 1 3 1 2

[13,] 1 3 1 1 1 2

[14,] 1 4 1 3 4 1

[15,] 3 1 4 4 3 1

[16,] 4 4 3 2 2 1

MEMEが見つけるモチーフの中身

2012-06-19T14:00:00.000+09:00

たとえリピートをマスキングしてから実行したとしても、MEMEは沢山のモチーフを提示してくれる。

ランダムにサンプリングしたプロモーター配列をMEMEに渡すと、MEMEが発見するモチーフの傾向がつかめる。

転写開始点付近のGCに富む領域、Alu配列の下流にあるA-richな領域などが目立つ。A-richな領域はリピートマスキングでもマスクされないのでどうしてもMEMEが拾ってしまう。

ハズレらしきモチーフは極端に長い傾向にある。特にGCに富む領域の場合、ひとつの配列に類似のモチーフが複数見つかる。もちろん、ハズレと断定できるわけではないのだが、このような傾向があることを知っておくのは重要だと思う。

PowerPointの低い解像度について

2012-06-19T13:53:00.002+09:00

macでKeynoteを使っていると、PowerPointが使い難くてたまらない。数ある問題のひとつ、「貼付けたPDFの画質が落ちる」という問題への対処法。せっかくPDFで貼付けているのに、どうしてガビガビにしてしまうのか。液晶プロジェクタで出力するなら高解像度は必要ないが、配布資料の画質が荒いのは困る。

使っているプレゼンソフト
Keynote'09
PowerPoint 2008

解決法
PowerPointの「環境設定」の「保存オプション」の「解像度の設定」で解像度を高める。「グラフィックファイルの圧縮」を解除する。

PDFをソースからコピーし、プレビューで「クリップボードから新規作成」し、PNGで保存する。このときのPNGの解像度が重要。解像度 300 ピクセル／インチ以上が望ましい。

.pngファイルをPowerPointに貼付ける。

Keynoteにはめ込んだPDFほどではないが、それなりにクリアになる。

for()ループを入れ子にして検索するより、Rならではの方法で

2012-06-15T22:31:00.002+09:00

IDと数値をもつ表があり、それとは別にIDのリストがあるとする。そのとき、リスト中のIDに対応する数値を表から抽出するにはどうしたらいいだろうか。

Rを使う。

for()ループでIDリストを１つずつ取り出し、そのループの中で更にfor()ループを使って表を１行ずつ取り出し、対応するIDをif()で探すという方法が思い浮かぶ。for()ループを入れ子構造にするやり方だ。方法を下に示す。

この入れ子構造でやってみるとものすごく遅い。２つ目のfor()ループは表の行数マイナス１行分も無駄に動いているわけだから、遅いのも頷ける。でも、perlだったらもっと速くに結果がでるような気がする。どうやらRはループが苦手らしい。

作業効率がとても悪いので、for()ループをひとつ減らして同じ結果を得る方法を考えた。if()で条件検討をせず、セルが一致する行をそのままRっぽく拾う。方法を下に示す。

ここに示す小さなデータではスピードの違いを実感できないが、数千行の表の中から数百のリストに合致する行を得ようとすれば、その違いは雲泥の差。前者ならかけっぱなしで散歩にでも行きたくなるが、後者ならあくびをしている間に終了する。

とりあえず、説明用にデータを用意する。

データフレーム（df.1）を用意する。

id1 <- c("A","B","C","D","E","F")
value1 <- rnorm(6,mean=3,sd=1)
df1 <- data.frame(id1,value1)

> df1
id1 value1
1 A 2.484163
2 B 1.539250
3 C 2.773143
4 D 1.041669
5 E 2.409677
6 F 2.532448

IDリスト（id2）を用意する。

id2 <- c("B","D","E")

> id2
[1] "B" "D" "E"

入れ子構造のやり方。
forループを２つ重ねて検索する。

found1 <- df1[0,]
for (i in 1:length(id2)){
for (j in 1:nrow(df1)){
if (id2[i]==df1[j,1]){
hit1 <- df1[j,]
found1 <- rbind(found1,hit1)
}
}
}

> found1
id1 value1
2 B 1.539250
4 D 1.041669
5 E 2.409677

Rの強みを活かすやり方。
forループをひとつだけにする。

found2 <- df1[0,]
for(x in 1:length(id2)){
hit2 <- df1[df1[,1]==id2[x],]
found2 <- rbind(found2,hit2)
}

> found2
id1 value1
2 B 1.539250
4 D 1.041669
5 E 2.409677

MEMEの前にリピートをマスク

2012-06-11T15:28:00.007+09:00

MEMEを使ってプロモーター中のモチーフ探しをすると、意味がありそうな、なさそうな、保存性の高いエレメントが多数ヒットすることがある。これらは通常とても長い。ゲノムにはAluなどのリピート配列が多く、それらが解析対象のプロモーター配列セットの中に沢山含まれてしまっているためだ。

MEMEでDNAモチーフを探すときは、あらかじめリピートをマスクしたFASTAを用いるべきだろう。

リピートのマスキングには、giri（Genetic Information Research Institute）のRepeat Maskingツール「CENSOR」が使える。生物種（Sequence source）を指定し、Report simple repeatsをチェックし、FASTAを与えて実行すると、リピート部分の塩基を「X」に置き換えたFASTAを返してくれる。もちろん、どのようなリピートがどこにヒットするかも教えてくれる。

都合のいいことに、CENSORが返してくれたFASTA（「X」でリピートがマスクされたもの）は、そのままMEMEに与えることができる。

.wigのtag countをmedianでノーマライズする

2012-06-07T19:02:00.001+09:00

MACSが出力したwiggleのtag countや、samtools mpileup -Dあるいはsamtools depthで得られたdepthを比較解析したいとき、ノーマライズのひとつのやり方として、medianで割るという方法が考えられる。

Rを使えば簡単。

wiggleを読み込む。ヘッダー２行は省く（skip=2）。

wig <- read.table("<in.wig>",skip=2)

> head(wig)

V1 V2

1 1 3

2 11 34

3 21 49

4 31 65

5 41 81

6 51 105

２列目のデータを、２列目のmedianで割った値で置き換える。

wig$V2 <- wig$V2/median(wig$V2)

> head(wig)

V1 V2

1 1 0.007575758

2 11 0.085858586

3 21 0.123737374

4 31 0.164141414

5 41 0.204545455

6 51 0.265151515

小数点以下２桁にする

wig$V2 <- round(wig$V2,digits=2)

> head(wig)

V1 V2

1 1 0.01

2 11 0.09

3 21 0.12

4 31 0.16

5 41 0.20

6 51 0.27

まとめてやるなら、

wig <- read.table("<in.wig>",skip=2)

wig$V2 <- round(wig$V2/median(wig$V2),digits=2)

必要に応じてread.tableで書き出したり、プロットしたりする。

IGVの調整によりwiggle等を開けるようにする

2012-06-07T18:36:00.004+09:00

IGVの調整によりwiggleを開けるようにする

MACSが作成したwiggle（.wig）のヘッダー情報にある染色体名（chrom=）がIGVで用いるアセンブリと一致しない場合、wiggleのヘッダーを書き換える以外に、IGVの設定を調整するという解決策もある。

macなら、ホームフォルダに「igv」というフォルダがあり（IGV 2.0のときは不可視化されていた？）、その中の「genomes」に、エイリアスファイルなるものを作って入れてやるだけで良い。

詳しい説明はこちら。

Creating a Chromosome Name Alias File

.wigと染色体名が一致しない問題は、BWA等によるアライメントの時に使用するリファレンスゲノム（.fasta）の染色体名をあらかじめIGVに合わせておくことで回避できるが、様々なソースから得たGFFやBEDなどもまた、染色体名の記述方法が異なる場合がある。その都度書類を書き換えるより、上記のごとくIGVの設定を変更しておくと楽。