ex301: Count data import and scatter plot

''arab2.txt''は、6 libraries (2 groups x 3 biological replicates) のシロイヌナズナRNA-seqのデータである。すでにマッピング済みで遺伝子毎のリードカウントがタブ区切りテキストとして提供されている。このexerciseでは、テーブルの中身を確認しデータの概要を把握する基本テクニックを習得する。

Data

(~/gitc/data/SS/以下にある。以下R等コマンド例ではファイルのパスは省略しているので、適宜適切なパスを指定すること。）

• arab2.txt : count table

シロイヌナズナに細菌 Pseudomonasを感染させてどのように遺伝子発現が変動するかを調べた実験。6 libraries (2 groups x 3 biological replicates) のデータ。すでにマッピング済みで遺伝子毎のリードカウント（ノーマライズはされていない）がタブ区切りテキストとして提供されている。m1, m2, m3 が mock (control), h1, h2, h3 が　Pseudomonas感染。

以下全てR環境。（受講生の皆さんは、ローカル環境でRを実行してください）。

Inspect table with R

R でテーブルの内容を確認し、さらにscatter plotを書いてみよう。

Data import

> dat <- read.delim("arab2.txt", row.names=1, head=T)  # read tab-delimited text
> head(dat)
          m1 m2 m3 h1 h2 h3
AT1G01010 35 77 40 46 64 60
AT1G01020 43 45 32 43 39 49
AT1G01030 16 24 26 27 35 20
AT1G01040 72 43 64 66 25 90
AT1G01050 49 78 90 67 45 60
AT1G01060  0 15  2  0 21  8

Inspect table

> dim(dat)
[1] 26221     6

Q：dimコマンドは何をするものですか？また、その結果得られた、26221 と 6　は何を意味しますか？

Hint: コマンドの詳細を調べたいときは、Rコマンドラインで、?に続けてコマンドを打てば、ヘルプ画面が表示されます。

> ?dim

Inspect data by column

それぞれのライブラリの、リードカウント合計を計算してみよう。

例としてm1カラムの合計を計算する。

> sum(dat$m1)  
[1] 1902032

他のカラムも合計を計算しよう。

カラムを一つずつコマンドを叩いて計算するのは面倒である。colSums を使ってみよう。

> colSums(dat)
     m1      m2      m3      h1      h2      h3 
1902032 1934029 3259705 2129854 1295304 3526579 

Q4: 約2万5千遺伝子にはカウント０のものから非常にたくさんのカウントをもつものがある。カウントの、i)最大値、最小値、平均値、中央値を計算しなさい。ii) ヒストグラムを書きなさい。

m1 を例に実行例を示す。

i) 最大値、最小値、平均値、中央値

> sum(dat$m1)
[1] 1902032
> max(dat$m1)
[1] 61791
> min(dat$m1)
[1] 0
> mean(dat$m1)
[1] 72.5385
> median(dat$m1)
[1] 9

ii) ヒストグラム

> hist(dat$m1)

しかし、このグラフではあまり特徴がつかめないと思う。対数をとってみよう。

> hist(log10(dat$m1 + 1), breaks="Scott")

RNA-seqの解析では、カウントやノーマライズドデータの対数をとって処理することが多いので覚えておこう。

Scatter plot

Q: m1 vs m2 をscatter plotで比較しよう。

> plot(dat$m1 + 1, dat$m2 + 1, log="xy")

それぞれ+1しているのは、log0は計算できないため。+1して下駄を履かせている。

Q: ほかのライブラリどうしもscatter plotを描いて比較しよう。

Play with scatter plot

plotなどのグラフィックス関数には、様々な引数を与えることによって非常に多くの描画パラメータを変更でき、グラフの見栄えを変更することが出来る。scatter plotの色、形、などを変更する練習をしてみよう。

以下のコマンドをテンプレートにして、col（色）, pch（点の形状）, cex（点の大きさ）, main（グラフのタイトル）, xlab|ylab （x軸やy軸のラベル）を変更して、変化を確認しよう。

> plot(log2(dat$m1 + 1), log2(dat$h1 + 1), col="DarkBlue", pch=16, cex=0.6, 
       main="scatter plot", xlab="m1", ylab="h1")

発展：透明度の指定

RNA-seqのカウントデータの散布図はデータ数が多いためにデータが集中しているところは真っ黒に塗りつぶされてしまい、データの分布がわかりにくい。透明度を指定するとデータの多い部分は濃く、少ない部分は薄く表示されるため、データの分布を視覚的により把握しやすい。

> plot(log2(dat$m1 + 1), log2(dat$h1 + 1), col=rgb(0,0,1, alpha=0.2), pch=16, cex=0.6, main="m1 vs h1", xlab="m1", ylab="h1")

色指定オプション col= でalpha を指定しているのがポイント。

（参考）Inspect table with spreadsheet software

Rに不慣れな場合、表計算ソフトウェア（MS Excel, OpenOfficeなど）でも簡単なデータ確認は可能である。

MS Excelで、m1とm2のscatter plot（散布図）を書いてみよう。このふたつは同一コンディションのbiological replicateなので、発現パターンは両者で良く似ているはずである。次にm1とh1を比較しよう。このふたつはコントロールと実験群の比較なので有る程度の発現パターンの違い（すくなくともm1 vs m2よりも大きい違い）が期待される。

ヒント：xy軸ともに対数をとること。

Tips: MS Excelの「条件付き書式」機能を使い、カウントデータの寡多をセルの色に反映させるように設定させると、ヒートマップに似た表示を実現できる。