まずは、以下の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と入力してください。 プロットを見てみると、、、
.png)
むう。なかなか面白いですね。一昔前なら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
}
}
2020