ベンチ屋のベンチ屋のためのベンチ屋によるRの解説〜ENCODEデータで遊んでみる

2017年09月19日（火）

ベンチ屋のベンチ屋のためのベンチ屋によるRの解説〜ENCODEデータで遊んでみる

投稿者: 中川真一

　恥ずかしながら、僕自身、ウェットの研究者がピンセットとピペットマンだけ握っていれば良い時代はとうの昔に終わったことに気づいたのは、実はつい最近。今更データサイエンティストを目指すのもなんだしなあ、、、と半ば諦めてたところはあるのですが、ちょろっと触ってみるとこれが意外と楽しいんですね。HelloとShitとGreatだけの英会話みたいなレベルでも、意思疎通できるというのはそれだけで喜びを与えてくれるものです。それに加えてI love youだけ覚えれば人生の90％のミュニケーションは事足ります。というわけで、I love you R！

まずは実験のメタデータが入ったファイル"metadata.tsv"をダウンロードしてきます。これはENCODEのshRNA-RNAseqのbatch downloadをするためのファイルの一行目にかかれているURLから取ってこれます。ファイルのダウンロードはR上でやっても良いのでしょうが、普通にブラウザのリンクから取ってくるのが手っ取り早いですね。これを、デスクトップ上に作ったENCODEフォルダにでも放り込んでおきます。次にRを立ち上げて、作業ディレクトリをそのフォルダに設定し、さらにメタデータをmetaというデータフレームに取り込んでおきます。また、その後の作業で使うライブラリggplot2とdata.tableも読み込んでおきます。

library(data.table)
library(ggplot2)
setwd("~/Desktop/ENCODE")
meta<-read.delim("metadata.tsv")

初回、まだパッケージがインストールされていない時はR studioのtools>install packagesからインストールが必要です。

　次に、このメタデータから、各遺伝子にマップされたリードをカウントしたデータだけ選択的にダウンロードするためのファイルを作ります。hg19とGRCh38にマップされたものの両方が用意されているようですが、ここでは最新版のGRC38の情報だけ取ってきます。データの種類はmetadata.tsvの中の”Output.type”という列に格納されているのでそこが"gene quantifications"となっているものを選び、さらに、ゲノムは"Assembly"という列に格納されているので、GRCh38のものを選べば良いだけですね。これをmeta1というデータフレームに入れておきます。ついついエクセルに手が伸びてしまうところですが、ここはじっと我慢してRでやってしまいます。クリック操作は直感的で便利ですが、どこをクリックしたかという記録が残りません。変なミスをしてしまった時に後から検証するのが大変。その点、Rなら安心です。何をやったのか全て記録が残るのがR、クリック操作の記録が残らないのがExel。実験というのは記録を残すのが基本ですから、Rと実験屋の相性は実はExelより良いのでは？と思う今日この頃です。

meta1=meta[meta$Output.type=="gene quantifications"&meta$Assembly=="GRCh38",]

　ちなみにですが、Rとかだと「等しい」は=一個でなくて==なんですね。一個の=だと代入になってしまうらしい。ちょっと前までこの辺りで1時間以上浪費していたのはヒ・ミ・ツ。むだ話はともかく、必要なカウントデータのダウンロードリンクは"File.download.URL"という列に格納されているので、それを取り出してきて、先頭にターミナルのダウンロードコマンドである"wget"をつけたdownload.commandという変数を作り、それをdownload.txtという名前のテキストファイルとして出力します。この時に何もオプションをつけないと各要素にクオーテーションマークが付いてしまうのと余分な行名が付いてしまうので、それを"quote="オプションと"row.names="オプションをつけて外しておきます。

download.command=paste("wget",meta1$File.download.URL)
write.table(download.command,file="download.txt",quote=FALSE,row.names=FALSE)

　このファイルをテキストエディットか何かで開き、一行目のXを取り除き、代わりにターミナルで実行するための一行目、おまじないの#!/bin/shを入れて、セーブしておきます。そうすると、最初の数行はこんな感じ。

#!/bin/sh
wget https://www.encodeproject.org/files/ENCFF663TDK/@@download/ENCFF663TDK.tsv
wget https://www.encodeproject.org/files/ENCFF168GNX/@@download/ENCFF168GNX.tsv
wget https://www.encodeproject.org/files/ENCFF209EGI/@@download/ENCFF209EGI.tsv

　次にRを一旦離れ、ターミナルを立ち上げてこのファイルを実行します。ファイルがいろいろできてごちゃごちゃしてくるので、ダウンロードしたファイルはENCODEのフォルダの下にcountというフォルダを作って入れておくことにします。そのフォルダに移動して、先程のdownloadファイルをshコマンドで実行！

cd /Users/Nakagawa_MacBook/Desktop/ENCODE
mkdir count
cd count
sh ../download.txt

　そうすると画面が次々と流れていって、いつもながらの何もしていないのに何かすごく仕事をしている感が得られる至福の時間がやってまいります。この画面さえ出しておけば、ボーッとしていても仕事をしているオーラを出すことができますね。ネット環境が太いところではちょっと長めの昼休みから帰ってくる頃には、自宅のネット環境ぐらいショボくても一晩経てばダウンロードが終わってずらりファイルがダウンロードされてきます。ちなみにRでも外部ファイルをダウンロードするdownload.fileというコマンドがあるようですが、そちらではエラーが出てうまくダウンロードできないようです。

　さて、これで、適当なファイルを開けば自分の興味のある遺伝子がどれぐらい発現しているか数値ですぐわかるわけですが、各RNA結合蛋白質のノックダウンで目的の遺伝子のFPKMがどれだけ変わったかがわかる一覧の表が欲しいところです。というわけで、それぞれのファイルからFPKMの値だけ取り出して、ENSEMBLEのgene IDごとに並べた表を作ってみます。

　エクセルだったらそれぞれのファイルを開いてFPKMの列をコピー、それをひたすらペーストしていくということになるんでしょうが、Rなら簡単！のはずが、、、実は結構はまってしまいました。非情なエラーメッセージが出るたびにそれをコピペしてググって、、、と苦闘していたら、ラボの小松くんがサクッと、こんなんでよければできましたよ、と。おおーっ！持つべきものは良き同僚、良き院生です。まずはlist.filesコマンドでダウンロートしてきたファイルリストの一覧を撮ってきてflという変数に突っ込みます。ここではcountフォルダの中には目的のファイルしか入っていませんが、他のファイルも多数ある時は、patternオプションを指定してやれば特定の拡張子や名前を持つものだけを取ってくることもできるようです。そして、読み込むべきファイルが何個あるかをlengthコマンドで撮ってきてnumという変数に入れておきます。なんとなく仕事をやっている感を出すために、読み込むべきファイルが何個あるかも出力してみます。

setwd("~/Desktop/ENCODE/count")
fl=list.files(pattern="*.tsv")
num=length(fl)
print(paste("number of files to be imorted =",num))

　次に、データを結合していくためのテンプレートを作ります。まずは、fnという変数にダウンロードするファイル名、この場合はflの1番目に入っているファイル名を突っ込み、爆速のファイル読み込みコマンドであるfreadでdataというデータフレームに取り込んで、gene_idという列とFPKMという列だけ抜き出します。このままだと列名がFPKMのままなので、ファイル名から.tsvを除いた名前に置き換えておきます。ちゃんと仕事をしている感を出すために、どのファイルを読み込んだかもログで出力してみます。

fn=fl[1]
data=fread(fn)
data=data[,c("gene_id","FPKM")]
colnames(data)[2]=sub(".txt","",fn)
print(paste("read",1,"/",num,fl[1]))

　あとはforループで次々とファイルを読み込んでは上で用意したテンプレートに結合していきます。ここでは、新たに良いこんだファイルはtempというデーターフレームに一時的に入れておいて、tempとdataの共通の列、すなわちgene_IDを用いてmergeコマンドで結合しています。この手の作業をする時にmergeコマンドを使うのが地味に大事なところで、cbindコマンドでも同じようなことはできますが、gene_IDが全てのファイルで同じ順番に並んでいるとは、、、限らない。欠損値がないとは、、、限らない。多分そんなことないんでしょうが、不測の事態も予測してよりrobustな系で作業を進める。これ、実験の鉄則ですね。ものすごく仕事が進んでいる感を出すために、読み込むべきファイルのうち何番目の何を読み込んだかをログで出力してみます。おーっ！進んでる進んでる。コンピューターと会話している感がすごい！

for (i in 2:num){
  fn=fl[i]
  temp=fread(fn)
  temp=temp[,c("gene_id","FPKM")]
  colnames(temp)[2]=sub(".tsv","",fn)
  data=merge(data,temp,by="gene_id")
  print(paste("read",i,"/",num,fl[i]))
}

　完成したdataファイルの中身をheadコマンドでちょっと見てみます。うんうん。これだ、欲しかったファイルは。ありがとう小松くん。行数も列数もバッチリ。一列目はgene_IDになってますので、列数は読み込んだファイル+1にちゃんとなっています。この作業は最初に一回やっておけば良いので、玉の子のデーターフレームのファイルをtsv形式で出力しておきます。

head(data)
setwd("~/Desktop/ENCODE")
write.table(data, file="shRNA_RNAseq_FPKM.txt",quote=FALSE, sep='\t',col.names = TRUE,row.names = FALSE)

　次からはこれを読み込んでやれば、手っ取り早いです。ここから先は、保存したデータをdata1に読み込んだもので進めてみます。データーフレームは、setkeyコマンドで参照したい列を指定しておくと、簡単に特定の列の情報を引っ張ってくることができます。ここでは試しにACTBで行ってみましょうか。ACTBのGENCODEv24のgene_IDはENSG00000075624.13です。

data1<-fread("shRNA_RNAseq_FPKM.txt")
setkey(data1,"gene_id")
data["ENSG00000075624.13"]

　実行させてみると、途中で切ってますが、目的の数値は確かに得られています。とはいえ、各実験のIDが何をKDしたのか情報がないので、何が何やらサッパリです。

              gene_id ENCFF003AYH ENCFF003BXW ENCFF004CAV ENCFF006XHX ENCFF007FIZ ENCFF008IBD
1: ENSG00000075624.13      725.49      730.23      931.73     1491.28      540.56      919.24
   ENCFF009QBJ ENCFF010VKI ENCFF012QOZ ENCFF013LUC ENCFF016MOB ENCFF016WUI ENCFF017FCR
1:      502.93      517.87      452.93      781.17      616.85      512.14       902.8

　というわけで、このデータを初っ端にダウンロードしたmetadataとmergeして、その他の情報が参照できるテーブルを作ってみます。この場合、各列にいろいろな情報、各行にそれぞれのKD実験が入る方がわかりやすいので、まずはdata1から目的の遺伝子(gi: gene of interest)を引っ張り出したデータテーブルFPKM_giを作り、行列をtコマンドで入れ替え、一行目には先程の列名が入ってしまっているのでそれをのぞいてやって、File.accessionとrnをキーにしてmergeしてます（rnはdata.tableが自動的につけた列名）。全ての情報が入ったデータテーブルはFPKM_mergeで、発現量はV1の列に、ノックダウンした遺伝子名はExperiment.targetの列に、細胞名はBiosample.term.nameの列に入っています。最後に、この後グラフを書くときにV1の列は数値データでなくてはならないので、文字列から数値にas.numericコマンドを用いて変換してます。

setkey(data1,"gene_id")
gi="ENSG00000075624.13" # set the ENCODEv24_ID for the gene of interest
FPKM_gi<-data1[gi] # select FPKM data for the gene of interest
FPKM_gi<-as.data.table(t(FPKM_gi),keep.rownames=TRUE) #exchange the row and column
FPKM_gi<-FPKM[-1,] #delete row containing extra information
FPKM_meta=merge(meta,FPKM_gi,by.x="File.accession",by.y="rn")
FPKM_meta$V1<-as.numeric(FPKM_meta$V1)

　さて、これで図を描くためのデータは揃いました。HepG2とK562の二つの細胞を使ったデータが入っていますから、それぞれ分けておきます。

FPKM_HepG2=FPKM_meta[FPKM_meta$Biosample.term.name=="HepG2",]
FPKM_K562=FPKM_meta[FPKM_meta$Biosample.term.name=="K562",]

あとはggplot2で書くだけです。ggplot2の便利なところは、見た目の綺麗さもありますが、なんでも良いから情報を色々突っ込んだ元データさえあれば、それをカテゴリー別に分類しなくても一発でグラフを作ってくれるところにあるような気がします。この例でも、各RNA結合蛋白質に対して各々4つの実験があるのですが、それらの平均値をとるセルを作って、後でまとめて、、、なんていうexel風の作業は全くしなくても良い！しかも、何をしたのかが非常にわかりやすく記録に残るので、人に教えるときも、人から教わる時も、Google先生に教わる時も、何をすれば良いかが明確なんですね。なんで今まで使ってこなかったんだろうと反省しきりですが、ここではboxplotを使ってみます。X軸とY軸に何を使うのかをaesで指定し、geomでどういうグラフを書きたいか指定します。あとは、X軸の名前を90度回転させたいなあ、と思ったら"R　X軸　90　回転"で検索すればQiitaのエントリーが教えてくれますし、Y軸の値の順にX軸の項目をソートした図を作りたかったら"R ソート　グラフ　X軸"でこれまた別のQiitaのエントリーが教えてくれます。Google先生にQiita先輩、みたいな感じでしょうか。バイオインフォの専門家のように高度な解析をするソースコードを書くレベルに到達するのには高いハードルがあるかもしれませんが、F1レーサーでなくても愛車のアルトくんには乗るわけですし、何よりもこんな簡単な操作でいろいろな実験がショートカットできるのであれば、これに乗らない手はない、です。

#plot data HepG2
ggplot(FPKM_HepG2,aes(x=reorder(x=Experiment.target,X=V1),y=V1)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  ylab(paste("FPKM",gi)) +
  xlab("KD gene")
ggsave(filename=paste("plot_HepG2_",gi,".png",sep=""),height=5,width=nrow(FPKM_HepG2)/15,limitsize=FALSE)
  
#plot data K562
ggplot(FPKM_K562,aes(x=reorder(x=Experiment.target,X=V1),y=V1)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  ylab(paste("FPKM",gi)) +
  xlab("KD gene")
ggsave(filename=paste("plot_K562_",gi,".png",sep=""),height=5,width=nrow(FPKM_K562)/15,limitsize=FALSE)