2017年09月18日(月)

ENCODE shRNA-RNAseqの遺伝子別全データを可視化するR script

投稿者:
興味のある遺伝子の転写産物を制御するRNA結合蛋白質がどれか一目でわかる! 興味のある遺伝子の転写産物を制御するRNA結合蛋白質がどれか一目でわかる!

もう10年ほど前になりますが、Gomafuの核内局在因子やXistの染色体極在化因子を探そうとして僕らがとったアプローチは、RNA結合蛋白質に対してデザインしたカスタムsiRNAライブラリーを用いて片っ端からノックダウンしていって、その効果を調べて真打を探す、というものでした。HGWさんのたぐいまれなる強運もあり、このアプローチは大成功を収めたわけですが、ENCODEのshRNA-RNAseqのデータが使える今日だったら、実験を始める前からデータは手許にある!ということで、かなりの部分をショートカットすることができるわけです。というわけで、R初心者によるR講座第ウン段。今回は、自分の興味のある遺伝子の転写産物が、特定のRNA結合蛋白質をノックダウンした時にどのような変化をするのか、ということを可視化する方法を紹介します。まあ、そのうちUCSC Genome Browserでもすぐに見れるようになるんでしょうが、ちょっとだけ先を行こうということで、、、

 

まずは、以下の4つのファイルをダウンロードします。zipファイル以外は下の方の「添付ファイル」と書いてある全くクリックしてもらおうという気力の感じられないリンクからもダウンロードできます(ただしなぜか下のリンクからダウンロードしたplot.Rには.txtの拡張子がつくみたいなので外してください)。

shRNA_RNAseq_FPKM.txt.zip
shRNA_RNAseq_meta.txt
Xref.txt
plot.R

zipファイルは解凍し、全てデスクトップに作成したENCODEフォルダに入れておきます。これは色々なRNA結合蛋白質をノックダウンした時に遺伝子発現量がどのように変化するのか、ENCODEのダウンロードページからデータを落としてきてまとめた表になります。この表を作った手順については次のエントリーに書いてますのでそちらを参照してください。また、ggplot2とdata.tableがRにインストールされていない場合は、パッケージマネージャーで事前にインストールしておいてください。 次に、興味のある遺伝子のGENCODEv24のIDを調べます。例えば、そうですね。昔なつかしの網膜幹細胞因子、WNT2bを見てみましょうか。身近なgene symbolはわかっているけどGENCODEv24のIDがわからない時は、Rを立ち上げて以下のようなコマンドで調べることができます。例えばWNT2Bでしたら、

library(data.table)
Xref=fread("Xref.txt")
Xref[Xref$HGNC_symbol=="WNT2B",]

結果は

> Xref[Xref$HGNC_symbol=="WNT2B",]
   HGNC symbol Gene stable ID version
1:       WNT2B     ENSG00000134245.17

ですので、ENSG00000134245.17であるということが分かります。まあ、Xrefのファイルはそんなに大きくないので、こんな七面倒くさいことしないでも普通にテキストエディットで開いて⌘Fで検索した方が早いかも、、、 IDが準備できたら、

setwd("~/Desktop/ENCODE")
source("plot.R")

と打ち込むだけ!。すると、GENCODEv24のIDを聞いてくるので、それを入力すると、ENCODEフォルダに"plot_HepG2_ID名.png"と"plot_K562_ID名.png"いうファイルが出力されてきます。プロットが終わるとまたIDを聞いてくるので、連続して図を作成することも可能です。図の作成が終わると次のIDを聞いてきますので、終了するときはnと入力してください。 プロットを見てみると、、、

むう。なかなか面白いですね。一昔前なら10年分ぐらいの仕事の結果が手に入った感じでしょうか。これを見るとWNT2BはK562ではほとんど発現しおらず、HepG2でも発現していてもかなり弱いようですが、RNA結合蛋白質のKDでいろいろ発現変動することが分かります。ふうむ。SRFBP1のKDで上がってくるのか。ということは、、、 頭の中でいろいろな妄想が浮かんできますね。どう活用するかはアイデア次第!!このデータが10年前にあったらもう少しWNT2Bの仕事を続けてた、、、かも。

ちなみにコードはこんな感じです。次回のエントリーでもうちと詳しく各ファイルの作り方のメモを書きます。

#setup
library(data.table)
library(ggplot2)
setwd("~/Desktop/ENCODE")
data=fread("shRNA_RNAseq_FPKM.txt")
setkey(data,"gene_id")
meta=read.delim("shRNA_RNAseq_meta.txt")

while (TRUE){
  
 #input gene of interest ID
 gi=readline("enter GENCODEv24 ID: ") # set the ENCODEv24_ID for the gene of interest
 
 #create data table for ggplot2
  FPKM_gi<-data[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_gi[-1,] #delete row containing extra information
  FPKM_meta=merge(meta,FPKM_gi,by.x="File.accession",by.y="rn") #assign metadata
  FPKM_meta$V1<-as.numeric(FPKM_meta$V1) #change data type for ggplot2
  
  FPKM_HepG2=FPKM_meta[FPKM_meta$Biosample.term.name=="HepG2",]
  FPKM_K562=FPKM_meta[FPKM_meta$Biosample.term.name=="K562",]
  
  #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)

  #continue?
  x=readline("continue? (y/n): ")  
  if (x=="n") {
    break
  }
  }

test

中川 真一

北海道大学 薬学研究院 教授
▶ プロフィールはこちら

コメントする

コメントはどなたでも歓迎します。メンバーの方はログインしてからコメントして下さい。

ブログアーカイブ

ログイン

サイト内検索