eCLIPはご存知の通り、多数のRNA結合蛋白質がどのRNAに結合しているのかを調べたデータベースで、ENCODEのホームページから全てデータダウンロードできます。ここではまず、UPA-Seq論文のデータで使用したテーブル、eCLIP_peak_counts.txtを使ってみます。各columnにはそれぞれのRNA結合タンパク質、各rowは各遺伝子のEnsembleのID、中に入っている数値は、それぞれの遺伝子にRNA結合蛋白質のピーク数がどれぐらいマップされていたのかのカウント数が入ってます。それぞれの遺伝子のgene symbolとEnsembleのIDの対応は、これもUPA-Seq論文で使ったv28_genelist.txtを使います。これらの元ファイルは、毎度おなじみ、この記事の一番下、クリックしてもらおうという覇気が全く感じられないリンクのところに置いときました。ダウンロードしてデスクトップのtSNEフォルダに入れときます。これらのテーブルをどうやって作るかは、、、ちょっとめんどくさいので次回のエントリーで。

まずは、eCLIP_peak_counts.txtとv28_genelist.txtをmergeしたデータテーブルm1を作っておきます。よくある作業ですね。Rのmergeコマンドでちょろちょろっと。

#tSNE analyses eCLIP
library(data.table)
setwd("~/Desktop/tSNE")

eCLIP_data <- fread("eCLIP_peak_counts.txt")
colnames(eCLIP_data)[1] <- "gene_id"
meta <- fread("v28_genelist.txt")
m1 <- merge(eCLIP_data,meta)

次に、いよいよtSNE解析です。R用のtSNEのパッケージは、Tools > Install packages...に行って、Repository(RCRAN)から簡単インストールできます。Rのバージョンによっては何やら文句を言われることがあるみたいですが、最新の状態にしておけば大丈夫です。ライブラリをインストールできたら、あとはデータテーブルの数値部分を、えいやっと投げるだけ。超簡単。こんな簡単でいいのかしら、、、みたいな。ここでは、mergeしたデーターテーブルm1の中から、遺伝子タイプ（80列目）と数値が入っている部分（2列目から79列目）だけを取り出したデータフレームdfを作り、さらにその行名にENCODEのIDを入れてます。また、このままだとちょっと僕のか弱いMacBookちゃんだとデータ多すぎるので、遺伝子にマップされている全ピーク数の合計が100未満の遺伝子には退場してもらってます。一行目は数値データではないので、一列目を除いたデータフレーム（df[,-1]）をmatrixに変換してからRtsneに投げてます。

#tSNE install
install.packages("Rtsne")
library(Rtsne)
#prepare data frame for tSNE analyses
df <- data.frame(m1[,c(80,2:79)])
row.names(df) <- m1$gene_id
df <- df[apply(df[2:79],1,sum) > 100,] #select abundant target genes
tSNE <- Rtsne(as.matrix(df[,-1]), check_duplicates = FALSE, verbose=TRUE)

で、コードを走らせると、

おおっ！なんかやってます。なんかやってくれてる感満載のログが次々と流れていく画面。タイムラプスとってる時とか、SIMの元画像集めてる時とかと似てますねこの感覚。全く働いてないのに、働いてる気になれる至福の時です。

ここではtSNEの結果をtSNEという変数に入れてますが、次元下げされたデータはtSNE$Yに格納されてます。X軸とY軸の数値はY$[,1]とY$[,2]で取り出すことができるますので、データフレームdfにtSNE_1とtSNE_2という列を追加して入れてます。ggplotの色分けができるように念のため遺伝子タイプをファクター型に変換し、あとはこれをおなじみggplotでpngファイルで書き出してやれば、、、

#plot tSNE data
library(ggplot2)
df$tSNE_1 <- tSNE$Y[,1]
df$tSNE_2 <- tSNE$Y[,2]
df$gene_type <- as.factor(df$gene_type)

g <- ggplot(df,aes(x=df$tSNE_1,y=df$tSNE_2,color = df$gene_type))+
  geom_point()
ggsave(filename = "eCLIP_tSNE.png",g,width = 15, height = 10)

お！なんか出てる出てる。tSNEの解釈にはいろいろ気をつけなければならないところがあるみたいですが、RNA結合蛋白質のパターンで分類すると、protein coding（発現量が多いものだけ見てるのでこれが最大のグループになってますね）がドカンとあるのに加えて、ノンコーディングRNAやsnoRNAなどがグループを作っている様子がよくわかります。しかもよく見るとprotein codingの中にもブロッコリーみたいに塊が見えてるような見えてないような、、、さらに個別の遺伝子をよくよく見てみると、なんだかNeat1とかMalat1とかFtxとかKcnq1ot1とかの機能性のlncRNAがクラスター作っているようにも見えます。本気でやれば、結講応用が利きそうですね。

NGS解析が広まるにつれて最近はこういったデータ解析ツールが本当に手軽に使えるようになりました。ちょうど組換えDNA技術の登場によって『發生学』やら『解剖学』をやっていた人がこぞって大腸菌を扱うようになったように、いわゆるウェットの実験屋がターミナルでかたかたコマンドを打つのもそう珍しい風景ではなくなってきました。ウェットとか、ドライとか、もう死語かもしれませんし、最近の学生さんなどはあまり意識もされてないのではないでしょうか。RNA研究の今後の展開が楽しみです！