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 

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

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

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もつくる

せっかく自動で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をつくる

いくつかの解析系のことを考えていると、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: 左側のカラムを優先して全カラムをソートする

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が見つけるモチーフの中身

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

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

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

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

PowerPointの低い解像度について

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


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


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


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


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


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

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

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の前にリピートをマスク

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でノーマライズする

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等を開けるようにする

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の設定を変更しておくと楽。

IGVでwiggleを開く

MACSが解析したChIP-seqのデータ（.wig）をIGVで眺める。

必要なデータ

　<in.wig.gz>　MACSで作成したwiggle（.gzを解凍してもしなくても）

IGV（Integrative Genomics Viewer）はJavaで動くゲノムブラウザ。IGVのダウンロードページから、Binary distributionをダウンロードする。

IGV_2.1.17.zip（14.6 MB）を解凍し、igv.jarを起動。

生物（アセンブリ）を指定する。染色体名を確認する。「chr1」や「I」など、アセンブリにより染色体名の記述の仕方が異なる。

MACSが作成した.wig.gzを解凍して.wigにし、エディタで開いてヘッダーの２行目の「chrom=」の後ろを確認する。ここがIGVの記述方法と一致している必要がある。一致していなければ書き換える。

染色体名が一致していれば、.wig.gzのままでもIGVで開くことができる。

もうひとつのやり方はこちらを参照。