1.Introduction

Microarray data enables the high-throughput survey of mRNA expression profiles in genomic level. At the same time, it offers a challenging statistical problem due to the large number of genes surveyed with small sample of sizes. Identification of differentially expressed (DE) genes between two or more user defined groups is an important task to reduce the dimensionality of microarray data. There are several classical and Bayesian or empirical Bayes approaches for identification of DE genes. However, given the complexity of the microarray data, there are no models those can explain the data fully. It is generally difficult to scrutinize the irregular patterns of expression or contaminated genes that are not expected by the models gene by gene.

A statistical framework to detect irregular patterns of expression or contaminated genes and diagnose the model may reduce this difficulty. Inference about deferential expression is a typical objective in analysis of gene expression data. Bayesian approaches have become increasingly popular for detection of DE genes from microarray data. However, most of these approaches are very sensitive to outlier and produces misleading results. Therefore, as an extension of empirical Bayes (EB) procedures, I developed β-EB approaches assuming (i) constant gene-specific variance and (ii) variable gene-specific variance for the identification of DE genes. Also an attempt is made to extend the β-EB LNN approaches for paired gene expression data analysis. The proposed β-EB approaches are unique parametric approaches because, not only it is robust against outliers, but it also detects contaminating genes and statistically diagnoses gene expression profiles.

To robustify classical EB-approach, I maximize β-likelihood function using EM like algorithm, where the β-likelihood function is induced from the β-divergence. The proposed robust β-EB approaches introduce weight function, which I call β-weight function. The weight of a transcript t is described as a power function of its likelihood, fβ(yt|θ). Genes with low likelihoods have unexpected expression patterns and low weights. By assigning low weights to outliers, the inference becomes robust. The value of β, which controls the balance between the robustness and efficiency, is selected by maximizing the predictive β0-likelihood by cross-validation. The distribution of the weights is used to scrutinize the irregular patterns of expression and diagnose the model statistically.

2.Outline of the thesis

In chapter 1, I introduced microarray gene expression data, problem of the study and objective of the study in details.

Chapter 2 discusses the microarray technology for generating the gene expression data. This chapter also discuss several classical and Bayesian methods for identification of DE genes, where EBarrays which I call EB in this thesis is one of the most popular approaches for identification of DE genes. However, most of the existing algorithms are not robust against outliers and none of them can detect contaminated genes from high-dimensional gene expression datasets. To overcome this problems, in this thesis, I proposed some modification of classical EB approaches by maximizing β-likelihood function using EM like algorithm for robust statistical inference those are discussed in the subsequent chapters of this thesis.

Chapter 3 introduced β-EB approach for robust identification of DE genes assuming constant gene-specific variance. Numerical simulation showed that the contaminated genes are detected by inspecting the values of the β-weights. In the absence of contaminated genes, β-EB and the other existing procedures had similar performance. When the data includes contaminated genes, β-EB and not the others were robust.

In chapter 4, I developed β-EB approach for robust identification of DE genes assuming variable gene-specific variance. To see the performance of the estimation and model diagnosis, I conducted two types of simulations. The first simulation compares the performance for different sizes of data and different levels of outliers. The second simulation generates gamma distributed expression profiles, whereas estimation procedures assume log-normal distributions. When the Gamma distribution has the shape parameter < 1, it has a large mass near the value of zero, and will not be approximated by the log-normal distributions. Again, it was shown that the β-EB approach, and not the others, is robust against outliers. By comparing the β-weight distribution with the predicted distribution, it was possible to detect the genes with expression profiles that contradict log-normal distribution.

In chapter 5, I analyzed three sets of real gene expression data (head and neck cancer, lung cancer and Arabidopsis thaliana) using both classical and proposed EB-LNNMV approach. In the analysis of head and neck cancer data, the β-EB approach detected six contaminating genes (LRP8, S100A8, S100A9, TRIM29, CSTA, ACP5) as outliers with the posterior probability of DE > 0.95; the posterior probability for these genes by the classical EB-LNN approach was < 0.5. Inspection of the expression profiles of these genes confirmed the presence of outliers. The classical EB approach had the low power due to over estimation of variances within the groups. In the lung cancer data, the β-weight distribution deviated largely from the predicted distribution, and implied the sign of model misspecification. The analysis of scatter plot showed that this is due to the genes with little expression and the genes with large within variance. By excluding the genes with extremely low expression levels, the β-weight distribution became consistent with the model accompanied by a few outliers. When applied to the eQTL analysis of Arabidopsis thaliana, the β-EB approach gave on average larger numbers of regulated genes compare with classical EB approach. Furthermore, the proposed β-EB approach identified some potential master regulators that were missed by the EB approach. They include markers on a telomeric region of chromosome 4. This region includes three transcription factors one of which is CYC1 (cyclin-dependent protein kinase regulator).

Chapter 6 presents the modification of β-EB approaches for identification of DE genes in the case of paired genes expression data. Simulation results show that the performance of the proposed method is good for identification DE genes from paired observations.

Chapter 7 presents the overall conclusion.

3. Conclusion

In my thesis, I have discussed the robustification of EB approach by β-divergence assuming both constant and variable gene-specific variance. The proposed method reduces to the standard EB approach for β→ 0. The performance of the proposed method in a comparison of the classical EB approach, t-test for identification of DE genes investigated using AUC and pAUC in the simulation study. From the simulation results, I observe that the proposed method significantly improves the performance in a comparison of the others in presence of outliers; otherwise, it keeps almost equal performance.

Simulation and real gene expression data analysis results show that the performance of the proposed method much better than the other existing methods in presence of irregular gene expression patterns. Otherwise, it shows almost equal performance.

I extended our proposed β-EB approach for identification of DE genes from the correlated expressions between two user-defined groups. I investigate the performance of this approach using simulated data only. Therefore, I would like to apply this update version to real paired gene expression dataset soon to investigate the performance from the robustness point of view.

1.問題の所在

マイクロアレイは数千ないし数万の遺伝子の発現プロファイルを鳥瞰することを可能とする。もっとも基本的な解析は、条件により発現量が異なる遺伝子を検出することである。調査する遺伝子数が標本サイズを格段に上回るため、発現量が遺伝子間で異なる様子を確率分布で記述し、階層ベイズモデルで表現することにより、検出力を高める方法も開発されてきている。遺伝子間の関連を表現するネットワークを推定する方法も開発されている。しかしながら、これらハイスループットな統計的推定・検定手法は、例外なく遺伝子発現のパターンに対して何らかの仮定をしている。その仮定が成立するときは有効な方法となる。他方、何らかの理由でデータが想定外の発現パターンを含むときは、検出力が低下し、また偽陽性を多く拾う危険性がある。ところが、数多くの遺伝子を分析の対象としているため、個々の遺伝子を詳細に調べるのはほぼ不可能である。データには分析の前提条件を満たさない想定外の発現プロファイルを持つ遺伝子が含まれているのか、含まれているとするとそれはどのくらいの数か、数遺伝子なのか数十遺伝子なのか、あるいは数百、数千のオーダーなのか、これまで調べる術がなかった。

2.β尤度法とβ荷重

本論文は、グループ間で発現に違いがある遺伝子を検出する方法として、近年頑健な推定法として提案されたβ尤度法を階層ベイズモデルに適用する。遺伝子のβ荷重を計算し、その分布をモデルから予測される分布と対比する方法を提案し、シミュレーションと実データの解析によりこの方法の有効性を証明する。最尤法は、データの背後にある確率構造の統計的モデリングを行い、その前提の下でデータが生成される確率(尤度)を最大にするよう、パラメータを推定する。すなわちデータに最もなじむ確率構造を推定する。本論文で提案するβ尤度法は、回帰分析において対数変換を一般化させたBox-Cox変換に着想を得て、対数尤度を一般化させたものである。微分をとったスコア関数を見ることにより、各遺伝子が異なる重みを持たせた荷重尤度法と同等であることがわかる。遺伝子tの重み(β荷重とよぶ)は尤度ftをべき乗したfβtとなる。従って、モデルの想定外の発現プロファイルを持つ遺伝子は尤度が小さく、小さな重みを持つため、解析結果にあまり影響を与えない。すなわち、β尤度法は、異常値に影響されない頑健な推定法であることが期待される。交差検証法により頑健性と検出力をバランスさせるβの値を決める。さらにβ荷重の遺伝子間の分布をモデルが妥当する場合に期待される分布と対比させることにより、異常値を統計的に検出し、モデルの妥当性を診断することが可能となる。

3.シミュレーションによる有効性の検討

標本サイズが60の中規模のデータと標本サイズが20の小規模なデータについて、2通りのシミュレーションを行い、提案手法の有効性を検証した。第一のシミュレーションでは異常値による頑健性を提案手法と既存の手法を比較した。偽陽性のコストを勘案した真陽性の検出力を表すAUCおよびpAUC、グループ間で発現に差のある遺伝子の割合(ここでは混合率という)の推定精度を調査したところ、モデルの前提が成立し、異常値を含む遺伝子がない場合には、どの手法も良好な成績を示した。しかし、異常値を含む遺伝子が存在する場合には、既存のどの手法もAUCおよびpAUCが減少し、混合率も大幅に過大推定された。これに対しβ尤度法は、中規模のデータ、小規模のデータいずれにおいても、異常値による性能の低下は見られず、混合率も偏りなく推定することが示された。第二のシミュレーションでは、β荷重の分布によりモデルの妥当性を診断することの可能性を調査した。ガンマ分布に従う発現プロファイルデータに対して、対数正規分布を当てはめる場合について実験を行った。形状パラメータが小さい遺伝子では微小な発現が優先するため対数正規分布で近似することができず、β荷重が極端に小さくなることが示された。

4.実データの解析

公開されている3つのデータについて分析を行った。第一のデータは頭頸部癌の患者22人の癌組織と正常組織を対比したデータで、12625の遺伝子中2.2%にあたる261の遺伝子がp<10-5で有意に異常な発現プロファイルを持つ遺伝子として検出された。従来法では95%の確率で差なしと判定するにもかかわらず、β尤度法は逆に95%の確率で差ありと判定している遺伝子が6つあった。これらはLRP8、S100A8、S100A9、TRIM29、CSTA、ACP5で、癌との関連が報告されている。いずれも有意にβ荷重が小さく、発現プロファイルをプロットしたところ、異常値が検出された。第二のデータは肺癌患者からの54675のRNA転写産物で、40人が腺癌、18人が扁平上皮癌である。β荷重の分布は予測分布に比し上側と下側に裾が重く、モデルの妥当性が疑われた。そこで平均発現量と分散を調査したところ、有意にβ荷重が小さい遺伝子は分散が大きく、逆に有意にβ荷重が大きい遺伝子は発現に変異がないことが示唆された。後者の遺伝子をはずして再解析したところ、概ねデータがモデルと異常値で説明されることがわかった。第三のデータはシロイヌナズナBay0×Sha組換え近交系211株からとられた22810プローブの発現および578のSFPマーカー遺伝子型の情報である。ここにおいてもβ荷重の分布は予測分布と大きくかい離し、異常値を伴う発現プロファイルを数多く検出した。従来法によるeQTL解析は、発現を制御される遺伝子の数を大幅に過小推定する可能性があることが示された。

システム生物学を支える技術の革新に伴い、質的にも量的にもデータ自身の持つ情報量が飛躍的に膨らんできている。トランスクリプトーム、メタボローム、プロテオーム、およびそれらの統合するデータ解析の手法も急速に進歩している。一方で、データの精査がかつてなく困難な問題として突きつけられている。本論文が提案するβ尤度法とそこから派生するβ荷重の分布によるアプローチは、膨大で複雑なデータを二段階接近法により精査する枠組みを提供しており、今後有効性が立証されることが期待され、学問的にも応用的にも貢献するところが大きい。よって審査委員一同は本論文が博士(農学)の学位を受けるに十分な価値があると認めた。