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 のテーブルページにするメソッドがあります。




もくじにもどる

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

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

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

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

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

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

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

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

GuestAcirmGuestAcirm2017/09/06 10:27guest test post
<a href="http://googlee.te/">bbcode</a>
<a href="http://googlee.te/">html</a>
http://googlee.te/ simple

MaucardMaucard2018/02/06 07:04Propecia Dependency <a href=http://tadalaffbuy.com>cialis online</a> Cialis En Ligne Pas Chere Cialis Generika Gefahrlich

KevinButKevinBut2018/04/16 10:15erectile issues
<a href=http://erectiledysfunctionpillsntx.com/>erectile pills</a>
is erectile dysfunction psychological
<a href="http://erectiledysfunctionpillsntx.com/">best erectile pills</a>

MarkchadoMarkchado2018/04/17 00:49can erectile dysfunction be psychological
<a href=http://edpils-gg.com/>best erectile dysfunction medication</a>
erectile pillole
<a href="http://edpils-gg.com/">erectile pills canada</a>

eddrugsgeneric.comeddrugsgeneric.com2018/04/17 04:03erectile strength diet
<a href=http://eddrugsgeneric.com/#>new drugs for ed</a>
when does erectile dysfunction happen
<a href="http://eddrugsgeneric.com/#">buy erectile dysfunction pills online</a>

LdrgeDrevaLdrgeDreva2018/04/17 06:50erectile dysfunction in children <a href=http://dc-shuttle.com/>best erectile dysfunction medication</a>
erectile help <a href="http://dc-shuttle.com/">buy erectile dysfunction medications online</a> http://dc-shuttle.com/

KoldrSiKoldrSi2018/04/18 05:17erectile stamina
<a href=http://edpilsmystery.com/>ed drugs</a>
erectile shots
<a href="http://edpilsmystery.com/">erectile dysfunction medications</a>

EarnestexcibEarnestexcib2018/04/18 06:53erectile disorder dsm v
<a href="http://erectiledysfunctionlq.com/">buy erectile dysfunction pills online</a>
youth erectile dysfunction
<a href="http://erectiledysfunctionlq.com/">erectile pills over the counter</a>

MarkchadoMarkchado2018/04/19 06:38trading binary options with reliable broker reviews
<a href="http://options-binary-platform-try-best.pw/trade-secret-australia-return-policy-1260.php#">autotrader vans west midlands</a>
discount brokerage promotions
<a href=http://options-try-best-platform-binary.pw/binary-option-wikihow-3090.php>der handel binaroperationen 60220</a>

JosephLOnJosephLOn2018/04/19 18:16erectile herbal supplements
http://toperectiledysfunctionpills.com/ - best erectile pills
erectile pill sponsored by doctor oz
<a href="http://toperectiledysfunctionpills.com/">best drugs for ed</a>

HilopirtyHilopirty2018/04/19 19:54Viagra Generico Controindicazioni Best Buy Generic Cod Only Pyridium Plymouth <a href=http://brandciali.com>online pharmacy</a> Healtyman Generic Viagra For Sale In Usa

MarkchadoMarkchado2018/04/19 20:08levitra bahtera pratama ptd
<a href="http://mylevitraok.com">levitra prices</a>
soldado de levitra wikipedia
<a href=http://myvardenafilok.com#vardenafil+20mg>levitra 10 mg prezzo</a>

MarkchadoMarkchado2018/04/20 04:25potenzmittel levitra cost
<a href="http://myvardenafilok.com">levitra 10 mg prezzo</a>
cialis generic levitra viagra php
<a href=http://levitra-gg.com#levitra+coupon>levitra 20 mg</a>

MarkchadoMarkchado2018/04/20 08:06dosierungen levitra
<a href="http://myvardenafilok.com">levitra coupon</a>
levitra 20mg preis deutschland alles
<a href=http://myvardenafilok.com#generic+levitra>levitra online</a>

MarkchadoMarkchado2018/04/20 11:56viagra cialis levitra comparison
<a href="http://myvardenafilok.com">levitra 10 mg prezzo</a>
kas geriau cialis ar levitra vs cialis
<a href=http://myvardenafilok.com#vardenafil+20mg>vardenafil 20mg</a>

JugeDrevaJugeDreva2018/04/20 13:27viagra vs cialis200 cialis couponcialis 20 mg best pricecialis in usa
п»ї<a href="http://jvrimages.com/">cialis 20 mg</a>
generic for cialiscialis dosage
http://jvrimages.com/#cialis+official+site

MarkchadoMarkchado2018/04/20 15:38levitra administrare vitamina
<a href="http://mylevitraok.com">levitra 10 mg prezzo</a>
generico levitra orosolubile significato
<a href=http://myvardenafilok.com#vardenafil+20mg>buy levitra</a>

MarkchadoMarkchado2018/04/20 19:43levitra cavitation technologies
<a href="http://myvardenafilok.com">buy levitra</a>
buy levitra on line
<a href=http://levitra-gg.com#levitra+20mg>levitra coupon</a>

JorgeabuthJorgeabuth2018/04/20 22:58cialis for sale in europa prezzo di cialis in bulgaria

MarkchadoMarkchado2018/04/20 23:43buy cheap generic levitra online
<a href="http://myvardenafilok.com">levitra rezeptfrei deutschland</a>
levitra very cheap
<a href=http://levitra-gg.com#levitra+prices>levitra online</a>

YtereDrevaYtereDreva2018/04/21 03:55Buy Cheap Cialiscialis for daily usecialis 30 day samplecialis for sale <a href="http://foreigncircles.com/#Trazodone">buy Trazodone</a> cialis couponscialis cost <a href=http://foreigncircles.com/>Trazodone 50 mg best price</a> 200 cialis couponcialis free trial http://foreigncircles.com/#Trazodone

MarkchadoMarkchado2018/04/21 06:08levitra vs viagra alcohol side
<a href="http://levitra-gg.com">levitra 20mg</a>
levitra falcons score
<a href=http://mylevitraok.com#levitra+10+mg+prezzo>levitra online</a>

YtereDrevaYtereDreva2018/04/21 08:14cialis for daily usecialis savings cardcialis purchasingcialis canada <a href="http://foreigncircles.com/#Trazodone">buy Trazodone</a> Buy Cheap Cialis in usviagra vs cialis <a href=http://foreigncircles.com/>Trazodone 50 mg online</a> cialis 20 mggeneric for cialis http://foreigncircles.com/#Trazodone-50-mg

MarkchadoMarkchado2018/04/21 09:35viagra levitra ou cialis qual o melhor navegador
<a href="http://myvardenafilok.com">levitra 10 mg prezzo</a>
ketoconazole 20 mg generic levitra
<a href=http://levitra-gg.com#levitra+coupon>levitra 20mg</a>

MarkchadoMarkchado2018/04/21 12:56como se dice levitra en ingles
<a href="http://levitra-gg.com">levitra 20 mg bayer prezzo</a>
viagra cialis ou levitra qual o melhor investimento
<a href=http://myvardenafilok.com#generic+levitra>levitra online</a>

MarkchadoMarkchado2018/04/21 16:18levitra siglo xviii literatura
<a href="http://mylevitraok.com">levitra online</a>
baixaki playback levitra flavio canto
<a href=http://mylevitraok.com#levitra+20+mg>vardenafil 20mg</a>

YtereDrevaYtereDreva2018/04/22 05:58cialis side effectscialis costcialis generic availabilityBuy Cheap Cialis <a href="https://cheapcialis.wixsite.com/cialis">cheap cialis</a>

tadalafil tabletstadalafil tablets2018/04/26 06:26
Truly a lot of terrific material!

JamesfefJamesfef2018/04/26 09:51cheap viagra cialis online
<a href="http://cialispaxl.com/">generic cialis</a>
viagra and cialis price comparison
<a href="http://cialispaxl.com/">cialis</a>
cialis dosing options
http://cialispaxl.com/
cialis dose maximale

RobertPrickRobertPrick2018/04/26 10:07viagra girl commercial blue dress
<a href="http://viagrapipls.com/">generic viagra</a>
compare price of viagra cialis and levitra
<a href="http://viagrapipls.com/">viagra generic</a>
como se debe tomar un viagra
http://viagrapipls.com/
date for generic viagra

LouisPefLouisPef2018/04/26 10:33cialis shop usa
<a href="http://cialisgrudj.com/">cialis generic</a>
cialis professional 20 mg
<a href="http://cialisgrudj.com/">generic cialis online</a>
long term side effects of cialis
http://cialisgrudj.com/
side effects of viagra and cialis

TrevorVerTrevorVer2018/04/26 17:06cialis vs viagra user reviews
<a href="http://cialispaxl.com/">buy cialis</a>
buy generic cialis online in canada
<a href="http://cialispaxl.com/">buy cialis online</a>
cialis generic 2018 expiration patent
http://cialispaxl.com/
main ingredient in cialis

BrucejaismBrucejaism2018/04/26 17:36fiat viagra commercial youtube
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra dosage instructions
<a href="http://viagrapipls.com/">generic viagra</a>
cost comparison viagra levitra cialis
http://viagrapipls.com/
compare cialis and viagra side effects

MarvinMupMarvinMup2018/04/26 23:03buy cheap generic viagra co uk kamagra oral jelly 100mg
<a href="http://viagrapipls.com/">generic viagra</a>
buy cialis viagra online
<a href="http://viagrapipls.com/">generic viagra online</a>
viagra price vs cialis
http://viagrapipls.com/
cost of generic viagra in india

PeterTitPeterTit2018/04/27 01:46cialis side effects loss of vision
<a href="http://cialisgrudj.com/">cialis generic</a>
preço cialis generico brasil
<a href="http://cialisgrudj.com/">generic cialis</a>
cialis dosagem diaria
http://cialisgrudj.com/
cvs pharmacy cialis cost

FrancisfireeFrancisfiree2018/04/27 03:26cialis tadalafil 20mg reviews
<a href="http://cialispaxl.com/">cialis generic</a>
cialis side effects dangers oral
<a href="http://cialispaxl.com/">buy cialis online</a>
cialis 5mg price walgreens
http://cialispaxl.com/
acquisto cialis generico in italia

MichaelMedayMichaelMeday2018/04/27 04:42viagra commercial actress name black
<a href="http://viagrapipls.com/">buy viagra online</a>
generic viagra prices at walgreens
<a href="http://viagrapipls.com/">viagra</a>
price comparison viagra cialis levitra
http://viagrapipls.com/
2016 viagra commercial actress

WilliamchodoWilliamchodo2018/04/27 08:31cialis when generic in us
<a href="http://cialispaxl.com/">buy generic cialis online</a>
cialis generico precio espaГ±a
<a href="http://cialispaxl.com/">generic cialis</a>
effectiveness of cialis vs viagra vs levitra reviews
http://cialispaxl.com/
levitra viagra cialis cost comparison

LucianohenLucianohen2018/04/27 09:51lady in viagra commercial 2016
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra effects on the liver
<a href="http://viagrapipls.com/">buy viagra</a>
viagra maximum daily dosage
http://viagrapipls.com/
viagra going generic date wikipedia

DwightcaxDwightcax2018/04/27 13:31cialis coupon 5mg cvs
<a href="http://cialispaxl.com/">cialis</a>
cialis going generic date
<a href="http://cialispaxl.com/">buy generic cialis</a>
generic cialis release date
http://cialispaxl.com/
cialis retail price at cvs

GregorySonGregorySon2018/04/27 15:30brand viagra online
<a href="http://viagrapipls.com/">viagra generic</a>
viagra professional 100mg ft myers florida
<a href="http://viagrapipls.com/">buy viagra</a>
cialis vs viagra cost comparison
http://viagrapipls.com/
precisa ter receita para comprar viagra

BillycabBillycab2018/04/27 16:44best price cialis generic
<a href="http://cialisgrudj.com/">generic cialis</a>
dosage for 20mg cialis
<a href="http://cialisgrudj.com/">generic cialis</a>
cialis price walmart pharmacy
http://cialisgrudj.com/
cialis professional dosage

RobertTiedeRobertTiede2018/04/27 19:24buy cialis in costa rica
<a href="http://cialispaxl.com/">cialis</a>
generic cialis safety
<a href="http://cialispaxl.com/">buy generic cialis online</a>
cialis generico 5 mg preço
http://cialispaxl.com/
buy cheap viagra and cialis

MorrisAssepMorrisAssep2018/04/27 21:35generic viagra from walmart
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra dosage chart
<a href="http://viagrapipls.com/">generic viagra</a>
female viagra pills
http://viagrapipls.com/
viagra super active 150mg for sale

AndrewImaniAndrewImani2018/04/28 00:29cialis generico in farmacia italia
<a href="http://cialispaxl.com/">buy cialis online</a>
onde comprar cialis generico no brasil
<a href="http://cialispaxl.com/">buy generic cialis</a>
cialis 5mg dosage
http://cialispaxl.com/
cialis price comparison walmart

CharlesTogCharlesTog2018/04/28 02:50viagra super active sildenafil citrate
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra vs cialis price in india
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra commercial football actors
http://viagrapipls.com/
viagra professional 100mg pills

JerryVaddyJerryVaddy2018/04/28 05:16viagra and cialis price comparison
<a href="http://cialispaxl.com/">buy generic cialis online</a>
viagra cialis levitra generici
<a href="http://cialispaxl.com/">buy generic cialis online</a>
cialis for daily use side effects
http://cialispaxl.com/
best generic cialis prices

KeithFeateKeithFeate2018/04/28 08:47levitra cialis or viagra which is better
<a href="http://cialisgrudj.com/">buy generic cialis</a>
precio cialis generico en farmacias
<a href="http://cialisgrudj.com/">buy generic cialis</a>
buy generic viagra cialis levitra
http://cialisgrudj.com/
physician samples of cialis

JamesTwityJamesTwity2018/04/28 10:25brand cialis 5mg
<a href="http://cialispaxl.com/">buy cialis</a>
cialis 20 mg vs 5mg daily
<a href="http://cialispaxl.com/">buy generic cialis</a>
generic cialis available in us
http://cialispaxl.com/
cialis side effects vs viagra

RichardBrastRichardBrast2018/04/28 13:35youtube viagra commercial 2015 date night
<a href="http://viagrapipls.com/">viagra generic</a>
best price viagra cialis
<a href="http://viagrapipls.com/">buy generic viagra online</a>
viagra dose riddim
http://viagrapipls.com/
female viagra pill review

SidneyAxionSidneyAxion2018/04/28 15:38cialis soft review
<a href="http://cialispaxl.com/">buy generic cialis online</a>
levitra vs cialis cost
<a href="http://cialispaxl.com/">buy cialis online</a>
cialis vs viagra vs levitra reddit
http://cialispaxl.com/
generic cialis available in us

PhillipteellPhillipteell2018/04/28 19:11name of woman in viagra commercial blue dress
<a href="http://viagrapipls.com/">buy viagra</a>
cialis dose recommendations vs viagra
<a href="http://viagrapipls.com/">buy generic viagra</a>
viagra tablet price in india
http://viagrapipls.com/
dosage for viagra in dogs

MarkchadoMarkchado2018/04/28 20:09erectile rehabilitation therapy
<a href=http://edpils-gg.com/>erectile dysfunction medications</a>
buy erectile dysfunction pills online india
<a href="http://edpils-gg.com/">erectile pills without a doctor prescription</a>

RogerHasseRogerHasse2018/04/28 20:51viagra cialis cost comparisons
<a href="http://cialispaxl.com/">buy cialis</a>
cialis with dapoxetine review
<a href="http://cialispaxl.com/">generic cialis online</a>
cialis side effects vision impairment
http://cialispaxl.com/
cialis super force dapoxetine review

DonaldWepDonaldWep2018/04/29 00:46when will viagra prices drop
<a href="http://viagrapipls.com/">buy viagra online</a>
what is generic viagra safe
<a href="http://viagrapipls.com/">buy generic viagra</a>
female viagra spoof commercial
http://viagrapipls.com/
viagra professional vs 200 mg viagra

JimmyCahJimmyCah2018/04/29 02:08side effects of cialis for daily use
<a href="http://cialispaxl.com/">cialis generic</a>
viagra vs cialis dose
<a href="http://cialispaxl.com/">buy generic cialis</a>
generic cialis capsules
http://cialispaxl.com/
cialis side effects dangers or levitra vs cialis vs

BrendanraxBrendanrax2018/04/29 07:08taking extra cialis 5 mg doses
<a href="http://cialispaxl.com/">cialis generic</a>
once daily use generic cialis
<a href="http://cialispaxl.com/">generic cialis online</a>
cialis vs viagra
http://cialispaxl.com/
cialis soft 20mg

MatthewGlyncMatthewGlync2018/04/29 11:26generic viagra 100mg sildenafil meds express
<a href="http://viagrapipls.com/">generic viagra</a>
viagra generic walmart cost
<a href="http://viagrapipls.com/">generic viagra</a>
prices for viagra at walmart
http://viagrapipls.com/
generic viagra soft tabs 100mg

StevencarStevencar2018/04/29 12:07cialis diario generico preço
<a href="http://cialispaxl.com/">cialis</a>
cialis 5 mg coupon lilly
<a href="http://cialispaxl.com/">buy cialis</a>
cialis super active plus
http://cialispaxl.com/
dosage strengths of cialis

HectorlabHectorlab2018/04/29 15:06cialis dosage 5mg or 10mg
<a href="http://cialisgrudj.com/">generic cialis online</a>
comprar cialis generico no brasil
<a href="http://cialisgrudj.com/">cialis</a>
cialis commercial black woman
http://cialisgrudj.com/
levitra cialis viagra comparison

CollindurneCollindurne2018/04/29 17:07generic viagra goodrx
<a href="http://viagrapipls.com/">buy viagra online</a>
cialis dose vs viagra dose
<a href="http://viagrapipls.com/">buy generic viagra online</a>
female viagra pills name in india
http://viagrapipls.com/
pfizer brand viagra 100mg sales

NelsonzixNelsonzix2018/04/29 22:36como hacer viagra natural para mujeres
<a href="http://viagrapipls.com/">buy viagra online</a>
maximum daily dosage of viagra
<a href="http://viagrapipls.com/">buy viagra</a>
viagra generico preço rj
http://viagrapipls.com/
cialis generic levitra viagra

KennethnuabeKennethnuabe2018/04/30 03:51cialis price versus viagra side effects
<a href="http://cialispaxl.com/">buy cialis online</a>
cialis 5mg best price
<a href="http://cialispaxl.com/">buy cialis</a>
price comparison viagra and cialis
http://cialispaxl.com/
cialis daily use side effects

WilliamPugWilliamPug2018/04/30 06:14cialis soft gel caps
<a href="http://cialisgrudj.com/">generic cialis</a>
price of cialis and viagra
<a href="http://cialisgrudj.com/">cialis generic</a>
cialis generic availability
http://cialisgrudj.com/
cialis dosage information

JosephBrataJosephBrata2018/04/30 08:55cialis dosage info
<a href="http://cialispaxl.com/">buy cialis</a>
cialis price at walmart pharmacy
<a href="http://cialispaxl.com/">cialis</a>
buy generic viagra and cialis online canada
http://cialispaxl.com/
cialis coupon rite aid

DavidadubsDavidadubs2018/04/30 15:03viagra precisa de receita controlada
<a href="http://viagrapipls.com/">generic viagra online</a>
viagra commercial actress asian
<a href="http://viagrapipls.com/">buy generic viagra online</a>
cost of viagra cialis and levitra
http://viagrapipls.com/
effects of viagra on females

ReteDrevaReteDreva2018/04/30 19:03Buy Cheap Cialis usacialis purchasinggeneric cialis <a href="http://chineseteacrafts.com/">Buy Cialis Without A Doctor Prescription</a>
Buy Cheap Cialisbuy cialis online <a href="http://chineseteacrafts.com/">Buy Generic Cialis Online Without Prescription</a>

VictorJagVictorJag2018/04/30 19:19cialis dosage and timing
<a href="http://cialispaxl.com/">cialis</a>
cialis 5 mg daily side effects
<a href="http://cialispaxl.com/">cialis generic</a>
generic viagra cialis levitra
http://cialispaxl.com/
viagra cialis cost

RobertBeelpRobertBeelp2018/04/30 20:40cialis vs viagra vs levitra dosage reviews
<a href="http://viagrapipls.com/">buy generic viagra</a>
viagra levitra cost comparison
<a href="http://viagrapipls.com/">buy viagra online</a>
cheap viagra in usa
http://viagrapipls.com/
generic viagra super active reviews

GeorgeFedGeorgeFed2018/04/30 21:14generic cialis coming out
<a href="http://cialisgrudj.com/">generic cialis online</a>
effectiveness of cialis vs viagra vs levitra alcohol
<a href="http://cialisgrudj.com/">buy generic cialis online</a>
walmart pharmacy cialis price check
http://cialisgrudj.com/
buy cialis in usa online

RandallStupsRandallStups2018/05/01 00:33cialis coupon rite aid
<a href="http://cialispaxl.com/">cialis generic</a>
cialis dosage strengths and weaknesses
<a href="http://cialispaxl.com/">buy generic cialis</a>
cialis side effects red eyes
http://cialispaxl.com/
cialis printable coupon 2017

RogerserRogerser2018/05/01 05:31cialis dosage 20mg
<a href="http://cialispaxl.com/">buy generic cialis</a>
cialis coupons for walgreens pharmacy
<a href="http://cialispaxl.com/">generic cialis</a>
cialis prescription price australia
http://cialispaxl.com/
cialis soft tabs uk

DannyIdeotDannyIdeot2018/05/01 10:53cialis generic release in us
<a href="http://cialisonlinq.com/">buy generic cialis</a>
cialis super active 40 mg
<a href="http://cialisonlinq.com/">cialis generic</a>
cialis tablets 20mg for sale
http://cialisonlinq.com/
is generic cialis available in europe

DavidnewDavidnew2018/05/01 11:20viagra retail price cvs
<a href="http://viagrapipls.com/">viagra generic</a>
viagra stendra cialis levitra side effects
<a href="http://viagrapipls.com/">viagra generic</a>
generic viagra 100mg sildenafil prescription
http://viagrapipls.com/
viagra vs. cialis vs. levitra which is best

HenrySopHenrySop2018/05/01 15:40cialis commercial song 2012
<a href="http://cialisonlinq.com/">cialis generic</a>
daily 5mg cialis price
<a href="http://cialisonlinq.com/">generic cialis online</a>
viagra cialis comparison
http://cialisonlinq.com/
generics for viagra and cialis at costco

PrestonSlopePrestonSlope2018/05/02 04:23best generic cialis pills price
<a href="http://cialisonlinq.com/">buy cialis</a>
cialis 5 mg generic best price india
<a href="http://cialisonlinq.com/">generic cialis</a>
cialis and viagra generic
http://cialisonlinq.com/
cialis when generic in us

DavidNarDavidNar2018/05/02 05:05viagra commercial bathroom
<a href="http://viagrapipls.com/">viagra generic</a>
fiat viagra commercial youtube
<a href="http://viagrapipls.com/">buy viagra</a>
viagra commercial brunette woman
http://viagrapipls.com/
brand viagra australia

DannyicokeDannyicoke2018/05/02 07:18buy generic cialis online reviews
<a href="http://cialisonlinq.com/">buy generic cialis online</a>
viagra price vs cialis
<a href="http://cialisonlinq.com/">generic cialis</a>
walgreens price for cialis 20mg
http://cialisonlinq.com/
walmart pharmacy cialis 5 mg cost

StevenwetStevenwet2018/05/02 22:30dosage for viagra cialis levitra
<a href="http://cialisonlinq.com/">cialis generic</a>
cialis 5 mg daily side effects
<a href="http://cialisonlinq.com/">buy cialis online</a>
cialis cena apoteka
http://cialisonlinq.com/
cost of generic cialis

MichaelbibMichaelbib2018/05/02 23:14viagra professional 100 mg pills kiwi
<a href="http://viagrapipls.com/">buy viagra online</a>
levitra vs viagra cost
<a href="http://viagrapipls.com/">buy viagra online</a>
como fazer viagra caseiro masculino
http://viagrapipls.com/
viagra casero para hombres instantaneo

WillieberWillieber2018/05/02 23:15generic cialis available yeti
<a href="http://cialisonlinq.com/">cialis generic</a>
is generic cialis safe
<a href="http://cialisonlinq.com/">buy cialis</a>
cialis side effects dangers or levitra
http://cialisonlinq.com/
cialis commercial british lady

RoyceTibRoyceTib2018/05/03 14:28cialis commercial 2013 actors
<a href="http://cialisonlinq.com/">cialis generic</a>
cialis super active plus rezeptfrei
<a href="http://cialisonlinq.com/">cialis generic</a>
professional samples of cialis
http://cialisonlinq.com/
20mg cialis vs viagra 100mg

RodolfoWedRodolfoWed2018/05/03 15:52cheap cialis and viagra
<a href="http://cialisonlinq.com/">buy cialis online</a>
cialis generico colombia
<a href="http://cialisonlinq.com/">generic cialis</a>
buy cialis online usa
http://cialisonlinq.com/
brand name cialis 20mg

EugeneRoogsEugeneRoogs2018/05/03 16:36viagra and cialis comparison
<a href="http://viagrapipls.com/">generic viagra</a>
cost per pill for generic viagra at walmart
<a href="http://viagrapipls.com/">generic viagra online</a>
viagra without a doctor prescription reddit
http://viagrapipls.com/
prices of cialis and viagra

JerrychideJerrychide2018/05/04 05:10cialis soft reviews
<a href="http://cialisonlinq.com/">buy generic cialis</a>
viagra and cialis best price
<a href="http://cialisonlinq.com/">generic cialis</a>
walmart pharmacy cialis price checklist
http://cialisonlinq.com/
cialis 30 lu tablet 20 mg

AndrewBusAndrewBus2018/05/04 08:30cost of cialis 5mg for daily use
<a href="http://cialisonlinq.com/">buy generic cialis</a>
how long does it take a 5 mg cialis to work
<a href="http://cialisonlinq.com/">cialis</a>
cialis vs viagra for bph
http://cialisonlinq.com/
cialis generic name

RobertoRakRobertoRak2018/05/04 09:13viagra maximum safe dosage
<a href="http://viagrapipls.com/">generic viagra</a>
viagra dosage strengths
<a href="http://viagrapipls.com/">buy generic viagra online</a>
female viagra pills price
http://viagrapipls.com/
generic viagra india 100mg trustworthy site

JustinCheftJustinCheft2018/05/04 20:06generic cialis lowest price
<a href="http://cialisonlinq.com/">generic cialis online</a>
order generic cialis online canada
<a href="http://cialisonlinq.com/">buy cialis</a>
cialis generico opinion
http://cialisonlinq.com/
cialis super active reviews

JosephAleceJosephAlece2018/05/05 01:59viagra and cialis dosage
<a href="http://cialisonlinq.com/">buy cialis online</a>
cialis 20mg tablets side effects
<a href="http://cialisonlinq.com/">cialis generic</a>
cialis generic 20mg price shop
http://cialisonlinq.com/
cialis 5 mg generico prezzo

HarrySeimeHarrySeime2018/05/05 02:38buy real viagra levitra cialis online edrugstore.md
<a href="http://viagrapipls.com/">buy generic viagra online</a>
generic viagra in united states
<a href="http://viagrapipls.com/">buy generic viagra</a>
viagra soft pill
http://viagrapipls.com/
viagra soft tabs 100mg review

MichaelDamMichaelDam2018/05/05 19:01generic cialis vs brand cialis reviews
<a href="http://cialisonlinq.com/">generic cialis</a>
lowest cost generic cialis
<a href="http://cialisonlinq.com/">buy cialis</a>
cialis generic india
http://cialisonlinq.com/
cialis 100mg dosage information

DonaldinfedDonaldinfed2018/05/05 19:39viagra natural chino
<a href="http://viagrapipls.com/">viagra</a>
viagra for women pink
<a href="http://viagrapipls.com/">buy generic viagra online</a>
viagra preço drogaria sp
http://viagrapipls.com/
cialis dosage 20mg vs 100mg viagra

MarkchadoMarkchado2018/05/05 23:41koncepcja produktu levitra coupon
<a href="http://myvardenafilok.com">levitra coupon</a>
levitra 20mg directions
<a href=http://mylevitraok.com#levitra+rezeptfrei+deutschland>levitra 20mg</a>

RusselRogRusselRog2018/05/06 02:25cialis or viagra price
<a href="http://cialisonlinq.com/">cialis generic</a>
levitra side effects levitra vs cialis
<a href="http://cialisonlinq.com/">buy generic cialis online</a>
generic cialis daily 5 mg
http://cialisonlinq.com/
cialis side effects long term

MarkchadoMarkchado2018/05/06 03:07nussberg levitra
<a href="http://myvardenafilok.com">levitra coupon</a>
andy levitra contract extension
<a href=http://myvardenafilok.com#buy+levitra>levitra rezeptfrei deutschland</a>

MarkchadoMarkchado2018/05/06 06:30viagra cialis levitra canadian pharmacy
<a href="http://mylevitraok.com">levitra 20mg</a>
s kym si rozumie levitra
<a href=http://levitra-gg.com#levitra+20+mg>levitra 20 mg</a>

MarkchadoMarkchado2018/05/06 09:43levitra usage information
<a href="http://levitra-gg.com">levitra 20 mg</a>
levitra bayer indonesia
<a href=http://levitra-gg.com#levitra+coupon>levitra prices</a>

MarkchadoMarkchado2018/05/06 12:53levitra reviews medication
<a href="http://mylevitraok.com">levitra prices</a>
cialas vs levitra
<a href=http://mylevitraok.com#levitra+generic>levitra 20 mg bayer prezzo</a>

RichardTycleRichardTycle2018/05/06 13:03compare viagra cialis levitra side effects
<a href="http://viagrapipls.com/">generic viagra online</a>
viagra commercial
<a href="http://viagrapipls.com/">buy viagra</a>
150 mg viagra dosage recommendations
http://viagrapipls.com/
effects of viagra on healthy men

JoddyabuthJoddyabuth2018/05/06 15:58<a href="http://bikecommunications.com/#">http://bikecommunications.com/</a> rx cialis para comprar order generic cialis online does cialis cause gout cialis rezeptfrei sterreich how do cialis pills work cialis taglich <a href=http://bikecommunications.com/#Coupon>cialis 30 day trial coupon</a>

MarkchadoMarkchado2018/05/06 15:59patentschutz levitra
<a href="http://levitra-gg.com">buy levitra</a>
levitra orodispersible contre indication total
<a href=http://myvardenafilok.com#levitra+rezeptfrei+deutschland>levitra 20 mg bayer prezzo</a>

BrendanBipBrendanBip2018/05/06 17:44buy real viagra levitra cialis online edrugstore.md
<a href="http://cialisonlinq.com/">buy generic cialis online</a>
generic names for cialis
<a href="http://cialisonlinq.com/">generic cialis online</a>
price of viagra and cialis per pill
http://cialisonlinq.com/
cialis generic

MarkchadoMarkchado2018/05/06 19:06levitra duodenum
<a href="http://levitra-gg.com">levitra rezeptfrei deutschland</a>
united states levitra commercial babe
<a href=http://myvardenafilok.com#generic+levitra>generic levitra</a>

MarkchadoMarkchado2018/05/06 22:17aty an tany levitra for women
<a href="http://levitra-gg.com">levitra online</a>
levitra and aleve
<a href=http://mylevitraok.com#generic+levitra>levitra 10 mg prezzo</a>

MarkchadoMarkchado2018/05/07 01:36levitra e levitinha cifras y
<a href="http://myvardenafilok.com">levitra 10 mg prezzo</a>
levitra 10 mg costochondritis
<a href=http://myvardenafilok.com#levitra+online>buy levitra</a>

MarkchadoMarkchado2018/05/07 04:59soy guarura de levitra oscar chavez parodias
<a href="http://mylevitraok.com">vardenafil 20mg</a>
levitra antiguo testamento catolica
<a href=http://levitra-gg.com#levitra+online>levitra 20mg</a>

RalpharcagRalpharcag2018/05/07 05:38cialis dosagem
<a href="http://cialisonlinq.com/">buy generic cialis online</a>
generic 5mg cialis for daily use
<a href="http://cialisonlinq.com/">generic cialis online</a>
cheap viagra cialis levitra
http://cialisonlinq.com/
cialis generic availability 2017

LamarCusLamarCus2018/05/07 06:10price viagra vs cialis
<a href="http://viagrapipls.com/">buy viagra online</a>
viagra levitra cialis dosage
<a href="http://viagrapipls.com/">buy generic viagra</a>
cialis viagra dosage comparison
http://viagrapipls.com/
viagra natural barrio chino

MarkchadoMarkchado2018/05/07 08:27levitra y sacerdote alli pasaron por
<a href="http://levitra-gg.com">levitra 10 mg prezzo</a>
einnahme cialis vs levitra
<a href=http://levitra-gg.com#levitra+generic>levitra coupon</a>

MichaelenadaMichaelenada2018/05/07 08:33black woman in cialis commercial
<a href="http://cialisonlinq.com/">cialis</a>
levitra vs sildenafil vs cialis
<a href="http://cialisonlinq.com/">generic cialis</a>
cialis generico preço no brasil
http://cialisonlinq.com/
cialis professional reviews

MarkchadoMarkchado2018/05/07 11:55effekten av cialis vs levitra
<a href="http://levitra-gg.com">levitra 20 mg bayer prezzo</a>
traje levitra 1900722
<a href=http://levitra-gg.com#levitra+20+mg>levitra 10 mg prezzo</a>

MarkchadoMarkchado2018/05/07 15:11viagra levitra cialis offers
<a href="http://levitra-gg.com">levitra prices</a>
levitra flavio rar files
<a href=http://myvardenafilok.com#levitra+generic>levitra prices</a>

DavidPiextDavidPiext2018/05/07 22:48levitra vs sildenafil vs cialis
<a href="http://cialisonlinq.com/">cialis</a>
viagra vs cialis effectiveness
<a href="http://cialisonlinq.com/">buy generic cialis</a>
cost of viagra cialis and levitra
http://cialisonlinq.com/
cialis vs viagra for bph

AnthonyLonAnthonyLon2018/05/07 23:38cialis dosage options
<a href="http://cialisonlinq.com/">cialis generic</a>
buy cheap cialis in canada
<a href="http://cialisonlinq.com/">buy cialis</a>
generic cialis release date
http://cialisonlinq.com/
buy cialis cheap prices fast delivery

HectorLizHectorLiz2018/05/08 15:32cialis generico preço brasil
<a href="http://cialisonlinq.com/">cialis generic</a>
generic cialis tadalafil 20mg best prices
<a href="http://cialisonlinq.com/">generic cialis online</a>
generic cialis brands india
http://cialisonlinq.com/
cytotechnology specialist

SteveJaxSteveJax2018/05/09 05:29walmart price for cialis 20mg
<a href="http://cialisonlinq.com/">cialis</a>
buy generic cialis online india
<a href="http://cialisonlinq.com/">generic cialis online</a>
cialis generic vs brand name
http://cialisonlinq.com/
viagra cialis cost

RichardSnaskRichardSnask2018/05/09 09:19viagra super active review
<a href="http://amsboatyard.com/">generic viagra</a>
generic viagra online united states
<a href="http://amsboatyard.com/">generic viagra online</a>
viagra commercial actress name football
http://amsboatyard.com/
viagra prices

EdwardExegeEdwardExege2018/05/09 12:26generic cialis in the united states
<a href="http://cialisonlinq.com/">buy generic cialis online</a>
generic cialis is it safe
<a href="http://cialisonlinq.com/">buy generic cialis</a>
cialis 5mg tablets bedwetting
http://cialisonlinq.com/
cialis coupons printable 2017

EdgarFefEdgarFef2018/05/09 19:57cialis 5 mg 30 tablets cost
<a href="http://airvietnamairline.com/">generic cialis</a>
cialis side effects pain in leg
<a href="http://airvietnamairline.com/">generic cialis online</a>
cialis 5 mg generico preço
http://airvietnamairline.com/
cialis shop usa

MichaelSubMichaelSub2018/05/09 23:41kamagra gold teeth
<a href="http://kamagraonl.com/">buy kamagra 100mg</a>
kamagra oral jelly 100mg sildenafil citrate
<a href="http://kamagraonl.com/">buy kamagra</a>
kamagra vs viagra
http://kamagraonl.com/
kamagra 100mg oral jelly review

JamesbitJamesbit2018/05/10 02:19pfizer viagra coupons walgreens
<a href="http://amsboatyard.com/">buy viagra</a>
name of generic viagra
<a href="http://amsboatyard.com/">buy generic viagra online</a>
viagra price drop to boost sales
http://amsboatyard.com/
viagra natural receta

JosephfepJosephfep2018/05/10 11:11date generic cialis is available
<a href="http://airvietnamairline.com/">buy cialis online</a>
cialis 20 mg 30 tablets
<a href="http://airvietnamairline.com/">cialis</a>
preço cialis generico brasil
http://airvietnamairline.com/
best prices for cialis 5mg

MarionloombMarionloomb2018/05/10 16:40kamagra kopen apotheek
<a href="http://kamagraonl.com/">kamagra 100 mg oral jelly</a>
kamagra uk
<a href="http://kamagraonl.com/">buy kamagra online</a>
kamagra chew tablets - 100 mg
http://kamagraonl.com/
kamagra oral jelly amazon nederland texas

MartinGafMartinGaf2018/05/10 19:12viagra without a doctor prescription india
<a href="http://amsboatyard.com/">buy viagra online</a>
generic viagra vs cialis
<a href="http://amsboatyard.com/">generic viagra online</a>
herbal viagra side effects
http://amsboatyard.com/
levitra or viagra which is better

JefferyesopyJefferyesopy2018/05/11 02:07cialis coupon cvs
<a href="http://airvietnamairline.com/">cialis</a>
cialis available generic in 2017
<a href="http://airvietnamairline.com/">generic cialis online</a>
cialis 20 mg price walmart
http://airvietnamairline.com/
cialis side effects reviews

JefferyLahJefferyLah2018/05/11 05:50the kamagra store scam
kamagra 100mg oral jelly suppliers indianapolis indiana
kamagra gold from ajanta pharma
buy kamagra oral jelly usa

JefferyLahJefferyLah2018/05/11 16:06kamagra gel opinie forum
kamagra 100mg chewable tablets
ajanta kamagra 100 chewable
kamagra oral jelly wirkungsdauer

JefferyLahJefferyLah2018/05/11 22:19kamagra forum romania
kamagra 100mg tablets uk
ajanta kamagra 100 chewable tablet
kamagra forum doctissimo

JefferyLahJefferyLah2018/05/12 04:21cost of kamagra
the kamagra store hoax
kamagra winkel utrecht
kamagra online bestellen nederland

JefferyLahJefferyLah2018/05/12 04:58kamagra oral jelly kaufen deutschland paypal
kamagra forum srpski
kamagra cost
kamagra kopen afhalen amsterdam

SteevchadoSteevchado2018/05/14 00:46levitra everyday
<a href="http://myvardenafilok.com">coupon levitra</a>
archive info levitra personal php remember
<a href=http://mylevitraok.com#levitra+coupon>cheap levitra 10 mg</a>
levitra professional 100mg

SteevchadoSteevchado2018/05/14 03:04levitra side effects duration
<a href="http://myvardenafilok.com">10 mg levitra</a>
janet aponte levitra guerrero viejo
<a href=http://myvardenafilok.com#20+mg+levitra>levitra generic</a>
levitra vagisaling

SteevchadoSteevchado2018/05/14 05:15levitra joke
<a href="http://myvardenafilok.com">20 mg levitra</a>
acquistare levitra generico online thesaurus
<a href=http://mylevitraok.com#vardenafil+buy>coupon levitra</a>
levitra sample coupon

SteevchadoSteevchado2018/05/14 07:35levitra spanish translation
<a href="http://mylevitraok.com">levitra generic</a>
000002b4.htm commercialization discussion levitra sfp.forprod.vt.edu
<a href=http://mylevitraok.com#levitra+online>levitra 20 mg</a>
levitra 20 mg instructions on how to play

SteevchadoSteevchado2018/05/14 09:56avanafil vs levitra coupons
<a href="http://myvardenafilok.com">levitra 10mg prezzo</a>
levitra antsanitia
<a href=http://mylevitraok.com#vardenafil+prices>buy 20mg levitra</a>
levitra pills look like

SteevchadoSteevchado2018/05/14 12:11levitra dosage strengths of requip
<a href="http://mylevitraok.com">levitra generic</a>
buy levitra online 7 news
<a href=http://myvardenafilok.com#buy+vardenafil+online>levitra 10 mg prezzo</a>
can you take percocet with levitra

SteevchadoSteevchado2018/05/14 14:18vardenafil levitra tadalafil cialis canada
<a href="http://mylevitraok.com">buy 10 mg levitra</a>
aty an tany levitra 20mg
<a href=http://mylevitraok.com#levitra+coupon>vardenafil prices</a>
levitra side effects youtube

SteevchadoSteevchado2018/05/14 16:18levitra flavio parenti
<a href="http://mylevitraok.com">cheap 10mg levitra</a>
levitra babyfoon camera repair
<a href=http://mylevitraok.com#levitra+buy>levitra prices</a>
stan mc levitra

SteevchadoSteevchado2018/05/14 18:22levitra price in pakistan iphone
<a href="http://mylevitraok.com">buy 10 mg levitra</a>
levitra flavio vine
<a href=http://mylevitraok.com#buy+levitra+10mg>cheap levitra 20 mg</a>
buy cheap levitra link online jixx de