R：文字列をつなぐ

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にやってもらう

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"ということになる。

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

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

たとえば、ある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"


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

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

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

例えば、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() で環境をつくる

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 にパッケージを投稿する際の問題

メモです。

パッケージのチェックには「--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系を使い慣れておかないと

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)

$a

[1] 0.358924

$b

[1] 0.1322028

$c

[1] 0.2111516

$d

[1] -0.06713568

$e

[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 

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