Hatena::Groupbioconductor

Bioconductorノート このページをアンテナに追加 RSSフィード

2007-06-18Bioconductorチュートリアル @ NIG

はじめに はじめに - Bioconductorノート を含むブックマーク はてなブックマーク - はじめに - Bioconductorノート はじめに - Bioconductorノート のブックマークコメント

目次

  • はじめに
    • このチュートリアルの概略
    • 参考資料
  • R の導入について
    • Windowsの場合
    • Mac OS Xの場合
    • Unix系OSの場合
    • R による DDBJ のアクセス方法(なかまさん)
    • Bioconductor の導入
  • Bioconductor の紹介
    • Bioconductor のゴール
    • Bioconductor の特徴
    • スクリーンショット
    • Bioconductor Task Views: BiocViews
    • Bioconductor にできること
    • Bioconductor にできないこと
    • もっと知りたいとき
      • R のヘルプ機能
      • Bioconductor ヴィネット
      • R のグラフサンプル集
      • R と Bioconductor のパッケージソースコード検索
      • 書籍やウェブサイト
      • サーチエンジンでRの話題を検索するコツ
    • Bioconductor 以外のバイオ系Rパッケージ集
      • CRAN Task Views
      • 遺伝学関連
      • 生態学関連
  • ケーススタディ:遺伝子発現データに基づいた分類
    • ストーリー
    • 腎細胞ガンデータの読み込み
    • 分類器の訓練と評価の準備をする
    • 変数選択で利用する遺伝子を選択する
    • 分類を実行
    • 複数回のランダム分割
    • 分類予測
    • まとめ
  • ケーススタディ:CELファイルから興味ある遺伝子リストの作成
    • ストーリー
    • データ読み込み
    • 前処理
    • 遺伝子の順位付け
      • 遺伝子の順位づけのための統計検定量の比較
      • カットオフの選択
    • 関連するアノテーションの付加
      • 統計検定量の付加
      • PubMedアブストラクトの付加
    • レポートの生成
    • まとめ

このチュートリアルは、生物学研究者向けにBioconductorのできることとできないことを概説し、遺伝子発現データからの解析のケーススタディを紹介します。各話題は一時間程度で終わるようになっています。また、R 2.4Bioconductor 1.9 で実行可能なコードを用いています。


このチュートリアルは、「RとBioConductorを用いたバイオインフォマティクス」(シュプリンガージャパン 7月出版予定)にもとづいています。


このチュートリアルの概略

  1. R の導入についてではRとBioconductorの導入方法について説明します。
  2. Bioconductor の紹介では、Bioconductorのできることとできないこと、簡単なBioconductorの使用例、そしてこのチュートリアルのあとに見るべき資料について説明します。
  3. ケーススタディ:遺伝子発現データの分類では、遺伝子発現データを用いた典型的な解析のひとつである、遺伝子発現からの細胞(フェノタイプ)の分類の研究について例を示して、データから分類結果までの流れとそこで用いられる Bioconductor の機能について説明します。
  4. ケーススタディ:CELファイルから興味ある遺伝子リストの作成では、Affymetrix DNA chip解析データCELファイル)から注目している条件下で発現変動を示す遺伝子リスト選別し、その遺伝子リストアノテーションメタデータ)を関連づけし、レポートを作成するまでの解析の流れとそこでもちいられる Bioconductor の機能について説明します。

参考資料

  1. 「RとBioConductorを用いたバイオインフォマティクス」(シュプリンガージャパン 7月出版予定)
  2. Bioconductor - Help
    1. http://www.bioconductor.org/mogr/chapter-code/class.R
    2. http://www.bioconductor.org/mogr/chapter-code/fromcels.R
    3. Bioconductor - Getting the Monograph via R
    4. Bioconductor - Data

screenshot screenshot screenshot


R の導入について R の導入について - Bioconductorノート を含むブックマーク はてなブックマーク - R の導入について - Bioconductorノート R の導入について - Bioconductorノート のブックマークコメント

目次

  • はじめに
    • このチュートリアルの概略
    • 参考資料
  • R の導入について
    • Windowsの場合
    • Mac OS Xの場合
    • Unix系OSの場合
    • R による DDBJ のアクセス方法(なかまさん)
    • Bioconductor の導入
  • Bioconductor の紹介
    • Bioconductor のゴール
    • Bioconductor の特徴
    • スクリーンショット
    • Bioconductor Task Views: BiocViews
    • Bioconductor にできること
    • Bioconductor にできないこと
    • もっと知りたいとき
      • R のヘルプ機能
      • Bioconductor ヴィネット
      • R のグラフサンプル集
      • R と Bioconductor のパッケージソースコード検索
      • 書籍やウェブサイト
      • サーチエンジンでRの話題を検索するコツ
    • Bioconductor 以外のバイオ系Rパッケージ集
      • CRAN Task Views
      • 遺伝学関連
      • 生態学関連
  • ケーススタディ:遺伝子発現データに基づいた分類
    • ストーリー
    • 腎細胞ガンデータの読み込み
    • 分類器の訓練と評価の準備をする
    • 変数選択で利用する遺伝子を選択する
    • 分類を実行
    • 複数回のランダム分割
    • 分類予測
    • まとめ
  • ケーススタディ:CELファイルから興味ある遺伝子リストの作成
    • ストーリー
    • データ読み込み
    • 前処理
    • 遺伝子の順位付け
      • 遺伝子の順位づけのための統計検定量の比較
      • カットオフの選択
    • 関連するアノテーションの付加
      • 統計検定量の付加
      • PubMedアブストラクトの付加
    • レポートの生成
    • まとめ

R は The Comprehensive R Archive Network からソースコードバイナリパッケージが配布されている。

screenshot


Windowsの場合

screenshot


Mac OS Xの場合

screenshot


UnixOSの場合

screenshot screenshot

R による DDBJ のアクセス方法(なかまさん)

BioConductor tutorialWindowsへのインストールについてくわしくかかれています。



Bioconductor の導入

Bioconductorパッケージ集であるので、R のもっているネットワーク機能をつかって導入することができます。

▼R を起動して、コンソール画面でつぎのスクリプト入力します。このスクリプトは、ネットワークからパッケージファイルダウンロードしつつ、パッケージの導入をおこなうので、実行にはPCネットワーク接続している必要があります。

source("http://bioconductor.org/biocLite.R")
biocLite()

このスクリプトによって次のパッケージインストールされる。

  1. affy
  2. affydata
  3. affyPLM
  4. annaffy
  5. annotate
  6. Biobase
  7. Biostrings
  8. DynDoc
  9. gcrma
  10. genefilter
  11. geneplotter
  12. hgu95av2
  13. limma
  14. marray
  15. matchprobes
  16. multtest
  17. reposTools
  18. ROC
  19. vsn
  20. xtable

もくじにもどる


Bioconductor の紹介  Bioconductor の紹介 - Bioconductorノート を含むブックマーク はてなブックマーク -  Bioconductor の紹介 - Bioconductorノート  Bioconductor の紹介 - Bioconductorノート のブックマークコメント

目次

  • はじめに
    • このチュートリアルの概略
    • 参考資料
  • R の導入について
    • Windowsの場合
    • Mac OS Xの場合
    • Unix系OSの場合
    • R による DDBJ のアクセス方法(なかまさん)
    • Bioconductor の導入
  • Bioconductor の紹介
    • Bioconductor のゴール
    • Bioconductor の特徴
    • スクリーンショット
    • Bioconductor Task Views: BiocViews
    • Bioconductor にできること
    • Bioconductor にできないこと
    • もっと知りたいとき
      • R のヘルプ機能
      • Bioconductor ヴィネット
      • R のグラフサンプル集
      • R と Bioconductor のパッケージソースコード検索
      • 書籍やウェブサイト
      • サーチエンジンでRの話題を検索するコツ
    • Bioconductor 以外のバイオ系Rパッケージ集
      • CRAN Task Views
      • 遺伝学関連
      • 生態学関連
  • ケーススタディ:遺伝子発現データに基づいた分類
    • ストーリー
    • 腎細胞ガンデータの読み込み
    • 分類器の訓練と評価の準備をする
    • 変数選択で利用する遺伝子を選択する
    • 分類を実行
    • 複数回のランダム分割
    • 分類予測
    • まとめ
  • ケーススタディ:CELファイルから興味ある遺伝子リストの作成
    • ストーリー
    • データ読み込み
    • 前処理
    • 遺伝子の順位付け
      • 遺伝子の順位づけのための統計検定量の比較
      • カットオフの選択
    • 関連するアノテーションの付加
      • 統計検定量の付加
      • PubMedアブストラクトの付加
    • レポートの生成
    • まとめ

Bioconductor は、統計グラフィックスフリーソフトウェア環境Rのためのゲノムスケールデータの解析パッケージ集です。ゲノムスケールデータとは、マイクロアレイSAGEなどによる大規模遺伝子発現プロファイル、SNPデータ質量分析データゲノム配列遺伝子アノテーションタンパク質相互作用データなどが含まれます。

screenshot screenshot

Bioconductor は R の共同開発者である Prof. Robert Gentlman (Fred Hunchinson Cancer Research Center) を中心にして 2001 年から開発されています。

screenshot screenshot


Bioconductor のゴール

Bioconductorプロジェクトのゴールには次のようなものがあります。

  1. ゲノムスケールデータの解析のために、さまざまな統計手法や可視化手法を使用する手段を提供する
  2. 実験データの解析のために、生物学メタデータPubMedからのデータなど)の統合を容易にする
  3. 拡張可能性、大量データへの対応性そして相互運用性のあるソフトウェアの迅速な開発を可能にする
  4. 品質ドキュメントと再現性のある研究を振興する
  5. ゲノムスケールデータの解析のための計算機的手法・統計手法のトレーニングを提供する

Bioconductor の特徴

Bioconductorプロジェクトの特徴には次のようなものがあります。

  1. R を利用している
  2. 実験の再現に十分なドキュメントを提供している
  3. 統計手法と可視化手法
  4. アノテーションメタデータ
  5. ショートコース
  6. オープンソース
  7. 開かれた開発

スクリーンショット

Bioconductor パッケージスクリーンショットをいくつか紹介します。

screenshot

image(affy) 関数
Affymetrixプローブレベルインテンシティの図示。(AffyBatchクラスのメソッド)

http://bioconductor.org/whatisit/screenshots/affy-image.jpg


mva.pairs(affy) 関数
三つのチップのスキャッターMAプロットMAプロットのM軸は差の対数値、A軸は和の対数値となる。

http://bioconductor.org/whatisit/screenshots/affy-maplot.jpg


ReadAffy(affy) 関数
Affymetrix CDF、CELファイル入力ウイジット

http://bioconductor.org/whatisit/screenshots/affy-widget.jpg


matbox.start(edd) 関数
直線化した行列の箱ヒゲ

http://bioconductor.org/whatisit/screenshots/edd-boxplot.jpg


geneplotter
二群の平均発現レベルを比較する

http://bioconductor.org/whatisit/screenshots/cColor.jpg


alongChrom(geneplotter) 関数
サンプル x 遺伝子発現行列染色体に沿ったカラーイメージ。それぞれの鎖は分離してある。

http://bioconductor.org/whatisit/screenshots/chrom-image.jpg


alongChromo(geneplotter) 関数
発現測定値の染色体に沿った累積図。

http://bioconductor.org/whatisit/screenshots/chrom-plot.jpg


plotChr(geneplotter) 関数
5'から3'へのセンス/アンチセンスの平滑曲線の図示。(ここではヒト21番染色体を使用している)

http://bioconductor.org/whatisit/screenshots/plotChr.png


maImage(marrayPlots) 関数
Cy3 バックグラウンド・インテンシティの可視化

http://bioconductor.org/whatisit/screenshots/marray-image.jpg


maBoxplot(marrayPlots) 関数
プリントチップグループにおける log_2(R/G) ログ比の箱ヒゲ図。

http://bioconductor.org/whatisit/screenshots/marray-boxplot.jpg


ll.htmlpage(annotate) 関数
LocusLink へのリンク付きの遺伝子リスト HTML レポートページ。

http://www.bioconductor.org/whatisit/Screenshots/genelist.jpg


pm.getabst(annotate) 関数
R で PubMed 検索をした結果アブストラクトを表示している。

http://bioconductor.org/whatisit/screenshots/pm.jpg


widget.marrayRaw(marrayInput) 関数
cDNA 画像解析ファイル(例:GenePix、Spot)のデータ入力用のウイジェット

http://bioconductor.org/whatisit/screenshots/marray-widget.jpg


hdiffplot(hexbin) 関数
二つの六角形ビンの間の個別の細胞プロットしたもの。

http://bioconductor.org/whatisit/screenshots/hexbinDiff.jpg


graph and Rgraphviz
TAPタンパク質複合体データスポークモデルグラフ。graphパッケージで作成し、描画にはRgraphvizパッケージを用いた。

http://bioconductor.org/whatisit/screenshots/spoke.jpg


graph and Rgraphviz
相互作用タンパク質ペアのグラフ細胞周期の発現プロファイルクラスタメンバーにそって色づけされている。graphパッケージで作成し、描画にはRgraphvizパッケージを用いた。

http://bioconductor.org/whatisit/screenshots/litGnc.jpg


graph and Rgraphviz
類似発現プロファイルをもつ遺伝子クラスタリングのための完全結合サブグラフ。graphパッケージで作成し、描画にはRgraphvizパッケージを用いた。

http://bioconductor.org/whatisit/screenshots/clustGallc.jpg


graph and Rgraphviz
インテグリンの介在した細胞接着パスウェイのグラフ表現。ノードには対応する遺伝子セットの発現レベルパイチャートで付加してある。

http://bioconductor.org/whatisit/screenshots/Tcell.png


vExplorer(tkWidgets) 関数
インタラクティブヴィネットブラウザ

http://bioconductor.org/whatisit/screenshots/vExplorer.jpg


Bioconductor Task Views: BiocViews

Bioconductor パッケージを関連する作業分野ごとに階層的に分類したものです。パッケージを探すときの整理棚として利用できます。

screenshot screenshot screenshot screenshot

Bioconductor にできること

  1. Affymetrix GeneChip のデータを簡単に扱えます
  2. KEGGデータを簡単に扱えます
  3. PubMedデータを扱えます
  4. GenBankデータを扱えます
  5. 典型的な解析は非常に簡単にできる
  6. よく使われるデータパッケージ化されているので簡単に利用できる。
  7. アノテーションデータネット越しに取得できます。

Bioconductor にできないこと

  1. かんたんに扱えない遺伝子発現データ形式があります
  2. マウスで操作してグラフ作成できない
  3. ワンクリックデータ解析
  4. 典型的でない解析は簡単にできない



もっと知りたいとき

このチュートリアルが終わったあとに参照するべき機能、書籍ウェブサイトを紹介します。

  1. R のヘルプ機能
  2. Bioconductor ヴィネット
  3. R のグラフサンプル集
  4. R と Bioconductorパッケージソースコード検索
  5. 書籍ウェブサイト
  6. サーチエンジンで R の話題を検索するときのコツ

R のヘルプ機能

R のヘルプ機能を利用してメソッドの説明を引き出したり、パッケージの説明を引き出したり、キーワードでメソッドを検索することができます。

▼たとえば、plot メソッドについて知りたいときは、Rコンソールでつぎのように入力します。

help(plot)

△実行結果。plot メソッドに関するヘルプドキュメントが表示されます。R のヘルプドキュメントはすべておなじ形式になっています。つぎに plot のヘルプドキュメントを示します。

plot                package:graphics                R Documentation

Generic X-Y Plotting

Description:

     Generic function for plotting of R objects.  For more details
     about the graphical parameter arguments, see 'par'.

Usage:

     plot(x, y, ...)

Arguments:

       x: the coordinates of points in the plot. Alternatively, a
          single plotting structure, function or _any R object with a
          'plot' method_ can be provided.

       y: the y coordinates of points in the plot, _optional_ if 'x' is
          an appropriate structure.

     ...: Arguments to be passed to methods, such as graphical
          parameters (see 'par'). Many methods will accept the
          following arguments:

     'type' what type of plot should be drawn.  Possible types are

             *  '"p"' for *p*oints,

             *  '"l"' for *l*ines,

             *  '"b"' for *b*oth,

             *  '"c"' for the lines part alone of '"b"',

             *  '"o"' for both "*o*verplotted",

             *  '"h"' for "*h*istogram" like (or "high-density")
                vertical lines,

             *  '"s"' for stair *s*teps,

             *  '"S"' for other *s*teps, see _Details_ below,

             *  '"n"' for no plotting.

          All other 'type's give a warning or an error; using, e.g.,
          'type = "punkte"' being equivalent to 'type = "p"' for S
          compatibility.

     'main' an overall title for the plot: see 'title'.

     'sub' a sub title for the plot: see 'title'.

     'xlab' a title for the x axis: see 'title'.

     'ylab' a title for the y axis: see 'title'.

     'asp' the y/x aspect ratio, see 'plot.window'. 

Details:

     For simple scatter plots, 'plot.default' will be used. However,
     there are 'plot' methods for many R objects, including
     'function's, 'data.frame's, 'density' objects, etc.  Use
     'methods(plot)' and the documentation for these.

     The two step types differ in their x-y preference: Going from
     (x1,y1) to (x2,y2) with x1 < x2, 'type = "s"' moves first
     horizontal, then vertical, whereas 'type = "S"' moves the other
     way around.

See Also:

     'plot.default', 'plot.formula' and other methods; 'points',
     'lines', 'par'.

Examples:

     plot(cars)
     lines(lowess(cars))

     plot(sin, -pi, 2*pi)

     ## Discrete Distribution Plot:
     plot(table(rpois(100,5)), type = "h", col = "red", lwd=10,
          main="rpois(100,lambda=5)")

     ## Simple quantiles/ECDF, see ecdf() {library(stats)} for a better one:
     plot(x <- sort(rnorm(47)), type = "s", main = "plot(x, type = \"s\")")
     points(x, cex = .5, col = "dark red")

「Usage:」と「Examples:」をざっとみてメソッドの使い方の雰囲気をつかみ、必要に応じて「Arguments:」や「Details:」、「See Also:」を参考しするとよい。「Examples:」の内容はRコンソールにそのままコピーペーストすれば実行可能なものです。

パッケージの内容を調べるには、library メソッドを使用します。次の例では annaffy パッケージについて調べます。

library(help=annaffy)

△実行結果


		パッケージ 'annaffy' の情報

記述:

Package:       annaffy
Version:       1.6.2
Date:          2007-02-20
Title:         Annotation tools for Affymetrix
               biological metadata
Author:        Colin A. Smith
               <annaffy@colinsmith.org>
Maintainer:    Colin A. Smith
               <annaffy@colinsmith.org>
Depends:       R (>= 2.1.0), methods, Biobase, GO (>=
               1.6.0), KEGG
Suggests:      hgu95av2, multtest, tcltk
Description:   Functions for handling data from
               Bioconductor Affymetrix annotation
               data packages. Produces compact HTML
               and text reports including
               experimental data and URL links to
               many online databases. Allows
               searching biological metadata using
               various criteria.
License:       LGPL
SaveImage:     yes
biocViews:     OneChannel, Microarray, Annotation,
               GO, Pathways, ReportWriting
Packaged:      Tue Feb 20 13:28:36 2007; biocbuild
Built:         R 2.4.1; ; 2007-02-21 06:14:43; unix

索引:

aaf.handler             Handle feching annotation data columns
aafChromLoc             Constructor for aafChromLoc objects
aafChromLoc-class       Class aafChromLoc, a class for gene chromosome
                        locations
aafChromosome           Constructor for aafChromosome objects
aafChromosome-class     Class aafChromosome, a class for gene
                        chromosome assignments
aafCytoband             Constructor for aafCytoband objects
aafCytoband-class       Class aafCytoband, a class for cytoband data
aafDescription          Constructor for aafDescription objects
aafDescription-class    Class aafDescription, a class for gene
                        descriptions
aafExpr                 Sample exprSet used for demonstration purposes
aafFunction             Constructor for aafFunction objects
aafFunction-class       Class aafFunction, a class for gene product
                        functions
aafGO                   Constructor for aafGO objects
aafGO-class             Class aafGO, a class for gene ontology ids
aafGOItem-class         Class aafGOItem, a class for gene ontology id
                        elements
aafGenBank              Constructor for aafGenBank objects
aafGenBank-class        Class aafGenBank, a class for GenBank accession
                        numbers
aafIntensity-class      Class aafIntensity, a class for gene expression
                        values
aafList-class           Class aafList, a specialized subclass of list
aafLocusLink            Constructor for aafLocusLink objects
aafLocusLink-class      Class aafLocusLink, a class for LocusLink ids
aafPathway              Constructor for aafPathway objects
aafPathway-class        Class aafPathway, a class for KEGG pathway ids
aafPathwayItem-class    Class aafPathwayItem, a class for KEGG pathway
                        id elements
aafProbe                Constructor for aafProbe objects
aafProbe-class          Class aafProbe, a class for Probe ids
aafPubMed               Constructor for aafPubMed objects
aafPubMed-class         Class aafPubMed, a class for PubMed ids
aafSearchGO             Find probe ids corresponding to GO ids
aafSearchText           Search metadata annotation text
aafSigned-class         Class aafSigned, a class for signed numerical
                        data
aafSymbol               Constructor for aafSymbol objects
aafSymbol-class         Class aafSymbol, a class for gene symbols
aafTable                Constructor for aafTable objects
aafTable-class          Class aafTable, a tabular microarray data class
aafTableAnn             Constructor for aafTable objects from
                        annotation data
aafTableFrame           Constructor for aafTable objects from data
                        frames
aafTableInt             Constructor for aafTable objects from exprSets
aafUniGene              Constructor for aafUniGene objects
aafUniGene-class        Class aafUniGene, a class for UniGene cluster
                        ids
getCSS-methods          Methods for function getCSS
getHTML-methods         Methods for function getHTML
getTD-methods           Methods for function getTD
getText-methods         Methods for function getText
getURL-methods          Methods for function getURL
is.annpkg               Determine if packages contain annotation
selectorWidget          Dialog to select items from a list

これ以上の情報はディレクトリ
'/Library/Frameworks/R.framework/Resources/library/annaffy/doc'
にある以下の vignettes 中にあります:

annaffy: annaffy Primer (source, pdf)

索引:」以下がこのパッケージに含まれるメソッドの一覧になります。


キーワードからメソッドを検索する。

help.search("MA plot")

Bioconductor ヴィネット

Bioconductor では R のヘルプドキュメントをより高品質化をおしすすめ、背景知識を含めたレポート形式のドキュメント(通称:ヴィネット、vignette)を提供しています。

ヴィネットの表示。調べたいトピック(例:grid)の文章をひらきます。

vignette("grid")

△実行結果PDF書類が表示されます。

f:id:nakao_mitsuteru:20070617155541p:image

vignette でつかえるトピックは、annaffy パッケージの場合は library(help=annaffy) の最後に書いてあります。

▼ vignette であつかっているトピックリストはつぎのようにして表示します。

vignette()

R のグラフサンプル集

no title では R のパッケージの Example にあるグラフとそのコードがすべて閲覧可能になっています。

screenshot screenshot


また、Statistics with R3.From Data to Graphics4.Customizing graphicsコードと図が豊富にあり参考になります。

screenshot screenshot screenshot


R と Bioconductorパッケージソースコード検索

すこしプログラミングに慣れてくると、ほかの人の書いたプログラムを読んで参考にすることができます。また、実際に動作する他人の書いたプログラムプログラミングの自習に最適な教材です。

バイオインフォマティクスで利用されるソフトウェアソースコード検索サーバ産業技術総合研究所CBRCサービスされています。

screenshot screenshot


Googleソースコード検索では、言語指定(lang:R)すると効果的な検索がおこなえます。

screenshot screenshot


中間さんによる R や Bioconductorソースコード検索サービス

screenshot screenshot screenshot


書籍ウェブサイト

7月出版予定の R と Bioconductor を用いたバイオインフォマティクスシュプリンガージャパン)はつぎの本の翻訳です。

プログラミング言語 R の基礎についてはつぎの本が良いです。

Rの基礎とプログラミング技法

Rの基礎とプログラミング技法

統計解析の基礎についてはつぎの本が良いです。

Rによる統計解析の基礎 (Computer in Education and Research)

Rによる統計解析の基礎 (Computer in Education and Research)

分子進化の解析にはつぎの本を読むと良いでしょう。

Analysis of Phylogenetics And Evolution With R (Use R)

Analysis of Phylogenetics And Evolution With R (Use R)

RjpWiki に幅広い話題がかかれています。ウェブで R の話題を調べるときのスタート地点にすると良いでしょう。

screenshot


サーチエンジンでRの話題を検索するコツ

R は文字列として短すぎるので、R だけでは Google などの検索エンジンでは意味のある検索結果は期待できません。

仮説:R を話題にしているページは R のオフィシャルウェブサイトリンクしている
これを期待するならば、link:www.r-project.org を検索に加えると良いでしょう。

Bioconductor 以外のバイオ系Rパッケージ

Bioconductor 以外にも生物学(遺伝学、進化系統解析など)に有用な R パッケージがあります。


CRAN Task Views

screenshot

膨大な R パッケージが代表的な研究テーマごとに整理されています。


遺伝学関連

screenshot


生態学関連

screenshot



もくじにもどる



ケーススタディ遺伝子発現データに基づいた分類  ケーススタディ:遺伝子発現データに基づいた分類 - Bioconductorノート を含むブックマーク はてなブックマーク -  ケーススタディ:遺伝子発現データに基づいた分類 - Bioconductorノート  ケーススタディ:遺伝子発現データに基づいた分類 - Bioconductorノート のブックマークコメント

目次

  • はじめに
    • このチュートリアルの概略
    • 参考資料
  • R の導入について
    • Windowsの場合
    • Mac OS Xの場合
    • Unix系OSの場合
    • R による DDBJ のアクセス方法(なかまさん)
    • Bioconductor の導入
  • Bioconductor の紹介
    • Bioconductor のゴール
    • Bioconductor の特徴
    • スクリーンショット
    • Bioconductor Task Views: BiocViews
    • Bioconductor にできること
    • Bioconductor にできないこと
    • もっと知りたいとき
      • R のヘルプ機能
      • Bioconductor ヴィネット
      • R のグラフサンプル集
      • R と Bioconductor のパッケージソースコード検索
      • 書籍やウェブサイト
      • サーチエンジンでRの話題を検索するコツ
    • Bioconductor 以外のバイオ系Rパッケージ集
      • CRAN Task Views
      • 遺伝学関連
      • 生態学関連
  • ケーススタディ:遺伝子発現データに基づいた分類
    • ストーリー
    • 腎細胞ガンデータの読み込み
    • 分類器の訓練と評価の準備をする
    • 変数選択で利用する遺伝子を選択する
    • 分類を実行
    • 複数回のランダム分割
    • 分類予測
    • まとめ
  • ケーススタディ:CELファイルから興味ある遺伝子リストの作成
    • ストーリー
    • データ読み込み
    • 前処理
    • 遺伝子の順位付け
      • 遺伝子の順位づけのための統計検定量の比較
      • カットオフの選択
    • 関連するアノテーションの付加
      • 統計検定量の付加
      • PubMedアブストラクトの付加
    • レポートの生成
    • まとめ

つぎのclass.Rを参考にして、遺伝子発現データに基づいた分類についてのケーススタディを説明します。

▼のついてソースコードをRコンソールにコピーペーストしながらすすめていきます。

ストーリー

  1. データ:腎細胞ガン遺伝子発現データ
  2. 目的遺伝子発現データから細胞の表現形(ガンの種類)を分類予測する
  3. 目的2:うまく分類できる分類法を模索する
  4. データの読み込み
  5. 変数選択で利用する遺伝子を選択
  6. 複数の分類方法(DLDA、kNN、SVM)を評価
  7. 複数回のランダム分割によるパラメータ最適化
  8. 細胞の分類予測をおこなう
  9. 出口:適切な分類方法とそれによる細胞の分類


細胞ガンデータの読み込み

GeneChip による腎臓ガン細胞遺伝子発現データは、kidpack データパッケージの eset データを使用します。

データの読み込み。

set.seed(32)
library("kidpack")
data(eset)

▼あとで利用するためにデータセットから一部取り除いておきます。

x        <- t(exprs(eset)) 
y        <- pData(phenoData(eset))$type
y        <- as.factor(y)
y.train  <- y[-c(1:10) ]
x.train  <- x[-c(1:10),]  

▼eset の中身

eset

△実行結果

ExpressionSet (storageMode: lockedEnvironment)
assayData: 4224 features, 74 samples 
  element names: exprs 
phenoData
  sampleNames: 1, 2, ..., 87 (74 total)
  varLabels and varMetadata:
    lfdnr: NULL
    patientid: NULL
    ...: ...
    22q: NULL
    (64 total)
featureData
  featureNames: 1, 2, ..., 4224 (4224 total)
  varLabels and varMetadata: none
experimentData: use 'experimentData(object)'
Annotation [1] ""

▼eset は Biobase パッケージに含まれる exprSet クラスデータです。

class(eset)

△実行結果

[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"

細胞型は $type でアクセスできます。

eset$type

△実行結果

 [1] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "pRCC"  "ccRCC"
 [8] "pRCC"  "chRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "chRCC"
[15] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "pRCC" 
[22] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[29] "pRCC"  "ccRCC" "ccRCC" "ccRCC" "pRCC"  "ccRCC" "chRCC"
[36] "chRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "pRCC"  "pRCC" 
[43] "ccRCC" "pRCC"  "pRCC"  "ccRCC" "pRCC"  "ccRCC" "ccRCC"
[50] "pRCC"  "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[57] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "pRCC" 
[64] "ccRCC" "chRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[71] "chRCC" "pRCC"  "ccRCC" "ccRCC"

▼三種類の細胞型のそれぞれの数の分布をみます。

table(eset$type)

△実行結果

ccRCC chRCC  pRCC 
   52     9    13

分類器の訓練と評価の準備をする

分類には、訓練データ学習データと評価データの三種類のデータを用意する必要があります。それぞれを学習データと評価データは訓練データ2:1 に分割して作成されます。

▼訓練データサイズ(t.size)、学習データサイズ(l.size)、評価データサイズ(v.size)

t.size   <- length(y.train)
l.size   <- round((2/3) * t.size) # 2/3 が l.size へ
v.size   <- t.size - l.size       # 1/3 が v.size へ

▼訓練データ(t.samp)、学習データ(l.samp)、評価データ(v.samp)

l.samp   <- NULL
K        <- nlevels(y.train)
props    <- round(l.size / t.size * table(y.train))
props[1] <- l.size - sum(props[2:K])
for (k in 1:K)
  {
    y.num  <- as.numeric(y.train)
    l.samp <- c(l.samp, sample(which(y.num == k))[1:props[k]])
  }
v.samp     <- (1:t.size)[-l.samp]

データ集合をみる

学習集合(l.samp)を表示してみる。

table(y.train[l.samp])

△実行結果

ccRCC chRCC  pRCC 
   31     5     7 

▼評価集合(v.samp)を表示してみる。

table(y.train[v.samp])

△実行結果

ccRCC chRCC  pRCC 
   14     3     4

このようにデータの比率は保存されています。


変数選択で利用する遺伝子を選択する

効果的な分類のために分類予測に使用する遺伝子を選別します。ここでは単純なF統計量(群間平方和と群内平方和の比)を用います。F統計量の上位200遺伝子を選択します。

library("multtest")
yl.num     <- as.numeric(y.train[l.samp]) - 1
xl.mtt     <- t(x.train[l.samp,])
f.stat     <- mt.teststat(xl.mtt, yl.num, test = "f")
best.genes <- rev(order(f.stat))[1:200]
x.learn    <- x.train[l.samp, best.genes]
x.valid    <- x.train[v.samp, best.genes]
y.learn    <- y.train[l.samp]
y.valid    <- y.train[v.samp]

遺伝子発現のF統計量が大きいとは、群内での分散が小さくかつ群間での分散が大きい、すなわち、群(ここでは分類対象:細胞タイプ)に特徴的な遺伝子であるとしています。


分類を実行

分類には対角線形判別(DLDA)、k最近傍分類(kNN)、サポートベクターマシンSVM)の三種類のアルゴリズムを試してみます。

▼対角線形判別(diagonal linear disriminant analysis: DLDA)の実行

library("sma")
yl.num      <- as.numeric(y.learn)
yv.num      <- as.numeric(y.valid) - 1
dlda.predic <- stat.diag.da(x.learn, yl.num, x.valid)
dlda.error  <- mean(dlda.predic$pred - 1 != yv.num)

▼k最近傍分類(kNN)の実行。kNNではパラメータ k を 1, 3, 5 に振って最適な分類を探索します。

library("class")
knn.error <- numeric(3)
for(k in c(1,3,5))
  {
    i              <- ((k - 1) / 2) + 1
    knn.predic     <- knn(x.learn, x.valid, y.learn, k = k)
    knn.error[i]   <- mean(knn.predic != y.valid)
  }

サポートベクターマシンSVM)の実行。SVM では二つのパラメータ(γと c)を振って最適なパラメータを探索します。パラメータγ(gamma)は -1, 0, 1に、パラメータc(cost) は 0, 1, 2 に振っています。実際のパラメータは,cost = 2^cost なので、{2^0, 2^1, 2^2}です。また、gamma = 2^gamma/ncol(x.learn) なので、{2^(-1)/200, 2^0/200, 2^1/200} です。

library("e1071")
svm.error <- matrix(0, nrow = 3, ncol = 3)
for(cost in 0:2)
  {
    for(gamma in (-1):1)
      {
        i              <- cost + 1
        j              <- gamma + 2
        svm.fit        <- svm(x.learn, y.learn, cost = 2^cost,    
                              gamma = 2^gamma/ncol(x.learn))
        svm.predic     <- predict(svm.fit, newdata = x.valid)
        svm.error[i,j] <- mean(svm.predic != y.valid)
      }
  }

分類器 DLDA、kNN、SVM の誤分類率を表示

▼DLDAの誤分類率を表示

dlda.error

△実行結果

[1] 0

▼kNNの誤分類率を表示

knn.error

△実行結果

[1] 0.04761905 0.00000000 0.04761905

パラメータ k = 3 の時がもっとも正確な分類をおこなっています。

SVMの誤分類率を表示

svm.error

△実行結果

           [,1]       [,2]      [,3]
[1,] 0.04761905 0.04761905 0.0952381
[2,] 0.04761905 0.04761905 0.0952381
[3,] 0.04761905 0.04761905 0.0952381

行がパラメータc {0,1,2} で列がパラメータγ{-1,0,1}です。

c が {0,1} かつ γ{-1,0} の組み合わせが、SVM のなかでもっとも正確な分類をおこなっています。

ここで三つの分類法を比較しました。その結果、DLDA と k = 3 の kNN がもっとも正確な分類をおこなうことがわかりました。今後の解析(診断)ではどちらかの方法をつかうとよいといえます。

しかし、今回の訓練データの分割がたまたま DLDA や k = 3 の kNN に有利ダッタのかもしれません。このような偶然性への依存を減らすには訓練データの分割を複数回おこなうことが有効です。


複数回のランダム分割

上記のようなパラメータ最適化をおこなう場合は、データの分割を複数回ためす必要があります。ここでは、訓練データの50 回のランダム分割によるパラメータ最適化をみてみます。

randiv.repr <- function(x, y)
  {
    t.size   <- length(y)
    l.size   <- round((2 / 3) * t.size)
    v.size   <- t.size - l.size
    l.samp   <- NULL
    K        <- nlevels(y)
    props    <- round(l.size / t.size * table(y))
    props[1] <- l.size - sum(props[2:K])
    for (k in 1:K)
      {
        y.num   <- as.numeric(y)
        l.samp  <- c(l.samp, sample(which(y.num == k))[1:props[k]])
      }
    v.samp      <- (1:t.size)[-l.samp]
    yl.num      <- as.numeric(y[l.samp]) - 1
    f.stat      <- mt.teststat(t(x[l.samp,]), yl.num, test = "f")
    best.genes  <- rev(order(f.stat))[1:200]
    x.learn     <- x[l.samp, best.genes]
    x.valid     <- x[v.samp, best.genes]
    y.learn     <- y[l.samp]
    y.valid     <- y[v.samp]
    yl.num      <- as.numeric(y.learn)
    yv.num      <- as.numeric(y.valid) - 1
    dlda.predic <- stat.diag.da(x.learn, as.numeric(y.learn), x.valid)
    dlda.error  <- mean(dlda.predic$pred-1 != yv.num)
    knn.error   <- numeric(3)
    for(k in c(1,3,5))
      {
        i              <- ((k - 1) / 2) + 1
        knn.predic     <- knn(x.learn, x.valid, y.learn, k = k)
        knn.error[i]   <- mean(knn.predic != y.valid)
      }
    svm.error <- matrix(0, nrow = 3, ncol = 3)
    for(cost in 0:2)
      {
        for(gamma in (-1):1)
          {
            i              <- cost + 1
            j              <- gamma + 2
            svm.fit        <- svm(x.learn, y.learn, cost = 2^cost,    
                                  gamma = 2^gamma / ncol(x.learn))
            svm.predic     <- predict(svm.fit, newdata = x.valid)
            svm.error[i,j] <- mean(svm.predic != y.valid)
          }
      } 
    list(dlda = dlda.error, knn = knn.error, svm = svm.error,
         lsamp = l.samp, vsamp = v.samp)
  }

データの分割と分類を50回くり返す

runs         <- 50
dlda.errors  <- numeric(runs)
knn.errors   <- matrix(0, nrow = 3, ncol = runs)
svm.errors   <- array(0, dim = c(3, 3, runs))
vsamp        <- matrix(0, runs, length(v.samp))
lsamp        <- matrix(0, runs, length(l.samp))
for(r in 1:runs)
  {
    results         <- randiv.repr(x.train, y.train)
    dlda.errors[ r] <- results$dlda
    knn.errors[, r] <- results$knn
    svm.errors[,,r] <- results$svm
    vsamp[r,]       <- results$vsamp
    lsamp[r,]       <- results$lsamp
    cat("This was run", r, "of", runs, "\n")
  }

△実行結果

This was run 1 of 50 
This was run 2 of 50 
This was run 3 of 50 
This was run 4 of 50 
This was run 5 of 50 
This was run 6 of 50 
This was run 7 of 50 
This was run 8 of 50 
This was run 9 of 50 
This was run 10 of 50 
This was run 11 of 50 
This was run 12 of 50 
This was run 13 of 50 
This was run 14 of 50 
This was run 15 of 50 
This was run 16 of 50 
This was run 17 of 50 
This was run 18 of 50 
This was run 19 of 50 
This was run 20 of 50 
This was run 21 of 50 
This was run 22 of 50 
This was run 23 of 50 
This was run 24 of 50 
This was run 25 of 50 
This was run 26 of 50 
This was run 27 of 50 
This was run 28 of 50 
This was run 29 of 50 
This was run 30 of 50 
This was run 31 of 50 
This was run 32 of 50 
This was run 33 of 50 
This was run 34 of 50 
This was run 35 of 50 
This was run 36 of 50 
This was run 37 of 50 
This was run 38 of 50 
This was run 39 of 50 
This was run 40 of 50 
This was run 41 of 50 
This was run 42 of 50 
This was run 43 of 50 
This was run 44 of 50 
This was run 45 of 50 
This was run 46 of 50 
This was run 47 of 50 
This was run 48 of 50 
This was run 49 of 50 
This was run 50 of 50 

DLDA、knn、svmにおける誤判別率の平均値を表示します。

▼DLDAの誤判別率の平均値を表示

mean(dlda.errors)

△実行結果

[1] 0.0467

▼kNNの誤判別率の平均値を表示

apply(knn.errors, 1, mean)

△実行結果

[1] 0.0590 0.0571 0.0600

SVMの誤判別率の平均値を表示

apply(svm.errors, c(1,2), mean)

△実行結果

       [,1]   [,2]  [,3]
[1,] 0.0657 0.0810 0.138
[2,] 0.0610 0.0705 0.128
[3,] 0.0562 0.0695 0.128

▼誤判別率の箱ヒゲ図と棒グラフを表示するメソッドを定義。

## Loading the required package
library("grid")

## Check the settings
if (!exists("pushViewport")) pushViewport <- push.viewport
if (!exists("popViewport"))  popViewport  <- pop.viewport

## The function for the barplots
grid.barplot <- function(x, xlab = FALSE, main = FALSE, xscale = NULL, yticks, ylabels)
{
  ## Customizing the data
  y      <- table(x)
  y.alt  <- table(x)
  x.alt  <- x
  x      <- sort(unique(x))
  if(!is.null(xscale))
    {
      index <- which(x > xscale[1] & x < xscale[2])
      x <- x[index]
      y <- y[index]
    }

  ## Horizontal bars
  grid.segments(x0 = unit(x, "native"), y0 = unit(0, "native"),
                x1 = unit(x, "native"), y1 = unit(y.alt, "native"))

  ## Drawing the lines
  grid.lines(unit(x, "native"), unit(y, "native"))

  ## Plotting the mean
  brks   <- as.numeric(attributes(table(x))$dimnames$x)
  start  <- max(which(brks < mean(x.alt)))
  ende   <- start + 1
  xinc   <- (mean(x.alt) - brks[start]) / (brks[ende] - brks[start])
  yh     <- ((y.alt[ende] - y.alt[start])*xinc) + (y.alt[start])             
  farbe  <- gpar(col="red")
  grid.segments(x0 = unit(mean(x.alt), "native"), y0 = unit(0, "native"),
                x1 = unit(mean(x.alt), "native"), y1 = unit(yh, "native"), gp = farbe)

  ## Annotation
  if(xlab) grid.xaxis()
  grid.yaxis(main = FALSE, at = yticks, name = ylabels)
  grid.rect()
}


## The function for the boxplots
grid.boxplot <- function(x, xlab = FALSE, main = FALSE, ylab = NULL, xscale = NULL)
{
  ## Generic boxplot
  x    <- boxplot(x, plot = FALSE)
  farb <- gpar(col="red")

  ## Transformation to the grid-system
  grid.lines(unit(x$stats[1],  "native"), unit(c(0.3, 0.7), "npc"))
  grid.lines(unit(x$stats[1:2],"native"), unit(0.5, "npc"), gp = gpar(lty = 2))
  grid.rect(unit(x$stats[2], "native"), y = unit(0.5, "npc"),
            height = unit(0.6, "npc"), width = unit(diff(x$stats[2:3]), "native"),
            just = c("left", "centre"))
  grid.rect(unit(x$stats[3], "native"), y = unit(0.5, "npc"),
            height = unit(0.6, "npc"), width = unit(diff(x$stats[3:4]), "native"),
            just = c("left", "centre"))
  grid.lines(unit(x$stats[4:5],"native"), unit(0.5, "npc"), gp = gpar(lty = 2))
  grid.lines(unit(x$stats[5], "native"), unit(c(0.3, 0.7), "npc"))
  grid.lines(unit(x$stats[3], "native"), unit(c(0.2, 0.8), "npc"), gp = farb)

  ## Outliers
  n <- length(x$out)
  if(n > 0) {
    if(!is.null(xscale)) index <- which(x$out > xscale[1] & x$out < xscale[2])
      else index <- 1:n
    if(length(index) > 0)
      {
        grid.points(unit(x$out[index],"native"),unit(rep(0.5,length(index)),
                    "npc"), size = unit(0.5, "char"))
      }
  }

  ## Annotation
  if(xlab) grid.xaxis()
  if(!is.null(ylab))
    {
      grid.text(ylab, x = unit(0, "npc") - unit(0.5, "lines"),
                gp = gpar(fontface = "bold"), just = c("right", "centre"))
    }
}


## Plotting the results
grid.boxNbar <- function(daten, mname)
  {
    grid.newpage()
    pushViewport(plotViewport(c(4, 6, 4, 4)))
    #pushViewport(viewport(layout = lyt))
    n         <- ncol(daten)

    ## Calibration of the x-axis
    xsc     <- range(daten) + c(-1, 1) * max(daten)/15
    xsc2    <- xsc          

    ## Calibration of the y-axis
    ysc     <- c(0, max(sapply(apply(daten, 2, function(x) table(x)), max)))
    ysc     <- ysc * 1.1
    ylabels <- c("0", "10", "20", "30", "40", "50")
    ynumb   <- floor((ysc[2] + 3) / 10) - 1
    yticks  <- (0:ynumb) * 10
    ylabels <- ylabels[1:(ynumb + 1)]
          
    ## The plot contains 3 pieces
    lyt <- grid.layout(1, 3, widths = unit(c(0.5 - 0.03 / 2, 0.03, 0.5 - 0.03 / 2), "npc"))
    pushViewport(viewport(layout = lyt))

    ## Left side: boxplots
    lyt <- grid.layout(n, 1, heights = unit(rep(1/n, n), "npc"))
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1, xsc = xsc))
    pushViewport(viewport(layout = lyt))
    grid.rect()
    for(i in n:1)
      {
        pushViewport(viewport(layout.pos.r = i, layout.pos.c = 1, xsc = xsc))
        main  <- (i < 2)
        grid.boxplot(daten[,i], m = main, xl = (i > (n - 1)), yl = mname[i], xsc = xsc)
        popViewport(1)
      }
    popViewport(1)
    grid.text("Boxplot", y = unit(-2.5, "lines"))

    ## Middle: an empty frame  
    popViewport(1)
    #schrift <- gpar(fontface = "bold")
    #grid.text(namen[j], y = unit(1, "npc") + unit(2, "lines"), gp = schrift)

    ## Right side: density plots
    lyt <- grid.layout(n, 1, hei = unit(rep(1/n, n), "npc"))
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 3, xscale = xsc))
    pushViewport(viewport(layout = lyt))
    grid.rect()
    for(i in n:1)
      {
        pushViewport(viewport(layout.pos.r = i, layout.pos.c = 1, xs = xsc2, ys = ysc))
        main  <- (i < 2)
        grid.barplot(daten[,i], main, xsc2, yticks, ylabels, xlab = (i > (n - 1)))
        popViewport(1)
      }
    popViewport(1)
    grid.text("Density", y = unit(-2.5, "lines"))
    popViewport(3)
  }

▼図を生成。

daten <- cbind(dlda.errors, knn.errors[2,], svm.errors[3,1,], svm.errors[2,1,])
mname <- c("DLDA", "kNN", "SVM1", "SVM2")
grid.boxNbar(daten, mname)

△実行結果

f:id:nakao_mitsuteru:20070616213644p:image

kNN は k = 3 の場合、SVM1 は c = 2^2 かつ γ = 2^1/200、SVM2 は c = 2^2 かつ γ = 2^0/200の場合です。赤いバーは平均値を示しています。左のパネルが箱ヒゲ図で分類精度の分布を表していて、DLDA の分類精度がほかのものよりよいことが読み取れる。右のパネルでは密度分布図で分類精度の分布を表している。こちらでも同様である。

補足:図をファイルとして保存する場合は、pdf() と dev.off() を利用します。

pdf(file = 'box_and_barplot.pdf')
grid.boxNbar(daten, mname)
dev.off()

分類予測

はじめにデータから取り除いていた 10 個のサンプルに対しての細胞型の予測をおこないます。

▼DLDA による予測を実行

test   <- (1:10)
train  <- (1:length(eset$type))[-test]  # needed!
trEset <- eset[,train]
library("MLInterfaces")
yt.num      <- as.numeric(factor(trEset$type)) - 1
xt.mtt      <- exprs(trEset)
f.stat      <- mt.teststat(xt.mtt, yt.num, test = "f")
best.genes  <- rev(order(f.stat))[1:200]
selEset     <- eset[best.genes,]
dlda.predic <- stat.diag.daB(selEset, "type", train)
dlda.predic@RObject$pred - 1

△実行結果

 [1] 0 0 0 0 0 2 0 2 1 0

▼DLDA による予測結果を混合行列で評価。

confuMat(dlda.predic)

△実行結果

       predicted
given   1 2 3
  ccRCC 7 0 0
  chRCC 0 1 0
  pRCC  0 0 2

ここでは10個すべてのサンプルで正しい細胞型を予測しました。正しく予測した場合、対角成分に数字がならびます。対角成分以外に数字がある場合は誤分類(偽陽性もしくは偽陰性)の数を示しています。


まとめ

ここでは、遺伝子発現データから細胞を分類予測する研究をみました。

このような解析はあたらしい診断、分類や予測方法の開発でもちいられます。

Bioconductor では数多くの分類法が利用できます。

Bioconductor には遺伝子発現データのサンプルデータ豊富にあります。

複雑なことをするにはそれなりに長いプログラムが必要です。



もくじにもどる


ケーススタディCELファイルから興味ある遺伝子リストの作成  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート を含むブックマーク はてなブックマーク -  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート  ケーススタディ:CELファイルから興味ある遺伝子リストの作成 - Bioconductorノート のブックマークコメント

目次

  • はじめに
    • このチュートリアルの概略
    • 参考資料
  • R の導入について
    • Windowsの場合
    • Mac OS Xの場合
    • Unix系OSの場合
    • R による DDBJ のアクセス方法(なかまさん)
    • Bioconductor の導入
  • Bioconductor の紹介
    • Bioconductor のゴール
    • Bioconductor の特徴
    • スクリーンショット
    • Bioconductor Task Views: BiocViews
    • Bioconductor にできること
    • Bioconductor にできないこと
    • もっと知りたいとき
      • R のヘルプ機能
      • Bioconductor ヴィネット
      • R のグラフサンプル集
      • R と Bioconductor のパッケージソースコード検索
      • 書籍やウェブサイト
      • サーチエンジンでRの話題を検索するコツ
    • Bioconductor 以外のバイオ系Rパッケージ集
      • CRAN Task Views
      • 遺伝学関連
      • 生態学関連
  • ケーススタディ:遺伝子発現データに基づいた分類
    • ストーリー
    • 腎細胞ガンデータの読み込み
    • 分類器の訓練と評価の準備をする
    • 変数選択で利用する遺伝子を選択する
    • 分類を実行
    • 複数回のランダム分割
    • 分類予測
    • まとめ
  • ケーススタディ:CELファイルから興味ある遺伝子リストの作成
    • ストーリー
    • データ読み込み
    • 前処理
    • 遺伝子の順位付け
      • 遺伝子の順位づけのための統計検定量の比較
      • カットオフの選択
    • 関連するアノテーションの付加
      • 統計検定量の付加
      • PubMedアブストラクトの付加
    • レポートの生成
    • まとめ

つぎのfromcels.Rを参考にして、CELファイルから興味ある遺伝子リストの作成についてのケーススタディを説明します。

▼のついたソースコード片をRコンソールにコピーペーストしてすすめます。

ストーリー

  1. データ:ある条件下で発現変動のある遺伝子をねらった Affymetrix GeneChip による遺伝子発現データ
  2. シグナル/ノイズ比を改善するための前処理
  3. 発現変動遺伝子の順位付けのための適切な統計量の選択
  4. 候補遺伝子選択のためのフィルタの決定
  5. 関連するアノテーションの付加
  6. 出口:ウェブブラウザで閲覧できる発現変動のある遺伝子リスト

データ読み込み

library("affy")
library("SpikeInSubset")
data(spikein95)

このケーススタディでは spikein95 という AffyBatch クラスデータを利用します。このデータスパイクイン実験データなので、発現変動している遺伝子あきらかな状態になっています。


前処理

▼このケーススタディ用にデータのラベルを操作します。

pd <- data.frame(population = c(1,1,1,2,2,2),
                 replicate = c(1,2,3,1,2,3))
rownames(pd) <- sampleNames(spikein95)
vl <- list(population = "1 is control, 2 is treatment",
           replicate = "arbitrary numbering") 
phenoData(spikein95) <- new("phenoData", pData = pd, varLabels = vl)

RMAを適用します。

eset <- rma(spikein95)

△実行結果

Background correcting
Normalizing
Calculating Expression

RMAバックグラウンド補正、正規化、プローブレベルデータの要約をおこないます。ここで eset は exprSet クラスデータです。


データ次元を確認

e <- exprs(eset)
dim(e)

△実行結果

[1] 12626     6

行が遺伝子、列がアレイを示しています。

▼それぞれのアレイがどの集団(population)に属しているかを示す

pData(eset)

△実行結果

                 population replicate
1521a99hpp_av06           1         1
1532a99hpp_av04           1         2
2353a99hpp_av08           1         3
1521b99hpp_av06           2         1
1532b99hpp_av04           2         2
2353b99hpp_av08r          2         3

▼集団毎のインデクス・ベクトルを作成

Index1 <- which(eset$population == 1)
Index2 <- which(eset$population == 2)

遺伝子の順位付け

▼各遺伝子の発現変動の対数変換された平均値を計算する

d <- rowMeans(e[, Index2]) - rowMeans(e[, Index1])

▼発現量の対数変換された平均値を計算する

a <- rowMeans(e)

▼発現変動が二倍以上ある遺伝子の数を計算する

sum(abs(d) > 1)

△実行結果

[1] 1

一つしか無かった。


遺伝子の順位づけのための統計検定量の比較

発現変動が興味深い遺伝子を選択するために、いくつかの基準を比較する。

MAプロットデータをみる

plot(a, d, ylim = c(-1, 1), main = "A) MA-plot", pch = ".")

△実行結果

f:id:nakao_mitsuteru:20070616221318p:image


▼t統計量を計算する

library("genefilter")
tt <- rowttests(e, factor(eset$population))

▼t検定の結果をVolcanoプロットでみる

lod <- -log10(tt$p.value)
plot(d, lod, cex =.25, main = "B) Volcano plot for $t$-test")
abline(h = 2)

△実行結果

f:id:nakao_mitsuteru:20070616221317p:image

Volcanoプロットは、(対数オッズ化した)p値と平均の差を同時にみるためのものです。

▼t検定の結果をVolcanoプロットで拡大してみる

o1 <- order(abs(d), decreasing = TRUE)[1:25]
o2 <- order(abs(tt$statistic), decreasing = TRUE)[1:25]
o <- union(o1, o2)
plot(d[-o], lod[-o], cex = .25, xlim = c(-1, 1), ylim = range(lod), main = "C) Close up of B)")
points(d[o1], lod[o1], pch = 18, col = "blue")
points(d[o2], lod[o2], pch = 1, col = "red")

△実行結果

f:id:nakao_mitsuteru:20070616221316p:image

対数スケールの倍数変化の平均によって上位25位までに順位付けられた遺伝子は青いひし形で示し,p値の小さい順に選択した25遺伝子は赤い円で示している。

倍数変化によって選択される遺伝子とt統計量によって選択される遺伝子には比較的大きな違いがあることが読み取れる。


改良されたt統計量のひとつ(moderated t検定)を試す。

limmaを実行する

library("limma")
design <- model.matrix(~factor(eset$population))
fit <- lmFit(eset, design)
ebayes <- eBayes(fit)

▼moderated t検定の結果をVolcanoプロットでみる

lod <- -log10(ebayes$p.value[,2])
mtstat<- ebayes$t[,2]
o1 <- order(abs(d), decreasing = TRUE)[1:25]
o2 <- order(abs(mtstat), decreasing = TRUE)[1:25]
o <- union(o1, o2)
plot(d[-o], lod[-o], cex = .25, xlim = c(-1,1), ylim = c(0,4), main = "D) Volcano plot for moderated $t$-test")
points(d[o1], lod[o1], pch = 18, col = "blue")
points(d[o2], lod[o2], pch = 1, col = "red")

△実行結果

f:id:nakao_mitsuteru:20070616221314p:image


カットオフの選択

順位付けられた遺伝子をどこまで利用するかをきめるカットオフについて考えます。

遺伝子数が多い(12,626個)ので、0.01の有意水準でも 126個の偽陽性が期待されてしまいます。

▼p値が 0.01 以下の遺伝子数を計算

sum(tt$p.value <= 0.01)

△実行結果

[1] 46

比較

データスパイクイン実験なので、発現変動遺伝子の数は既知です。それを取り出して比較します。

data(spikein95)
spikedin <- colnames(pData(spikein95))
spikedIndex <- match(spikedin, geneNames(eset))

スパイクインの16遺伝子がそれぞれの統計量でどのような順位になっているかを示します。

d.rank <- sort(rank(-abs(d))[spikedIndex])
t.rank <- sort(rank(-abs(tt$statistic))[spikedIndex])
mt.rank <- sort(rank(-abs(mtstat))[spikedIndex])
ranks <- cbind(mt.rank, d.rank, t.rank)
rownames(ranks) <- NULL
ranks

△実行結果

      mt.rank d.rank t.rank
 [1,]       1      1      1
 [2,]       2      2      3
 [3,]       3      3      5
 [4,]       6      4      8
 [5,]       7      5     13
 [6,]       8      6     17
 [7,]       9      9     19
 [8,]      10     11     27
 [9,]      11     12     28
[10,]      12     14     29
[11,]      16     16     45
[12,]      25     53     70
[13,]      28     86     71
[14,]      48    226     93
[15,]      77    331    390
[16,]     465    689    900

比較的、modarated t検定が良い結果でした。


関連するアノテーションの付加

統計検定量の付加

▼modareted t統計量をつかって、上位10遺伝子をとりだします。

tab <- topTable(ebayes, coef = 2, adjust = "fdr", n = 10)

▼tabの上位5つを表示します。

tab[1:5,]

△実行結果

           ID  logFC AveExpr      t  P.Value adj.P.Val    B
779   1708_at -7.061    7.95 -73.53 7.82e-17  9.87e-13 8.65
6252 36202_at  0.853    9.37   9.98 4.94e-07  3.12e-03 4.59
6362 36311_at  0.832    8.56   8.36 3.02e-06  1.27e-02 3.57
3284 33264_at  0.712    4.92   7.43 9.67e-06  2.71e-02 2.84
2674 32660_at  0.655    8.68   7.36 1.07e-05  2.71e-02 2.77

▼Affymetrix IDを取り出します。

genenames <- as.character(tab$ID)
PubMedアブストラクトの付加

▼eset が必要とするアノテーションパッケージを調べる。

annotation(eset)

△実行結果

[1] "hgu95a"

▼必要なパッケージを読み込み、関連した論文PubMed アブストラクトを取得する。

library("hgu95av2")
library("XML")
library("annotate")
absts <- pm.getabst(genenames, "hgu95av2")

△実行結果

Read 8255 items
Read 4181 items
Read 4277 items
Read 1329 items
Read 1157 items
Read 8379 items
Read 22358 items
Read 3816 items
Read 11971 items
Read 3653 items

アブストラクトの一部をみる。

absts[[1]][[4]]

△実行結果

An object of class 'pubMedAbst':
Title: Absence of excitotoxicity-induced apoptosis in
     the hippocampus of mice lacking the Jnk3 gene.
PMID: 9349820
Authors: DD Yang, CY Kuan, AJ Whitmarsh, M Rincón,
     TS Zheng, RJ Davis, P Rakic, RA Flavell
Journal: Nature
Date: Oct 1997

▼関連した論文タイトルを表示する。

titl <- sapply(absts[[2]], articleTitle)
strwrap(titl, simplify = FALSE)

△実行結果

[[1]]
[1] "Isolation and characterization of cDNA clones for an"
[2] "inhibitor protein of cAMP-dependent protein kinase." 

[[2]]
[1] "Inhibition of protein kinase-A by overexpression of"
[2] "the cloned human protein kinase inhibitor."         

[[3]]
[1] "Structure of a peptide inhibitor bound to the"
[2] "catalytic subunit of cyclic adenosine"        
[3] "monophosphate-dependent protein kinase."      

[[4]]
[1] "Tyrosine kinase catalyzed phosphorylation and"
[2] "inactivation of the inhibitor protein of the" 
[3] "cAMP-dependent protein kinase."               

[[5]]
[1] "Identification of an inhibitory region of the"      
[2] "heat-stable protein inhibitor of the cAMP-dependent"
[3] "protein kinase."                                    

[[6]]
[1] "Primary-structure requirements for inhibition by the"
[2] "heat-stable inhibitor of the cAMP-dependent protein" 
[3] "kinase."                                             

[[7]]
[1] "The expression and intracellular distribution of the"
[2] "heat-stable protein kinase inhibitor is cell cycle"  
[3] "regulated."                                          

[[8]]
[1] "Glutamic acid 203 of the cAMP-dependent protein" 
[2] "kinase catalytic subunit participates in the"    
[3] "inhibition by two isoforms of the protein kinase"
[4] "inhibitor."                                      

[[9]]
[1] "Evidence for the importance of hydrophobic residues"
[2] "in the interactions between the cAMP-dependent"     
[3] "protein kinase catalytic subunit and the protein"   
[4] "kinase inhibitors."                                 

[[10]]
[1] "HIV Rev uses a conserved cellular protein export"    
[2] "pathway for the nucleocytoplasmic transport of viral"
[3] "RNAs."                                               

[[11]]
[1] "PrKX is a novel catalytic subunit of the"      
[2] "cAMP-dependent protein kinase regulated by the"
[3] "regulatory subunit type I."                    

[[12]]
[1] "Protein inhibitor of neuronal nitric oxide synthase"
[2] "interacts with protein kinase A inhibitors."        

[[13]]
[1] "Generation and initial analysis of more than 15,000"
[2] "full-length human and mouse cDNA sequences."        

[[14]]
[1] "Mutants of protein kinase A that mimic the" 
[2] "ATP-binding site of protein kinase B (AKT)."

[[15]]
[1] "Development of human protein reference database as an"
[2] "initial platform for approaching systems biology in"  
[3] "humans."                                              

[[16]]
[1] "Human protein reference database as a discovery"
[2] "resource for proteomics."                       

[[17]]
[1] "Human protein reference database--2006 update."

▼R のなかで PubMed アブストラクトキーワード検索をおこなう。ここではキーワードが Protein もしくは protein です。

pro.res <- sapply(absts, function(x) pm.abstGrep("[Pp]rotein", x))
pro.res[[2]]

△実行結果

 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[10]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE

PubMed アブストラクトを閲覧するために HTML レポートを作成する。

pmAbst2HTML(absts[[2]], filename = "pm.html")

△実行結果

pm.html

f:id:nakao_mitsuteru:20070616231151p:image


レポートの生成

▼genenames 遺伝子リストに対応する LocusLink ID遺伝子名(Symbol)を取り出します。

ll <- getLL(genenames, "hgu95av2")
sym <- getSYMBOL(genenames, "hgu95av2")

▼つぎに、LocusLink ID (gene ID) を付与した HTML 形式のレポートを作成します。HTML 形式のレポートウェブブラウザで閲覧できるので便利です。

tab <- data.frame(sym, tab[, -1])
htmlpage(ll, filename = "report.html", 
                title = "HTML report",
           othernames = tab,
           table.head = c("Locus ID", colnames(tab)),
         table.center = TRUE)

△実行結果

report.html

f:id:nakao_mitsuteru:20070616231150p:image

▼report.html は tab の内容にリンクをつけたものです。

tab

△実行結果

            sym  logFC AveExpr      t  P.Value adj.P.Val     B
1708_at  MAPK10 -7.061    7.95 -73.53 7.82e-17  9.87e-13  8.65
36202_at   PKIA  0.853    9.37   9.98 4.94e-07  3.12e-03  4.59
36311_at  PDE1A  0.832    8.56   8.36 3.02e-06  1.27e-02  3.57
33264_at ENOSF1  0.712    4.92   7.43 9.67e-06  2.71e-02  2.84
32660_at   LBA1  0.655    8.68   7.36 1.07e-05  2.71e-02  2.77
38734_at    PLN  0.747    6.26   7.19 1.35e-05  2.83e-02  2.62
1024_at  CYP1A1  0.843    9.70   6.73 2.50e-05  4.40e-02  2.20
36085_at  GNAT1  0.645   12.19   6.65 2.79e-05  4.40e-02  2.12
33818_at    VCP  0.532   12.29   6.45 3.70e-05  5.19e-02  1.92
39058_at    ABR  0.609    7.53   6.28 4.77e-05  5.69e-02  1.74

KEGGやGOへのリンクを付与するには annaffy パッケージを使用します。

library("KEGG")
library("GO")
library("annaffy")
atab <- aafTableAnn(genenames, "hgu95av2", aaf.handler())
saveHTML(atab, file = "report2.html")

△実行結果

report2.html

f:id:nakao_mitsuteru:20070616231149p:image

report2.html は atab の内容をテーブル状に表示しています。atab は annaffy パッケージの aafTable クラスオブジェクトです。

▼カラム名を表示

colnames(atab['1708_at',])

△実行結果

 [1] "Probe"               "Symbol"             
 [3] "Description"         "Function"           
 [5] "Chromosome"          "Chromosome Location"
 [7] "GenBank"             "LocusLink"          
 [9] "Cytoband"            "UniGene"            
[11] "PubMed"              "Gene Ontology"      
[13] "Pathway" 

▼1708_at プローブ遺伝子シンボルを表示

atab['1708_at',]$Symbol

△実行結果

[[1]]
[1] "MAPK10"
attr(,"class")
[1] "aafSymbol"

attr(,"class")
[1] "aafList"

まとめ

CELファイルから、有意に発現変動する遺伝子を列挙する方法をみました。

発現変動のさまざまな基準が利用できます。

どの方法が良いかはデータ依存しています。

データパッケージがあると、さまざまな解析が簡単におこなえます。

PubMedアブストラクトタイトルはAffymetrix IDには簡単に対応づきます。

LocusLink や GO や KEGG などのメタデータを簡単に付与できます。

行列HTML のテーブルページにするメソッドがあります。




もくじにもどる

unjgoigpbcunjgoigpbc 2013/07/27 00:33 mxuftcjpdpoevdups, <a href="http://www.qfdljtypod.com/">lmcioritxd</a> , [url=http://www.gzzikngyus.com/]bdfczpsmri[/url], http://www.zftisrmfht.com/ lmcioritxd

nwnmrjzdklnwnmrjzdkl 2013/07/30 09:34 chgdqcjpdpoevdups, <a href="http://www.dtpjtpmcwz.com/">nwnjygkzzb</a> , [url=http://www.kdfgrbwcvr.com/]bypzwjsiji[/url], http://www.igffaiytve.com/ nwnjygkzzb

cdmzgyfdmrcdmzgyfdmr 2013/08/27 20:48 jjpbfcjpdpoevdups, <a href="http://www.sykbutykif.com/">ukuuydefwr</a> , [url=http://www.uhpsolyxko.com/]fzsrybsnxm[/url], http://www.ltbqizswiw.com/ ukuuydefwr

rzneokixsarzneokixsa 2013/12/19 06:07 ojooxcjpdpoevdups, <a href="http://www.mbhcuaxduh.com/">kcoetjdwiz</a>

lfenndqlkxlfenndqlkx 2013/12/21 07:18 qqkzkcjpdpoevdups, http://www.ezysqwcddb.com/ ahraopcipd

rcnnhmqvoercnnhmqvoe 2013/12/23 12:07 bxpfacjpdpoevdups, <a href="http://www.szrshwiogp.com/">cxeugadcxs</a> , [url=http://www.zzekrrgrud.com/]dxqnzipezs[/url], http://www.dydxvalriy.com/ cxeugadcxs

gbzpsabhzdgbzpsabhzd 2014/01/12 06:32 svcyrcjpdpoevdups, <a href="http://www.fmvxfxjkfd.com/">tnaexwiklw</a> , [url=http://www.gtxtbpmhhq.com/]rvsokjanyr[/url], http://www.hymjvzsvzt.com/ tnaexwiklw

JasonCoofeJasonCoofe 2017/01/25 04:28 печать журналов http://wkrolik.com.ua/products/otkrytki

ゲスト