1. カウントデータを整形して R に読み込ませるこ とができる
2. R を使って簡単なグラフを書くことができる
3. R のパッケージを使って発現変動遺伝子抽出が できる
セクション 1
発現変動遺伝子抽出
ターミナル上でUNIXを使っておこなってきたのはテキスト データの解析である。次世代シーケンサー使うことの大きな 利点として、定量的であることが挙げられる。これを生かす ためには、数の解析が必要である。ここではRを使って数の 解析を行う。Rには非常に多くの統計パッケージがそろって いて便利である。
Rのホームページは http://www.r-project.org/ である。イ ンストール方法については、このサイト
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/01.html で非常 に詳しく日本語で書かれている。筑波大学のミラーサイト http://cran.md.tsukuba.ac.jp/bin/macosx/ よりダウンロード する2012/03/24現在での最新バージョンは2.14.2であるよう だ。私の環境では2.13.0を使っている。
Rを用いた次世代シーケンスデータ解析については、門田 先生のホームページが丁寧でわかりやすい。http://
www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html かなり勉強させ て頂いたし、パッケージのインストールの不具合については 直接教えてもいただいた。ここではデータの整形を中心に一 通りの発現変動遺伝子抽出をデモンストレーションする。今 回でもデータとするのはSRAに登録されているマウスの肝臓 RNA-seq データ
http://trace.ddbj.nig.ac.jp/DRASearch/run?acc=SRR001358 と、同じくマウスの骨格筋データ
http://trace.ddbj.nig.ac.jp/DRASearch/run?acc=SRR001361 である。SRR001358, SRR001361 のそれぞれのアクセッショ ン番号を検索していただくとすぐに出てくる。これをbowtie を用いてcDNA(
ftp://ftp.ensembl.org/pub/release-66/fasta/mus_musculus/cdna/ 内の
Mus_musculus.NCBIM37.66.cdna.all.fa.gz) にマッピングし、-S オプ ションで.sam 形式の出力を得てそれを awk -F "\t" '{print
$3}' | sort | uniq -c の一連のコマンドラインでカウントデー タに変換した。作業は/Users/nor/ ディレクトリ以下に/ Users/nor/otherPJ/macde/ ディレクトリを作って行った。
このようになっている。骨格筋側もあるので、
これら２つのファイルを合体、整形してRの読み込めるファイ ルを作ろう。今のファイルには、どちらか一方でしか発現し ていない遺伝子の記述に問題がある。ここで使うのはjoin と いうコマンドである。早速作業を始める。
まずはどちらか一方でも発現した遺伝子についてのリストを つくる。
nor$ cat gene_counts_liver.txt gene_counts_skemuscle.txt | awk '{print $2}' | sort | uniq > genename.txt
たくさんのデータがある場合には、リファレンスファイルか らこのリストを作ってしまった方が楽だ。その場合には
$ more Mus_musculus.NCBIM37.66.cdna.all.fa | grep ">" | awk -F " " '{print $1}' | tr -d ">" | sort > genename.txt
とすればよい。次に、カウントデータファイルを整形してリ ストファイルとのjoin に備える。
$ more gene_counts_liver.txt | awk ‘{print $2”\t”$1}’ > temp1
カラムの順序を入れ替えただけである。join は一列目を見る ので、このように整形しておくことが重要である。また、
sort されていることも重要で、今回の場合には発現量を推定
する操作で遺伝子名がsort された状態になっているのでこの まま使える。tophat やCASAVA など他の方法で得た発現量 データを扱う際にはこの点をよく注意して欲しい。もしsort が必要な場合には more filename | sort > filename2 とすれば よい。引き続き骨格筋のファイルを整形する。
nor$ more gene_counts_skemuscle.txt | awk '{print
$2"\t"$1}' > temp2
整形したそれぞれの発現量ファイルを遺伝子リストとjoin さ せる。
コマンドは
nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp1 | awk '{print
$1"\t"$2+0}' | more
これをtemp11に書き出して
nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp1 | awk '{print
$1"\t"$2+0}' > temp11
とする。同様に骨格筋の方も
nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp2 | awk '{print
$1"\t"$2+0}' > temp21
ここで、wc コマンドを使って行数が一致していることを確か めよう。
temp1 とtemp2 は行数が違うのにtemp11 とtemp21の行数が 一致していることがわかる。UNIX上のテキスト操作では、狙 い通りにものが進んでいないことがあるので、こうやって処
理の度に処理の完了を確認していくことが重要である。最後 にこれらのファイルを一つにまとめてヘッダーをつけよう。
まとめるのにはpaste コマンドを用いる
左側（肝臓）と右側（骨格筋）の遺伝子名が一致しているこ とを確かめてほしい。片方でだけ発現している遺伝子の項目 には、もう片方で０が入っているのがポイントだ。
３列目はいらないのでawk で削除。これをtemp3として書き 出す。
nor$ paste temp11 temp21 | awk '{print $1"\t"$2"\t"$4}' >
temp3
次にヘッダーを作ります。
$ vi header.txt
とするとvimというテキストエディタが起動するのでi を押し てINSERTモードにして
genenameタブliverタブskmus、と書き込んで[esc]ボタンを押
し、[:][w][q]と入力しエンター。vim の使い方はここに詳し
い。http://archiva.jp/web/tool/vim_basic.html
もちろんテキストエディタを開いて同様の入力をしたファイ
ルをheader.txtという名前で保存しても同じことであるが、改
行コードがCRになってしまったりするので、mi
http://www.mimikaki.net/ を使って改行コードをLFに直す必
要が出てくるかもしれない。
では、これを先のtemp3と合体させる。
これが正解。これに好きな前をつけて書き出して完成だ。
nor$ cat header.txt temp3 > for_R_liv_vs_skm.txt
次にこれをRに読み込ませよう。Rを起動すると、 こういった感じである。データを読み込ませるにはデータま でのパスを書けばよいので、
data<-read.table("/Users/nor/otherPJ/macde/for_R_liv_vs_
skm.txt",header=T,row.names=1)
とする。簡単に説明すると、”data”として（名前は好きなも ので何でもよい）/Users/nor/otherPJ/macde/ ディレクトリ
に入ってるfor_R_liv_vs_skm.txt を読み込め、ヘッダーがつ いていて、それぞれのデータの名前は１列目に書いてある。
という意味がある。続いて、data をattachしてliver とskelm を読み込ませる。summary を使うとデータの確認ができるの で習慣としてここまでを一気に行うとよい。
試しにプロットしてみよう。
plot(log(liver),log(skelm),pch=20)
この図では一つの遺伝子をひとつの点として、liver での発現 量を水平軸に、skelmでの発現量を垂直軸にプロットしてい る。左上にあるのと、右下にあるのがなんとなくの発現変動 遺伝子(Differentially Expressed Genes/ DEG) である。これに
DESeqという解析パッケージを用いてDEG抽出を行う。
DESeqのインストールには門田先生のHP
http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#install の 通りにやればうまくいく。
source("http://www.bioconductor.org/biocLite.R") biocLite("DESeq", dependencies=TRUE)
としてDESeqを手に入れる。DESeqを選んだ理由は、これが
反復のないデータ(n=1)でも解析できるからである。以下にコ マンドラインを示す。
out_f <- "DESeq_out_liv_vs_skelm.txt"
param1 <- 1 param2 <- 1 library(DESeq)
data.cl <- c(rep(1, param1), rep(2, param2)) cds <- newCountDataSet(data, data.cl)
cds <- estimateSizeFactors(cds) sizeFactors(cds)
cds <- estimateVarianceFunctions(cds, method="blind") out <- nbinomTest(cds, 1, 2)
logratio <- out$log2FoldChange FDR <- out$padj
rank_FDR <- rank(FDR, ties.method="min") tmp <- cbind(rownames(data), data, logratio, FDR,
rank_FDR)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
以上である。
結果は、/Users/nor/ディレクトリ以下に
DESeq_out_liv_vs_skelm.txtという名前で格納されている。
これをエクセルなどに貼り、フィルタを使ってFDR順に並び 替えるとDEGが抽出できる。FDRはFalse Discovery Rate の 略で、ここではp値と同じような判断の指標として扱ってい る。
今回 の比較ではFDRが十分に小さいものが得られなかった。私の 手元のデータ（同じ方法でDEG55個(FDR<0.05)）のプロッ トと比較すると発現量の小さい遺伝子がかなり多いように見 える。ここまででDEG抽出は終了である。
続いて、もっとおおまかにトランスクリプトームをつかむ方 法を紹介する。
par(mfrow=c(1,2))
barplot(matrix(sort(liver/sum(liver))))
barplot(matrix(sort(skelm/sum(skelm))),ylim=c(0,1)) とすると以下の図が描ける。
左が肝臓で、右が骨格筋のものである。それぞれの遺伝子の 発現量（リード数）を多い順に並び替えてプロットした。も うひとつ、ハエの脂肪体のデータを載せる。こちらの図は先 のふたつよりもかなりかたよっているようだ。