バイオ系のためのR覚書：Heatmap.2の使い方覚書

Rでheatmapを作るのに,gplotsパッケージのheatmap.2を使うのがつかいよい。
一番わかりやすいのが、Mannheimiaさんのブログである。

基本的には
> install.packages("gplots")
>library(gplots)

で読み込んで、お好みのdataについてheatmapを書かせればよい。

>heatmap.2(as.matrix(data), col=greenred(75), scale="row", key=T, keysize=1.5,
density.info="none", trace="none",cexCol=0.9, cexRow=0.5)

個人的にはメモリのカラーはgreenredでRowも提示したほうがいいので、以上のようなコマンドが気に入っている。


あと入れるデーターについては、彼の説明がうまいのでいかに引用しておく。
また追記で補足します！

However this is not the way of dealing with the data. When we work with high throughput data, the first step is to log-transform the intensities and then apply a normalization method. We can do that with the following lines:

data = log2(data)
boxplot(data, cex.axis=0.8, las=2, main="Original distribution of data",
ylab="Log2(Intensity)") # Draw a boxplot.


# Normalization
library(preprocessCore)
data2 = normalize.quantiles(as.matrix(data)) # A quantile normalization.


# Copy the row and column names from data to data2:
rownames(data2) = rownames(data)
colnames(data2) = colnames(data)


boxplot(data2, cex.axis=0.8, las=2, main="Distribution after normalization",
ylab="Log2(Intensity)")


# t-test using the limma package:

library(limma)
design = cbind(Cell_A = c(1,1,1,0,0,0,0,0,0), # First 3 columns->Cell_A
Cell_B = c(0,0,0,1,1,1,0,0,0),
Cell_C = c(0,0,0,0,0,0,1,1,1)) # Last 3 columns->Cell_C
fit = lmFit(data2, design=design) # Fit the original matrix to the above design.
# We want to compare A vs. B, A vs. C and B vs. C
contrastsMatrix = makeContrasts("Cell_A-Cell_B","Cell_A-Cell_C","Cell_B-Cell_C",
levels = design)
fit2 = contrasts.fit(fit, contrasts = contrastsMatrix) # Making the comparisons.
fit2 = eBayes(fit2) # Moderating the t-tetst by eBayes method.

バイオ系のためのR覚書:How to make Venn diagrams by R(ベン図の書き方）

まずRでベン図を描いてみる (Making Venn diagrams)

データベースにディポジットされている各種遺伝子変異とそれに関連してCpGがHypermethylationされている遺伝子のリストDNAmet.csvをexcelで作成。

その後は
＞　x <- read.table(file="DNAmet.csv",header=T,sep=",")　
で読み込み

＞head(x)
# a b c d e
#1 PYROXD2 HPSE2 LOXL4 ABCC2 KAZALD1
#2 ENTPD7 ENTPD7 HPSE2 SH3PXD2A BTRC
#3 BLOC1S2 SCD ENTPD7 VWA2 FBXW4
#4 FAM178A FGF8 ABCC2 ADARB2 NOLC1
#5 FBXW4 ABLIM1 SCD CHST15 ELOVL3
#6 KCNIP2 ECHDC3 NDUFB8 DOCK1 SH3PXD2A

こんな感じのファイルである(カラム名のa,b,cは各遺伝子変異に相当）。要素の遺伝子数がそれぞれ異なるので、変異遺伝子ごとに遺伝子リストマトリックスを作成する。

> a <- as.matrix(x[,1])
> b <- as.matrix(x[,2][1:1135])
> c <- as.matrix(x[,3][1:4064])
> d <- as.matrix(x[,4][1:761])
> e <- as.matrix(x[,5][1:2450])

a,b,c遺伝子変異の間で共通するメチル化遺伝子のベン図をかかせるためには、

これをリスト型のオブジェクトに変え

> data <- list(a,b,c)

gplotsのvenn関数を利用すれば簡単にベン図が書ける

> library(gplots)
> venn(data)




４つ以上のベン図もVennでかける（making venn diagram with more than 4 elements）
同様に４つのベン図もVennでかける

> data2 <- list(a,b,c,d)
> venn(data2)



同様に5つのベン図もOK

> data3 <- list(a,b,c,d,e)
> venn(data3)



Venn図のデーター要素をとるためには？ (How to make a gene list from a venn diagram)

例えばこの図の真ん中の269に含まれるgene nameを出したい場合は、



intersect関数をつかって、a&b&cをt1に入れるといい。

> t1 <- intersect(a,intersect(b,c))
> t1 <- as.matrix(t1)

t1の列数（に含まれる遺伝子数）を見ると
>nrow(t1)
[1] 269
となっている。

実際には
>head(t1)

[1] "ENTPD7" "FAM53B" "DOCK1" "FAM196A" "EBF3" "BNIP3"

といった遺伝子が含まれている。

あとは
>write.table(t1,file="overlap.txt")
でテキストファイルに書き出せばエエやろ！

histgram by R

Rでhistgramは簡単に生成できるが、結構厄介な点もある。
まず基本的なhist()関数での生成を簡単に書いておく。

デモデーターを生成

> a <- rnorm(100,20,5)
> b <- rnorm(60, 10,3)
> c <- rnorm(80, 5,3)
> d <- rnorm(100,20,5)
> e <-rnorm(60,10,4)
> f <- rnorm(80,5,5)

まずはテストで
>hist(a)
と入力すると白黒のヒストグラムを簡単に作ってくれる。

これだと芸がないので、色をつけて、X軸のレンジとピッチサイズをbreaksで指定する。

>hist(a, col="blue", breaks=seq(-5,40,2.5))



同様にデーターb、ｃについても

＞hist(b, col="red", breaks=seq(-5,40,2.5))

で、


＞hist(c, col="green", breaks=seq(-5,40,2.5))
で

とここまでは簡単に生成できる。

次に重ね合わせ、二つ目以降の式でadd=Tとすればよいから、

>hist(a, col="blue", breaks=seq(-5,40,2.5))
>hist(b, col="red", breaks=seq(-5,40,2.5), add=T)
>hist(c, col="green", breaks=seq(-5,40,2.5), add=T)



Y軸のレンジをylim=c(1,40)で指定してちょい修正。

>hist(a, col="blue",ylim=c(1,40),breaks=seq(-5,40,2.5))
>hist(b, col="red", breaks=seq(-5,40,2.5), add=T)
>hist(c, col="green", breaks=seq(-5,40,2.5), add=T)



やはり重ね合わせはうまくいかない。
特に６つのグラフが重なるような状況になると悲惨である。

>hist(a, col="blue",ylim=c(1,40),breaks=seq(-15,40,2.5))
>hist(b, col="red", breaks=seq(-15,40,2.5), add=T)
>hist(c, col="green", breaks=seq(-15,40,2.5), add=T)
>hist(d, col="orange", breaks=seq(-15,40,2.5), add=T)
>hist(e, col="yellow", breaks=seq(-15,40,2.5), add=T)
>hist(f, col="purple", breaks=seq(-15,40,2.5), add=T)



とまあひどい感じの重ね合わせになる。

一応biostaticsのページにreshape2とggplot2をつかってやる綺麗にhistgramを重ね合わせる方法がでているのであるが、このページと同じ内容をおこなってもmelt()がうまく働いてくれない(*)

(*)なんとなくバグっぽい！！！

ANOVA by R

あまり詳しくないのですが、多変量解析もRでわりと簡単にできるのでやり方を書いておきます。フリーソフトでなんでもできるっていうのは、ある意味驚愕ですねぇ。

まずデモデーターをこれにて
#乱数で、サンプル数100, 平均20, 標準偏差5のデーターをつくる
> a <- rnorm(100,20,5)
#乱数で、サンプル数60, 平均10, 標準偏差3のデーターをつくる
> b <- rnorm(60, 10,3)
#乱数で、サンプル数80 平均5, 標準偏差3のデーターをつくる
> c <- rnorm(80, 5,3)
#乱数で、サンプル数100, 平均20, 標準偏差5のデーターをつくる
> d <- rnorm(100,20,5)
#乱数で、サンプル数60, 平均10, 標準偏差4のデーターをつくる
> e <-rnorm(60,10,4)
#乱数で、サンプル数80 平均5, 標準偏差5のデーターをつくる
> f <- rnorm(80,5,5)

ここではa-fはそれぞれ別々の実験条件、データーは例えば各細胞ごとのある遺伝子の発現とする（シングルセルの遺伝子解析みたいな）。

ANOVAをおこなうには、基本的には各実験条件a-fごとの数値データの複数のベクトルから、数値データの一続きのカラムとそれに対応した実験条件のカラムが入ったデーターシートをつくればよい。

詳しくはRにもともとはいっているChickWeightのデーター構造を参考にすればよい。
ChickWeightをみてみると

＞ChickWeight
weight Time Chick Diet
1 42 0 1 1
2 51 2 1 1
3 59 4 1 1
4 64 6 1 1
5 76 8 1 1
6 93 10 1 1
（中略)
574 175 14 50 4
575 205 16 50 4
576 234 18 50 4
577 264 20 50 4
578 264 21 50 4

とこのようになっており数値データー（weight)と実験条件（Diet)からなる２カラムのデーターフレームをつくるような感じである。

まずデーターカラムの作成をおこなう。
> 　z <- c(a,b,c,d,e,f)

次に実験条件カラムの作成をおこなう。a,b,c,d,e,fそれぞれのデーター数ごとに、実験条件を割り振ることが目標である。
まず各実験条件ごとのデーター数の把握。
＞　times <- c(length(a),length(b),length(c),length(d),length(e),length(f))
次に繰り返し関数rep()を用いて、各種事件条件の割り振り。
＞　condition <-factor(rep(c("control","treatment1","treatment2","treatment3","treatment4","treatment5"),times))


そしてANOVA解析

＞Data <- data.frame(z,condition)
でデータシートを作って
>head(Data)
で確認

z condition
1 23.246773 control
2 3.831205 control
3 27.163325 control
4 26.225653 control
5 22.088505 control
6 21.728285 control

こんな感じ。。

＞fm1 <- aov(z~condition, data=Data)
＞anova(fm1)

Analysis of Variance Table

Response: z
Df Sum Sq Mean Sq F value Pr(>F)
condition 5 20578.7 4115.7 217.81 < 2.2e-16 ***
Residuals 474 8956.9 18.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

有意差ありとでる。

そこでサブ解析も行うと、

＞pairwise.t.test(z, condition, p.adjust="bonferroni")

Pairwise comparisons using t tests with pooled SD

data: z and condition

control treatment1 treatment2 treatment3 treatment4
treatment1 < 2e-16 - - - -
treatment2 < 2e-16 6.4e-10 - - -
treatment3 1 < 2e-16 < 2e-16 - -
treatment4 < 2e-16 1 1.2e-07 < 2e-16 -
treatment5 < 2e-16 1.3e-08 1 < 2e-16 1.7e-06

P value adjustment method: bonferroni

実験条件のどれとどれの間に差があるのかわかるようになる。なお使用する統計量は何も書かなければ、適当にソフトが選んでくれるようだ。

box plotとdot plotの重ね合わせ(How to overlay your box plot with dot plot by R?)

次にやりたいのは、boxplotとdotplotの重ね合わせでしょうか？

#乱数で、サンプル数100, 平均20, 標準偏差5のデーターをつくる
> a<-rnorm(100,20,5)
#乱数で、サンプル数60, 平均10, 標準偏差3のデーターをつくる
> b<-rnorm(60,10,3)
#乱数で、サンプル数80 平均5, 標準偏差3のデーターをつくる
> c<-rnorm(80,5,3)

ここで問題となるのは、dotplotを生成する関数stripchart()が複数のオブジェクトを入れられないことです。
そこでdata.frame()をつかってデーターフレームを作りたいところですが、このケースのようにサンプル数が違う場合data.frameが働きません。そこでb, cのあまりの部分に空の値である。NAを放り込んで、サンプル数を同じにしてデーターフレームをつくってみます。

#b[61]～b[100]にNAを入れる
> b[61:100]<-NA
#c[81]～c[100]にNAを入れる
> c[81:100] <-NA
#ベクトルa,b,cからなるデーターフレームを作る
> data<- data.frame(a,b,c)

#確認
> data

ここまでくればだいたいOKでしょう。




といったグラフが得られるでしょう。
このままだとやはりみにくいので、

> boxplot(data, ylim=c(0,35),outline=F,col=c("Blue", "Red", "Green"),names=c("control","Treatment1", "Treatment2"),plot=T,ylab="% of CD34+ cells")

> stripchart(data, vertical = TRUE, pch = 21, col = "black", bg =NA,method = "jitter", add = TRUE)
として


とか

> boxplot(data, ylim=c(0,35),outline=F,col=c("Blue", "Red", "Green"),names=c("control","Treatment1", "Treatment2"),plot=T,ylab="% of CD34+ cells")

> stripchart(data, vertical = TRUE, pch = 21, col = "maroon", bg ="orange",method = "jitter", add = TRUE)
として



とかグラフをつくってみるとよいのではないでしょうか？

*Rスクリプトをうまくブログに乗せる方法がわからず、以前の表記に戻しました。