＊ここに掲載するのは、エジンバラ大学のLawrence Mitchell博士によるHECToRレポート「A parallel implementation of the Rank Product method for R, Lawrence Mitchell, April 30, 2011」を要約したものです。
概要
RライブラリSPRINTの一部であるrank product法の並列化を行いました。これはマイクロアレイ実験において異なる調節を受けた遺伝子の検出に用いられます。これによりHECToRスーパーコンピューティング設備や他の並列マシン上の利用が可能になります。
Rパッケージ (R Development Core Team, 2010)は、統計解析と可視化のオープンソフトウェアパッケージです。その実装の多くはR言語で記述されています。
The Simple Parallel R Interface：SPRINT (Hill et al., 2008; Petrou et al., 2010)は、エジンバラ大学のDivision of Pathway MedicineおよびEPCCによるプロジェクトとして、Rを用いたマイクロアレイ解析に対する並列ワークフローを提供することを目的としています。その狙いは、できるだけ既存のコードに近い機能を提供することで、エンドユーザである生物統計学研究者が自身のワークフローを大きく変えることなく、HPC環境を利用可能にすることです。
rank product法
これは、マイクロアレイ反復実験において異なる調節を受けた遺伝子を検出する方法です。例えば（処置済みと非処置といった）2つの異なる実験条件下の遺伝子発現レベルを見る場合です。初期のやり方では、サンプルAとBの発現レベル比である倍率変化に基づく手法（fold-change criterion, DeRisi et al., 1997）が用いられました。しかしながらこの方法は有意水準の計算が出来ず、fold-changeのカットオフをどうとるべきかも明らかではありません。rank product法はfold-change criterionを用いて構築されていますが、それを反復実験へ適用します。つまり一つのサンプルが各実験条件（クラス）に対応するのではなく、各クラスに複数のサンプルを持つことになります。
rank product法は2つの異なる実験条件（これをクラスAとクラスBとします）の比較実験に応用することが出来ます。各遺伝子について、クラスBに対するクラスAの全てのペアごとの比較により遺伝子のfold-change値をランキングし、全てのサンプルに渡ってこれらランク値の積を計算します。次のステップはこのrank productの帰無分布の計算です。これは遺伝子間とサンプル間どちらにも差が生じない場合に期待される分布です。実験で観測される各遺伝子のrank productを帰無分布と比較して、有意水準の正確な計測とカットオフ値の推定を行います(Breitling et al., 2004)。
·アルゴリズムの詳細
Ng個の遺伝子とNs個のサンプルを考えます。全サンプルは、Na個のクラスAとNb個のクラスBから構成されるとします（Ns= Na+Nb）。ここで、クラスAとクラスBの全ての対比較のfold-changeの（総数はNa×Nb）行列を構築します。
ここで行列要素fi,jは、遺伝子iの対比較jでのfold-changeです。ただしこのステップはシングルチャネルの場合に必要な物です。二色の場合は（発現レベルでなく）fold-changeを直接計測します。この場合は2種のクラスでなく一つのクラスで表現されますが、同じやり方で解析が可能です。入力データを直接fold-change行列として扱うかどうかのみが異なります。
次に各サンプルのfold changeを、最大から最小へ（アップレギュレーションの場合）あるいは最小から最大へ（ダウンレギュレーションの場合）へランキングします。最後に、適切に正規化された全てのサンプルに渡る特定の遺伝子のランクの積を計算すれば、その遺伝子のrank productが得られます。
実験的に観測されたrank productの有意水準を得るには（つまり強くアップレギュレーションされた遺伝子を決定する）、この実験的に観測れた値と帰無仮説下のrank productの推定分布を比較する必要があります。帰無仮説は差異的な遺伝子発現がないこと、つまり個々の遺伝子の発現レベルは互いに独立で同じ分布から生じ、全サンプルは独立である（つまりサンプルAの遺伝子発現レベルはサンプルBでの発現レベルと相関しない）事を意味します。残念ながら、帰無分布の解析的表現は困難ですが、ここではブートストラップ法で数値的に構築します。各サンプルの遺伝子発現ベクターを独立に交換したランダム実験を行い、このランダムデータから全遺伝子のrank productを計算します。これを何度も繰り返して、帰無仮説に対するrank productの分布を構築します。
シリアル実行でrank product計算を行うのは2つの理由で低速です。最初の理由は解析対象データが極めて大きいことです。実行時間を低減するためには、一つのデータセットに対するrank product解析計算を高速化する必要があります。2番目の理由として、多くのブートストラップ標本が要求されるため、この場合に実行時間を削減するにはブートストラップ解析計算を高速化しなくてはなりません。こうした観点から異なる並列手法に自然と至ります。
·並列化
rank productの高速化には、明らかな2種の並列化手法があります。最初のものは、一つのデータセットに対するrank product解析へデータ並列を実装するものです。第2の手法は、rank productはシリアル実行し、ブートストラップ標本毎に並列実行する（タスク並列手法）ものです。
·データ並列
マイクロアレイ実験では、通常は発現観測数は1万〜数千万ですが、サンプル数は1から1000程度です。
サンプル上の並列を行って入力データ行列を分散させると、膨大な通信が発生します。
遺伝子上の並列化の場合は、データの移動は少ないですが、fold-change行列を入力データから生成しなくてはならないため新たなボトルネックとなります。遺伝子上の並列化を行うと、fold changeはローカルに計算されますが、fold changeのランキング処理に並列ソートが必要です。並列ソートには様々な手法があります（初期の概説はAkl (1985)を参照してください）。しかしながらクイックソートのような再帰的手法は、適切なピボットの選択には通信が欠かせないため負荷バランスが困難です。
サンプルベースの並列化よりも遺伝子上の並列化実装の方が複雑性が大きくなります。
·タスク並列
単純な並列化手法は、ブートストラップ標本生成とそのrank productの帰無分布生成を並列に行うやり方です。その後、観測したrank productの統計的性質を測るためにこれら部分的な分布を全分布へ結合します。
タスク並列にしたうえでデータ並列を行う理由の一つに、データが一つのプロセスのメモリーに収まらないことが有ります。一見、ここでのアルゴリズムではfold-change行列とブートストラップの分布によりメモリーのほとんどが占められているように見えます。しかしながら、fold-changeを再計算しながらそのランキングを計算するよりも、一つのfold-changeのランキングだけを行って、各遺伝子についてのrank productの逐次積を計算していくことも可能です。この場合は入力データに対する追加のメモリーはNgオーダーのテンポラリー領域のみです。
ブートストラップ標本について興味がある統計量は、ブートストラップ帰無分布を遺伝子ｉの観測rank productまで積分したものです。計算の最後で全分布を集めて積分をするのではなく、ブートストラップ標本毎に遺伝子を入れ替えて計算していけば可能です。こうしてメモリーもNgオーダーへ削減できます。
大きなマイクロアレイデータでは、100サンプルとサンプル当り100万の発現データを持つ場合があります。そのような場合にはデータ並列手法の必要性は認められないと判断しました。その代わりにここでは並列ブートストラップ手法を実装しました。
性能
ここではHECToRスーパーコンピューティングサービスのXT4構成でベンチマークを行いました。前回のプロジェクト(Mitchell, 2011)では、XT4とXE6構成に対してSPRINTの通信パターンに性能差が示されていました。使用したテストデータは典型的なマイクロアレイデータで、2つのクラスが含まれる62サンプル（クラスAが35、クラスBが27）の23292遺伝子の発現レベル観測値で構成されます。2つのクラスのrank product計算ではそのランキングに945の対比較が必要です。1クラスのrank productの振舞いを調べる際には、クラスAサンプルのみを用いました。アルゴリズムの複雑性はランキング対象サンプル数に比例するため、1クラスrank product計算は2クラスに比べて約27倍高速であると予想されます。
·推定性能
ブートストラップの並列処理は、スレーブプロセスへ入力データをブロードキャストするところから開始され、rank productの部分分布の計算を通信を用いず行い、最後に部分的な結果をマスタープロセスへ収集する際に通信を行います。よって、スケール性は良好であることが推測されますが、プロセス毎のブートストラップ標本数が小さい場合はスケール性能も落ちると推測されます。
この並列実装コードの経過時間には2つのシリアル計算が含まれます。このシリアル計算は、アップ及びダウンレギュレーションに対するrank product実験値の計算です。この実装では、アップおよびダウンレギュレーション遺伝子に対するブートストラップ標本分布計算を並列に実装します。このため理論並列効率はその分の事前のペナルティを持ちます。
·実際の性能
最初に1クラスを用いた場合のコード性能を検証します。この場合のサンプル数はクラスAの数で、今は35です。1024および16384サンプルをブートストラップで解析しました。シングルコアにおいて、1024サンプルは1100秒、16384サンプルは17300秒掛かりました。既存のRankProdコードは1024サンプルで2100秒掛かっていました。結果は、1024サンプルでは512プロセッサまで、16384サンプルでは4096プロセッサまで、理論性能にスケールしました。
また2色データセット（35および27）では、対比較数が増加するため大幅に時間が掛かります。この場合はrank product解析において945個の対比較が必要です。その実行時間の増加は27倍と推定できます。これも同様に1024および16384サンプルをブートストラップで解析しました。シングルコアにおいては、1024サンプルは26000秒、16384サンプルは434000秒掛かりました。既存のRankProdコードは1024サンプルで42300秒掛かっていました。結果は、1024サンプルでは512プロセッサまで、16384サンプルでは8192プロセッサまで、理論性能にスケールしました。
ここで、プロセスを1〜4へ増やす際に性能の劣化が存在ます。この原因はシングルノード上のメモリー競合によるものです。
バイオインフォマティクス研究者がブートストラップ標本数を選ぶ理由の一つが、その計算時間です(Dunbar, 2011)。しかしながら高い精度でp値を計算するには、より大きなブートストラップ標本数が必要です。実験P値の計算においてN個のブートストラップ標本を用いた場合は、その誤差は典型的にはNの平方根の逆数に比例します。よって例えばp値を3桁まで正確に得たい場合は、百万個のブートストラップ標本比較が必要です。ただしこれは百万個のブートストラップ標本数を意味しません。各ブートストラップ標本は実験値と比較すべきNg個のrank productを生成するため、もし1万の遺伝子を持つ場合は百万の比較を得るのに100個のブートストラップ標本で足ります。
·巨大データセットを用いた場合の性能
将来のマイクロアレイ実験データは、エクソンのマーカーや一塩基多型(SNPs)といった各サンプルにおいて105の発現レベル数以上に達します。もし発現レベル数を倍にすると、実行時間は2(log2Ng)/(logNg)倍に増加します。これはアルゴリズムの複雑性が遺伝子数に対してO(Nglog Ng)であるためです。
ベンチマーク計測結果から、我々のコードは遺伝子数に対してアルゴリズムの複雑性に正確にスケールします。このことから、このコードは並列プロセス数の増加に対して実行時間を妥当なレベルへ抑えることが出来ます。
謝辞
このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | S. G. Akl. Parallel Sorting Algorithms. Academic Press, Orlando, Florida, 1985. |
| [2] | R. Breitling, P. Armengaud, A. Amtmann, and P. Herzyk. Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments. FEBS Letters, 573:83.92, 2004. doi: 10.1016/j.febslet.2004.07.055. |
| [3] | J. L. DeRisi, V. R. Iyer, and P. O. Brown. Exploring the metabolic and genetic control of gene expression on a genomic scale. Science, 278(5338): 680.686, 1997. doi: 10.1126/science.278.5338.680. |
| [4] | D. Dunbar. Personal communication, 2011. |
| [5] | D. R. Helman, J. Jájá, and D. A. Bader. A new deterministic parallel sorting algorithm with an experimental evaluation. Journal of Experimental Algorithsm, 3, 1998. |
| [6] | J. Hill, M. Hambley, T. Forster, M. Mewissen, T. Sloan, F. Scharinger, A. Trew, and P. Ghazal. SPRINT: A new parallel framework for R. BMC Bioinformatics, 9(1):558, 2008. |
| [7] | F. Hong, R. Breitling, C. W. McEntee, B. S. Wittner, J. L. Nemhauser, and J. Chory. Rankprod: a bioconductor package for detecting differentially expressed genes in meta-analysis. Bioinformatics, 22(22):2825.2827, 2006. doi: 10.1093/bioinformatics/btl476. |
| [8] | L. Mitchell. A parallel random forest implementation for R. Technical report, EPCC, 2011. |
| [9] | S. Petrou, T. M. Sloan, M. Mewissen, T. Forster, M. Piotrowski, and B. Dobrzelecki. Optimization of a parallel permutation testing function for the SPRINT R package. In Proceedings of the 19th ACM International Symposium on High Performance Distributed Computing, HPDC ’10, pages 516.521, 2010. doi: 10.1145/1851476.1851551. |
| [10] |
R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2010. URL http://www.R-project.org. ISBN 3-900051-07-0. |
| [11] | H. Shi and J. Schaeffer. Parallel sorting by regular sampling. Journal of Parallel and Distributed Computing, 14:361.372, 1992. |