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