講習会のページに戻る

RおよびBioconductorを使ったバイオデータ解析

講習では【実習】を皆さんとやっていきます。早くできてしまった人は、【発展】をやってみたりしてください。 Rのインストールは以下の統合TVを参考に行っておいてください。

• 参考：統合TV
  • 統計解析ソフト「R」の使い方 導入編
  • 統計解析ソフト「R」の使い方 導入編(MacOSX)
  • 統計解析ソフト「R」の使い方 正規化編
  • 統計解析ソフト「R」の使い方 ヒートマップ編
  • 統計解析ソフト「R」での立廻り

目次

_ 0. はじめに

• 【実習0】Rをインストールしなさい。

  ←こたえは左の+マークをクリックすると出てきます

  http://cran.ism.ac.jp/ からダウンロードしてインストールしてください

• 【発展0】Rを快適に使うにはRstudioというソフトウェアがある。Rstudioをググッて見つけ、インストールしなさい。

  ←こたえは左の+マークをクリックすると出てきます

  http://www.rstudio.com/ からダウンロードしてインストールしてください

• 【発展1】今回使うRのコードやデータはgithubのサイトにあります。そこのデータをまとめてダウンロードするには、UNIXのコマンドラインで

  git clone https://github.com/bonohu/AJACS51

  を実行しなさい。現在居るディレクトリ(current working directory)にどういった変化が起きたか?

  ←こたえは左の+マークをクリックすると出てきます

  AJACS51というディレクトリが作成され、それ以下にgithubにあるファイルがすべてダウンロードされた

_ 1. パイチャート

• 【実習1】Rを起動し以下のファイルに記述されたRのコマンド(01piechart.r)を実行しなさい。データファイルとして必要なsrabystudy.txtをgithubのサイトからダウンロードし、現在作業しているディレクトリにおいておく必要があります。実行後、そのディレクトリに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()

_ 2. 階層的クラスタリング

• 【実習2】Rを起動し以下のファイルに記述されたRのコマンド(02hclust.r)を実行しなさい。データファイルとして必要なmatrix.txtをgithubのサイトからダウンロードし、現在作業しているディレクトリにおいておく必要があります。実行後、そのディレクトリにhclust.pngというファイルが新たに生成されていることを確認しなさい。

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

_ 3. Affymetrixデータ正規化

このパートは以下のブログを参考に。 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(Robust Multichip Average)正規化してみましょう。自らのAffymetrix Genechip データ(CELファイル)がある人はそれを、ない人はGSE17264 Comparative transcriptome analysis of dedifferentiation in porcine mature adipocytes and follicular granulosa cellsにある生データ(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 と呼ばれる(カタログ)マイクロアレイ(プラットフォーム)を使ったデータでした。同じマイクロアレイ(プラットフォーム)のデータであればjustRMAを実行することが出来ます。同じプラットフォームのデータの中からiPS細胞のものを探しだして上記のデータに混ぜてjustRMAを実行しなさい。

  ←こたえ(の一つ)は左の+マークをクリックすると出てきます

  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を実行してみましょう

_ 4. 主成分分析(PCA)

PCA(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の結果可視化)

cufflinksパッケージにcuffdiffという発現データの差分を計算するプログラムがあります。それを実行したあとのデータを可視化する手段として使うBioconductorのパッケージにcummeRbundがあります(cuffdiffの使用例)。

参考：Mishima.syk#4のbonohuの発表資料「NGS解析(RNAseq)」

• 【発展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のすゝめ

どんなライブラリがあって、どういった可視化ができるかはR Graphical Manualが大変便利です。

• 【実習5】ウェブブラウザを起動し、インターネット検索でR Graphical Manualを検索してサイトに辿り着き、サイト内の検索フォームからキーワードcummeRbundで検索しなさい。どのような結果が返ってきたか?前節で紹介したcsDensityのエントリを見つけてみなさい。
• 【発展4】csDensityのエントリにあるサンプルコードを参考に、csVolcanoやdispersionPlotなどのエントリ中のサンプルコードを現在の例に合うように改変して実行してみなさい。