MotDB


AJACS51/bono の変更点


[[講習会のページに戻る>AJACS51]]

&size(36){RおよびBioconductorを使ったバイオデータ解析};

講習では【実習】を皆さんとやっていきます。早くできてしまった人は、【発展】をやってみたりしてください。
Rのインストールは以下の統合TVを参考に行っておいてください。
- 参考:統合TV[[http://togotv.dbcls.jp/ja/wp-content/uploads/2013/03/togotv_banner.png>http://togotv.dbcls.jp/]]
-- [[統計解析ソフト「R」の使い方 導入編>http://togotv.dbcls.jp/20090313.html]]
-- [[統計解析ソフト「R」の使い方 導入編(MacOSX)>http://togotv.dbcls.jp/20090618.html]]
-- [[統計解析ソフト「R」の使い方 正規化編>http://togotv.dbcls.jp/20090319.html]]
-- [[統計解析ソフト「R」の使い方 ヒートマップ編>http://togotv.dbcls.jp/20091219.html]]
-- [[統計解析ソフト「R」での立廻り>http://togotv.dbcls.jp/20111107.html]]

----
~目次
#contents
----
* 0. はじめに [#d16282c7]
- 【実習0】Rをインストールしなさい。
#fold(←こたえは左の+マークをクリックすると出てきます,http://cran.ism.ac.jp/ からダウンロードしてインストールしてください)
- 【発展0】Rを快適に使うにはRstudioというソフトウェアがある。Rstudioをググッて見つけ、インストールしなさい。
#fold(←こたえは左の+マークをクリックすると出てきます,http://www.rstudio.com/ からダウンロードしてインストールしてください)
- 【発展1】今回使うRのコードやデータは[[githubのサイト>https://github.com/bonohu/AJACS51]]にあります。そこのデータをまとめてダウンロードするには、UNIXのコマンドラインで
 git clone https://github.com/bonohu/AJACS51
を実行しなさい。現在居るディレクトリ(current working directory)にどういった変化が起きたか?
#fold(←こたえは左の+マークをクリックすると出てきます,AJACS51というディレクトリが作成され、それ以下にgithubにあるファイルがすべてダウンロードされた)

* 1. パイチャート [#ib0f08d6]
-【実習1】Rを起動し以下のファイルに記述されたRのコマンド(''01piechart.r'')を実行しなさい。データファイルとして必要な''srabystudy.txt''をgithubのサイトから[[ダウンロード>https://raw.githubusercontent.com/bonohu/AJACS51/master/srabystudy.txt]]し、現在作業しているディレクトリにおいておく必要があります。
-【実習1】Rを起動し以下のファイルに記述されたRのコマンド(''01piechart.r'')を実行しなさい。データファイルとして必要な''srabystudy.txt''をgithubのサイトから[[ダウンロード>https://raw.githubusercontent.com/bonohu/AJACS51/master/srabystudy.txt]]し、現在作業しているディレクトリにおいておく必要があります。実行後、そのディレクトリに''pie1.png''というファイルが新たに生成されていることを確認しなさい。

 png("pie1.png")
 dat <- read.delim2("srabystudy.txt", header=F)
 names(dat) <- c("Study", "freq")
 dat <- cbind(dat, serial=seq(dim(dat)[1]))
 dat2 <- dat[ dat$freq != 0, ]
 pie(dat2$freq, labels=paste(dat2$Study, dat2$freq), col=rainbow(dim(dat[1])[dat2$serial]), main="SRA by Study")
 dev.off()

#fold(←左の+マークをクリックするとこのスクリプトの経緯が…,Bono H et al. "Systematic Expression Profiling of the Mouse Transcriptome Using RIKEN cDNA Microarrays" Genome Res 2003の[[図1>http://genome.cshlp.org/content/13/6b/1318/F1.expansion.html]]を作成するために使ったものです。実際には、このRコードを生成するPerlのプログラムを作成してからそれをRで実行していました。現在よく使われているGSEA(Gene Set Enrichment Analysis)のハシリで、遺伝子にアノテーションされたGene Ontologyの高次(根っこに近い)タームで集計してその分布を見ていました)


* 2. 階層的クラスタリング [#z73b23bc]
-【実習2】Rを起動し以下のファイルに記述されたRのコマンド(''02hclust.r'')を実行しなさい。データファイルとして必要な''matrix.txt''をgithubのサイトから[[ダウンロード>https://raw.githubusercontent.com/bonohu/AJACS51/master/matrix.txt]]し、現在作業しているディレクトリにおいておく必要があります。
-【実習2】Rを起動し以下のファイルに記述されたRのコマンド(''02hclust.r'')を実行しなさい。データファイルとして必要な''matrix.txt''をgithubのサイトから[[ダウンロード>https://raw.githubusercontent.com/bonohu/AJACS51/master/matrix.txt]]し、現在作業しているディレクトリにおいておく必要があります。実行後、そのディレクトリに''hclust.png''というファイルが新たに生成されていることを確認しなさい。

 png("hclust.png")
 d <- read.table('matrix.txt')
 c <- hclust(as.dist(d), method = 'average')
 plot(c, hang=-1)
 dev.off()

#fold(←左の+マークをクリックするとこのスクリプトの経緯が…,このプログラムはかつての蛋白質・核酸・酵素に寄稿したレビューで紹介したもので、当時階層的クラスタリングをフリーウェアで実行する手段として重宝された(はず)。[[Bono H and Nakao MC, PNE 2003>http://www.ncbi.nlm.nih.gov/pubmed/12638180]])

* 3. Affymetrixデータ正規化 [#uca2d586]
このパートは以下のブログを参考に。 http://bonohu.jp/blog/2013/06/17/justrma/ 

-【実習3】Bioconductorを利用しましょう。Rを起動して以下のコマンドでaffy libraryをインストールしなさい。
 source("http://bioconductor.org/biocLite.R")
 biocLite("affy")
そして、affy library中のjustRMA関数を利用して、[[RMA>http://allie.dbcls.jp/short/exact/Any/RMA.html]](Robust Multichip Average)正規化してみましょう。自らのAffymetrix Genechip データ(CELファイル)がある人はそれを、ない人は[[GSE17264 Comparative transcriptome analysis of dedifferentiation in porcine mature adipocytes and follicular granulosa cells>http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17264]]にある生データ(CELファイル)で実行しましょう。実行する際、Session → Set Working Directory を、CELファイルをダウンロードしてきたディレクトリに指定することがポイントです。それができたら、以下のコード(''03justrma.r'')を実行します。
 source("http://bioconductor.org/biocLite.R")
 biocLite("affy")
 library(affy)
 write.exprs(justRMA(), file="RMA.txt")
その結果生成される''RMA.txt''ファイルがRMAによって正規化された遺伝子発現データになります。

-【発展2】上記の実習のサンプルデータは [[GPL3533 [Porcine] Affymetrix Porcine Genome Array>http://www.ncbi.nlm.nih.gov/gds?term=GPL3533%5BAccession%5D&cmd=search]] と呼ばれる(カタログ)マイクロアレイ(プラットフォーム)を使ったデータでした。同じマイクロアレイ(プラットフォーム)のデータであればjustRMAを実行することが出来ます。同じプラットフォームのデータの中からiPS細胞のものを探しだして上記のデータに混ぜてjustRMAを実行しなさい。
#fold(←こたえ(の一つ)は左の+マークをクリックすると出てきます,GSE15472 Induced Pluripotent Stem Cells from the Pig Somatic Cells の3つのサンプル(GSM388154,GSM388155,GSM388156)がそれですhttp://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE15472 CELファイルをダウンロードして同じdirectoryに置き、再度justRMAを実行してみましょう)

#fold(←左の+マークをクリックするとこの節の経緯が... ,このマイクロアレイデータは以下の論文で我々が発表したもので、全く同じ手段で正規化を行いました。また、発展課題にある公共遺伝子発現データを足したデータ解釈はやはりこの論文で実際に発表したもので、我々のマイクロアレイデータの生物学的解釈(DFATがMA(Mature AdipocyteよりもiPSに近い)に役だっています [[Ono H, Oki Y, Bono H, Kano K. BBRC 2011>http://www.ncbi.nlm.nih.gov/pubmed/21419102]])

* 4. 主成分分析(PCA) [#od88c1b7]
[[PCA>http://allie.dbcls.jp/short/exact/Any/PCA.html]](Primary Component Analysis)。
このパートは以下のブログを参考に。http://bonohu.jp/blog/2013/07/13/pcaonr/

-【実習4】Rを起動し以下のファイルに記述されたRのコマンド(04pca.r)を実行しなさい。データファイルとして必要なRMA.txtは実習3で作成したものを利用し、現在作業しているディレクトリにおいておく必要があります(うまく行かなかった人はUSBメモリでデータをお渡しします)。
 data <- read.table("RMA.txt", header=TRUE, row.names=1, sep="\t", quote="") 
 data.pca <- prcomp(t(data))
 names(data.pca)
 plot(data.pca$sdev, type="h", main="PCA s.d.")
 data.pca.sample <- t(data) %*% data.pca$rotation[,1:2]
 plot(data.pca.sample, main="PCA")  
 text(data.pca.sample, colnames(data), col = c(rep("red", 3), rep("blue",3),rep("green",3),rep("black",3)))

図に''.CEL.gz''の文字がいっぱい入っていて見づらい?そう思った方はこちらを参照 http://bonohu.jp/blog/2014/09/10/afterjustrma/


* 5. NGS解析(cuffdiffの結果可視化) [#xbad7386]
cufflinksパッケージにcuffdiffという発現データの差分を計算するプログラムがあります。それを実行したあとのデータを可視化する手段として使うBioconductorのパッケージにcummeRbundがあります([[cuffdiffの使用例>http://dx.doi.org/10.1371%2Fjournal.pone.0104283]])。

参考:[[Mishima.syk#4>http://connpass.com/event/8978/]]の[[bonohuの発表資料「NGS解析(RNAseq)」>http://figshare.com/articles/NGS_RNAseq_/1216717]]

-【発展3】Rを起動し以下のファイルに記述されたRのコマンドを実行しなさい。データファイルとして必要なcuffdiffディレクトリ以下のファイルは手持ちのものか、なければサンプルをUSBメモリでお渡しします。現在作業しているディレクトリにおいておく必要があります。
 source("http://bioconductor.org/biocLite.R")
 biocLite("cummeRbund")
 library("cummeRbund")
 cuff.dir <- "cuffdiff"
 cuff <- readCufflinks(dir=cuff.dir)
準備はここまでで、以下のコマンドで発現密度分布のプロット
 dens <- csDensity(genes(cuff))
 dens
CSV形式でFPKM値を出力などができます。
 gene.matrix <- fpkmMatrix(genes(cuff))
 write.csv(gene.matrix, file="fpkm.csv")

* 6.最後に: R Graphical Manualのすゝめ [#sc4355c8]
どんなライブラリがあって、どういった可視化ができるかは[[R Graphical Manual>http://rgm3.lab.nig.ac.jp/RGM]]が大変便利です。
-【実習5】ウェブブラウザを起動し、インターネット検索で''R Graphical Manual''を検索してサイトに辿り着き、サイト内の検索フォームからキーワード''cummeRbund''で検索しなさい。どのような結果が返ってきたか?前節で紹介したcsDensityのエントリを見つけてみなさい。
-【発展4】csDensityのエントリにあるサンプルコードを参考に、''csVolcano''や''dispersionPlot''などのエントリ中のサンプルコードを現在の例に合うように改変して実行してみなさい。