マイクロアレイ解析によって網羅的な遺伝子発現データを得ることができる。また https://www.omicsdi.org/ Omics DI などのデータベースで様々なデータを得ることができる。GSEA ( Gene Set Enrichment Analysis) は、それらのデータから生物学的に有用な知見を得るために有効に使える方法としてよく知られている。この方法は、遺伝子一つ一つが持つ値よりも、一群の遺伝子（共通した性質を持つ遺伝子グループ）が形成する分布の違いに注目するという点で、マイクロアレイの利点を生かした情報が得られる。
GSEA に関しては、既にすばらしい解説が公開されている。
門田幸二 博士による解説 http://www.iu.a.u-tokyo.ac.jp/~kadota/
GSEA に近い方法には様々なバリエーションがある。新しい方法もさらに開発されている。 http://www.riken.jp/pr/press/2016/20160510_1/ 「遺伝子発現から転写因子を予測」では、理化学研究所の川上英良特別研究員らのグループによって、従来のGene Set Enrichment解析に確率的な制御関係を考慮したwPGSA法（weighted Parametric Gene Set Analysis） という方法が開発されたことが紹介されている。 http://wpgsa.org/
Gene-set distance analysis (GSDA): a powerful tool for gene-set association analysis. Cao X, Pounds S. BMC Bioinformatics. 2021 Apr 21;22(1):207. doi: 10.1186/s12859-021-04110-x. PMID: 33882829 という論文も発表された。
一般に GSEA 法と普通に呼ばれるのは、Broadinstitute による
http://www.broadinstitute.org/gsea/index.jsp Gene Set Enrichment Analysis のページ
このページで扱われている方法が世界的な標準になっているらしい。Broadinstitute は世界的に名高い権威ある超一流大学であるハーバード大学とマサチューセッツ工科大学の共同研究施設だそうである。だから普通の研究者は Broadinstitute のページを見るべきで、以下に書いてあることを見る必要は全くないし参考にもならない。
一応自分でやっていることも記録しておく。GSEA に近い方法には様々なバリエーションがある。私は、門田博士が公開されていたかなり以前の資料を参考にして、
PAGE: parametric analysis of gene set enrichment. Kim SY, Volsky DJ. BMC Bioinformatics. 2005 Jun 8;6:144. PMID: 15941488
という論文に基づいて分析している。この論文の計算法はとてもシンプルで簡単に計算できる。それで私にとっては特に問題なく分析できているように思える（もちろん今ではもっと高度な改良された方法が多数発表されている）。 この論文の被引用数はとても大きいので、手法に問題がなく有用であることは一応コンセンサスとなっているのだろう（間違いがあり役に立たないなら誰も引用しない）。 メタボロームデータにも同様に適用することができる (Metabolite set enrichment analysis)。株) ヒューマン・メタボローム・テクノロジーズ によって資料、論文が公開されている。 http://humanmetabolome.com/rd/proceedings
Statistical hypothesis testing of factor loading in principal component analysis and its application to metabolite set enrichment analysis. Yamamoto H, Fujimori T, Sato H, Ishikawa G, Kami K, Ohashi Y. BMC Bioinformatics. 2014 Feb 21;15:51. doi: 10.1186/1471-2105-15-51. PMID: 24555693
Comparative metabolomics charts the impact of genotype-dependent methionine accumulation in Arabidopsis thaliana. Kusano M, Fukushima A, Redestig H, Kobayashi M, Otsuki H, Onouchi H, Naito S, Hirai MY, Saito K. Amino Acids. 2010 Oct;39(4):1013-21. doi: 10.1007/s00726-010-0562-y. Epub 2010 Mar 31. PMID: 20354740 という論文でも触れられている。
この PAGE 法は、Parametric analysis と呼ばれるタイプに属する（名前がそうなっている）。これは、Zscore が標準正規分布するということがこの方法の基盤、根拠になっていることによるらしい。そのため素人にもわかりやすい。 Broadinstitute の GSEA は、そうではないらしい。http://arxiv.org/ftp/arxiv/papers/1110/1110.4128.pdf という文書の最初の方に「GSEA’s nonparametric nature 〜」と書いてある。Broadinstitute の方が高度で複雑で専門的で、得られる結果の意味も PAGE と異なる。
Kim ら (2005) の論文では、predefined gene sets に対する Zscore を計算する。この場合の Zscore は入試などでよく出てくる偏差値に相当する、偏りを表す値である。試験の場合になぞらえると、『試験を受けた多数の学生の成績から「すべての学生の成績」と「ある大学を志望する学生のグループの成績」を取り出す。両者の分布を比較して「ある大学を志望する学生のグループ」では成績がどれくらい偏っているかを数値化する』というようになる。遺伝子の場合は「すべての遺伝子の値」と、「恣意的に決めた Gene set に含まれる遺伝子の値」を取り出す。
（Zscore と名前がついていても、その定義、計算法はいつも同じとは限らない。「一つの分布を平均０，標準偏差１に変換（Z 変換）して、その分布に含まれるあるサンプルの値が平均からどれくらい離れているかを表す値をそこから得る http://www.env.nagasaki-u.ac.jp/sekkei/faculty/nakayosi/excel/yougo1.html の 標準化」場合や、PAGE 法のように「二つの分布（母集団の分布、そこから何らかの基準に基づいて取り出された、標本となる一群のデータの分布）を比べてどれくらい偏っているかを表す値を得る http://www.env.nagasaki-u.ac.jp/sekkei/faculty/nakayosi/excel/yougo1.html の Z検定その１」という場合などがあるらしい。Z検定その２と言うのも解説されているので、すくなくても３種類の Zscore があることになる。専門家は「どれも本質的には同じものだ」と言うのだろうが、紛らわしくて困る。しかしそういうことになっているらしいのでしかたがない。）
Gene set の方は遺伝子の機能、発現様式などを元に研究者が arbitrary 恣意的に決める。そこに、その研究者が持っている生物学的な知識、実験結果からの情報を生かす・取り込むことができる。 既に公開されているアレイデータや発現データベースから得られる知識を用いることができる。「既に蓄積されている実験データ、知識を取り込むことができるように理論や解析法ができている」ことは大切なことらしい。（とはいっても、自分にとって都合がよい結果になるように Gene set をつくるようではいけない。もっともらしい基準は必要である）。生物学に詳しくない人は GO などの出来合いのリストを使えばよい。しかしそういうリストを見ていると「何でこの遺伝子がこのリストに入っているのか？ どうしてこの遺伝子が入っていないのか？」と思うことがかなり多い。生物学者から見れば Gene set をいかに作るかが一番重要なポイントで、計算法は簡単ですむならそれが一番よい。
すでに Gene set をデータベース化することが行われている。２０１２年の分子生物学会の講演で荒木博士が発表されている。 http://genesetdb.auckland.ac.nz/haeremai.html 動物のデータしか入っていないらしい。
3W2III-5 Gene Set Analysis のための 統合データベース, GeneSetDB, の構築 荒木 啓充（The University of Auckland）
http://www.broadinstitute.org/gsea/index.jsp から、Molecular Signatures Database v4.0 に進むと、gene set のデータベースを見ることができる。Gene set を作成する基準は gene ontology に限るわけではなく、様々な見地から作成された gene set が define されているらしい。
c1 positional gene sets for each human chromosome and cytogenetic band. 染色体上の遺伝子の位置が近い遺伝子でセットを作る。
c2 curated gene sets from online pathway databases, publications in PubMed, and knowledge of domain experts. 代謝やシグナル伝達の経路、ネットワークに基づいてセットを作る。
c3 motif gene sets based on conserved cis-regulatory motifs from a comparative analysis of the human, mouse, rat, and dog genomes. プロモーターにあるシスエレメントモチーフを共通して保持する遺伝子でセットを作る。
c4 computational gene sets defined by mining large collections of cancer-oriented microarray data. マイクロアレイデータからセットを作る。
c5 GO gene sets consist of genes annotated by the same GO terms. GO terms に基づいてセットを作る。
c6 oncogenic signatures defined directly from microarray gene expression data from cancer gene perturbations. がんと関係しそうな発現をしている遺伝子でセットを作る。
c7 immunologic signatures defined directly from microarray gene expression data from immunologic studies. 免疫と関係しそうな発現をしている遺伝子でセットを作る。
やり方： エクセルと R 言語を用意する。
まず、分析したいマイクロアレイデータを用意する。 https://www.omicsdi.org/ Omics DI は検索しやすい。ここでサーチして GEO database に進むとやりやすい。
処理区、対照区がセットになっている必要がある。cel ファイルを元に R言語の RMA 法で処理したデータにするとよい結果が得られるような気がする。RMA 法で処理したデータは、発現量が対数で表現されていて 3〜12 くらいの値になっている。
処理区の値を対照区の値で割り算し、２を底とした対数にしてそれぞれの遺伝子の fold change を出す。RMA 法の結果なら最初から対数なので引き算すればよい。
データ全体の fold change の平均と標準偏差を計算しヒストグラムを書く。Kim らの論文では Fig. 1 にヒストグラムがある。 普通のデータなら、０を中心として左右対称に近い分布を示す。０で鋭くとがったピークになっているように見える。データ全体での fold change の平均値は０に近い値になっているはずである。そうなっていなかったら何かおかしいので分析できない（この段階で偏りがあってはいけない）。これが母分布 (parent distribution) になる。Fold change を指標とする際は MA plot というプロットをして、fold change の値に発現量による偏りが少ないことを確認することが必要であると、門田先生が書かれている。これはとても重要なことらしい。RMA 法で処理した結果では偏りが比較的小さくなりやすいらしい（発現が最小のあたり、最大のあたりで fold change が０に偏りやすいようだが、これはそれほど問題にならない。最小値のあたりはデータとして信頼できないので０に偏っていてもよい。最大のあたりで０に偏りやすいのはあまりよくないが、植物だと発現量が高いものは光合成関係というように特定の機能に関連したものが多い。そういう特定の機能と関係ないことを探っているのなら影響は少ない）。偏りを指標として分析するので、元のデータ自体が偏っているとよくない。
ここで出てくる分布は、指数分布を０で折り返したような形に見える。こういう分布をラプラス分布と呼ぶそうである。 http://ja.wikipedia.org/wiki/%E3%83%A9%E3%83%97%E3%83%A9%E3%82%B9%E5%88%86%E5%B8%83 何か意味があるのかもしれない。しかしデータの前処理の方法（RMA 法) の影響が大きいらしい。 データを RMA 法で処理する際には、複数の .cel ファイル（例えば処理区から得たデータファイル３つ、対照区から得たデータファイル３つ）を一つのフォルダに入れる。そこに対して RMA 関数を適用すると複数あるデータファイルの分布の形ができるだけ一致するように補正される。そのため、処理区／対照区の値は 1 になっている場合が多くなる。 画像処理は各ピクセルの値が二次元に並んだデータを分析することで行われる。マイクロアレイデータは一次元だが、データ分析という観点からは画像処理と似ている。画像処理について書かれている文章を眺めていたら、「レナ画像の Xi - Xj のヒストグラム」が掲載されていた。そのヒストグラムは、０の部分でかなり尖っていた。画像に、正規分布するノイズが重なった場合に原画像を推定する方法が解説されていた。理屈的には、マイクロアレイデータでも似たことができる。RMA 法はそれに相当するのだろうが、まだよく理解していない。 画像データとマイクロアレイデータで同じ分布が出現するのなら、普遍性が高いということで、いろいろなことを考える基盤にするのによいかもしれない。 指数分布のタイプの式が出てくる分布に、カノニカル（ボルツマン）分布と言うものがある。 気体のなかの一個の分子がエネルギー E の状態にある確率 P(E) を表すカノニカル（ボルツマン）分布に関する、わかりやすい説明 http://www.h5.dion.ne.jp/~antibody/boltzmann.htm （抗体科学研究所 http://www.h5.dion.ne.jp/~antibody/index.htm 内のページ） 容器内の気体分子に重力が働いていると考えると、そうでないときの一様な分布から、指数分布に変化する。 遺伝子になぞらえると「処理区／対照区の値が１からずれると、ずれればずれるほど１へ戻そうとする力が強くかかる」状態では指数分布が出てくることになる。ここで出てくる分布の幅がとても広かったりしたらそれは「温度が高い」というように見立てることもできる。それはそれで何か情報を含んだ値として使えるのかもしれない。しかし元になっている RMA 関数の中身をよく理解していないので説得力はない。 「１からずれればずれるほど、１に戻そうとする力が強くかかる」というのが本当なら、バネの伸びによって生じる力と似ている。バネの先端におもりをつけた物を考える。その運動を考えるときは、おもりが運動する加速度におもりの質量をかけ算した力（加速度による力）、(速度に比例する) 摩擦力、バネの伸びに比例する力を考える。 「微分方程式で数学モデルを作ろう」という本に丁寧に説明してある。遺伝子も、それぞれの遺伝子にバネが付いていて「１からずれればずれるほど、１に戻そうとする力が強くかかる」ようになっていると見立てることもできる。 「すべてのものはバネである」ということは普遍的に成り立つそうである（一定の条件が満たされていれば）。 http://www008.upp.so-net.ne.jp/takemoto/index.htm 竹本先生のページで、「普遍的に成り立つフックの法則」という記事で紹介されている。 遺伝子の状態もその一つに入るように考えてよいのなら、いろいろなことを考える切り口になるかもしれない。問題点として、「本当にラプラス分布なのか」「前処理に使っている RMA法の中身を理解していないじゃないか」などがある。
ここで、母分布（データ全体）からランダムに幾つかの遺伝子を抽出し、それらに対応する値を取り出すとどうなるかを考えてみる（母集団から標本を抽出）。
Kim らの Fig. 1C では、10 個の遺伝子をランダムに抽出してそれらの fold change を平均した値（標本平均）を計算することを繰り返し（何回？）、出てきた値の分布を示している（ように思える）。
その結果は、母分布の平均（この場合０）を中心として狭い範囲に収まった正規分布になっている。正規分布は にいろいろくっついた様な形をしている。これの傾きを求めるために微分する。これは合成関数の微分と呼ばれるパターンになっている。
=
=
そのため x = 0 では傾きが０になりなだらかになる。
遺伝子数は 10〜50 で同様な結果になったと書かれている。
次に、研究者が何らかの基準によって用意した遺伝子リストを用い、同様に対応する値を取り出してみる。 遺伝子リストに含まれる遺伝子の数は、Kim らの結果によると最低で 10 個は必要である。私は 30〜50 個くらいにしている。スコアを計算するだけでなくヒストグラムを書くことで偏りを確認しないといけない。
遺伝子が20000 種類あるとする。20000 種類の遺伝子から30 種類を取り出す組み合わせは、「順番のない n個から、r個を取り出す（これも順番は考えない）方法」の公式を当てはめて、20000 C 30 = 20000! / ((20000-30)! * 30!) 通り存在する。 Excel の関数なら =COMBIN(20000,30)、R の関数なら choose(20000,30)で、答えは 3.96E+96 という大きな数になる。
GO のリストなら「XX機能、経路に関連する遺伝子のリスト」のようになる。 マイクロアレイのデータを元にするなら、「＠＠処理によって強く誘導される遺伝子のリスト」などが考えられる。それらは共通した性質を持つ遺伝子の一群と考えられる。
遺伝子というものはたくさん種類があり、それぞれ独自の性質をもつ。それをきちんと考えようとすると大変である。 「共通した性質を持つ遺伝子」というグループ、リストをつくれば、そのグループ内ではそれぞれの遺伝子を区別する必要が小さくなる。 「粗い情報」にすることができる。 それによって、気体を構成する分子の運動に関する統計的な考え方（一つ一つの分子を区別することはない）と似た考え方が使えるのかもしれない。 朝永振一郎先生の「物理学とは何だろうか」（岩波新書）に、生物学者向けの解説が詳しく記されている。
それらの遺伝子に対応する値はどのような分布を取るだろうか。もし遺伝子を選んだ基準と、元のデータの処理条件に全く関連がないならば、それはランダムに遺伝子を選んだのと同じことになる。 ゆえにその場合、分布の平均値（標本平均）は母分布の平均値とほとんど同じ値になる。
一方、遺伝子を選んだ基準と元のデータの処理条件に強い関連がある場合がある。 例えば、「低酸素処理を行った細胞／対照の細胞」 が元のデータ、選んだ基準が 「HIF で強く誘導される遺伝子 50 個」だとする。 その場合、得られる分布の平均値は母分布の平均とは異なる、高い値に偏る。 ある基準で選んだ遺伝子リストで、得られた平均値が母分布の平均と大きく異なり偏っていれば、その遺伝子リストを選んだ基準は、元のデータの処理条件と何らかの関連があると考えることができる。遺伝子リストをたくさん用意してそれぞれ計算すれば、元のデータの処理条件と関連がある物事はどんなものか、偏りの度合いを用いて判断できる。
しかし、母分布自体がとても分散が大きいなら、その偏りはそのせいかもしれない。そこで平均値の違いを母分布の標準偏差で補正する。 また遺伝子リストに含まれる遺伝子数が少ないと偶然の偏りが起きやすいので、それも補正する。それらは単純な割り算、かけ算で行える。この操作によって、「母分布における平均と、遺伝子リストのメンバーでの分布における平均が一致する」と言う条件では、Z score は平均０，標準偏差１の標準正規分布に従った分布を示すようになる。
言い換えると、「もし遺伝子リストのメンバーをランダムに選んで Zscore を計算することを何回も繰り返すと、求められた Zscore の分布は、平均値で一番山が高く左右になだらかに下がっていく正規分布になる。Zscore の定義がうまく作ってあるので、ただの正規分布ではなく平均＝０、標準偏差＝１の標準正規分布になる。 この場合 Z score がある値よりも偶然大きくなる確率は標準正規分布の数値から求められ、例えば 1.96 よりも大きくなる確率は 2.5% （-1.96 より小さくなる確率も 2.5% で、合わせると 5%）というように決まっている。」ということになる。 これを逆に解釈して、「ある遺伝子リストのメンバーから計算した Z scoreがとても大きい。遺伝子がランダムに選ばれた場合、このようなことが起きる確率は 1% にも満たない。 ということは、この遺伝子リストに含まれる遺伝子は、元のデータの処理条件と強い関連がある（そのせいで母分布からランダムに選ばれたものではなくなっている）のではないか」というように考えることができる。
上に書いた考えにはかなり重大な問題がある。遺伝子リストを例えば100通り用意したとする。それら全部について Zscore を計算する。そうすると、どれか一つか二つくらいは偶然に Zscore の絶対値が 2.6 を超えてもおかしくないことになる。だから「Zscore が大きくなるリストをたくさんの候補から選ぶ」のはよいとしても、それは偶然のせいではないことを確かめないといけない。Zscore が十分大きければ(3.29 なら p-value = 0.001 1000回に1回）問題ないだろうし、生物学的な事実や実験結果と一致していれば Zscore がやや低くても支持されるのではないか。また相関を見るなどの他の方法と組み合わせることも考えられる。
遺伝子発現の分析には相関係数がよく使われ、その有効性が実証されている。相関係数が高くなるには、二つのセットの間に全体を覆うような直線的な関係がある必要がある。それがあることは、関連があることと強く結びつく。しかし Zscore では偏りがあればスコアが高くなるので、Zscore が高くても相関係数が０になることはほとんどである。その分、本当に関連があると言ってよいか、気をつけないといけない。
B ｜ *
｜ * B と A に何も関係がない場合、B の平均値の期待値は０ しかし
｜ *
｜ *
-―――――――十―――――――― A が正なら、B も正 しかも相関がある
｜ A
｜
｜
｜
B ｜ *
｜ * B と A に何も関係がない場合、B の平均値の期待値は０ しかし
｜ *
｜ *
-―――――――十―――――――― A が正なら、B も正 しかし相関はない
｜ A
｜
｜
｜
下の場合でも、偏りはあるので Zscore は高くなる。
相関の場合でも擬似相関という問題があることが知られている。遺伝子の場合は生物学的に遺伝子の役割、上位下位が推定できる場合が多いので、擬似相関はあまり問題にはならない（本当は違うのかもしれないが）。Zscore の場合もそうだろう。いずれにせよ、生物学的な事実と組み合わせないといけない。
Z score を計算する数式は Kim らの論文の Methods に書かれている。 エクセルではこんな感じになる
=(E4-$B$4)*SQRT(E6)/$B$5 この場合、$B$4 に母分布の平均、$B$5 に母分布の標準偏差、E4 に遺伝子リストから得られた値の平均、E6 に遺伝子リストに含まれる遺伝子の数が入っている。
Z score は p-value に換算できる。そのやり方も Kim らの論文、門田幸二 博士による解説に書かれている。Z score は標準正規分布に従うことを利用する（Z score の Z という文字の元になっている Z 変換というのは、標準正規分布にあてはめるための変換）。
=(1-NORMDIST(ABS(C7),0,1,TRUE))*2 C7 に Z score
Z score が 2.6 を超すと、p-value は 0.01 より小さくなる。偶然によって平均値に偏りが生じた確率が 1% より小さいと解釈できる。よく生物学のデータでは ｔ 検定を行う。ｔ 表というのがあり、「自由度が 10 で基準が 1% なら 3.169」などと書いてある。その表で「自由度∞」というところの 1% の値が2.576 で、同じ値である。t 検定でサンプル数（自由度）が増えると正規分布に近づくので一致する。t 分布の自由度 30、1% では 2.75 で、正規分布の 2.576 とたいして変わらない。PAGE 法では遺伝子リストのメンバーが 10 あればあまり問題がないと論文には書かれている。Z score の計算法と t 値の計算法は基本的によく似ている。標本数がとても少ないときは t 値を使うことががふさわしくなり、その値が示す分布は標準正規分布ではなく、もっと拡がった t分布になる。 マイクロアレイでは多数の遺伝子の値を用いることができるのでサンプル数は十分多くなる。また母分布の標準偏差をデータから求めることができる。そのため Z score で問題がない。
この方法では遺伝子リストとして arbitrary 恣意的に選んだリストを用いる。リストされる遺伝子を選ぶ基準として、マイクロアレイ実験などの結果を用いることができる。そのときに、測定ノイズなどの影響で何の関係もない遺伝子がリストに紛れ込むことがありうる。しかしそれらの関係ない遺伝子から得られる発現レシオの値はほとんどの場合０に近い値になることが期待できる。そのため、それらの関係ない遺伝子による平均値の偏りはおきにくい。この仕組みは情報をうまく取り出すことに寄与できる。
「寺田寅彦随筆集(第二巻)」（小宮豊隆編 岩波文庫）の171ページから、「比較言語学における統計的研究法の可能性について」 という文章がある（昭和三年）。
「偶然による（統計的に期待されるべき）符合率と、実際の（実測したデータによる）符号率の比のみが意味を持つ、ここではそれを問題にした」と書かれている。 この考え方は GSEA 法の考え方と似ている。
言語 A と言語 B について、語彙に含まれる子音に注目する（この目の付け所がよい）。 子音が n 種類あって、すべて出現頻度が同じ (1/n) とする（近似： 近似しても問題ない部分と、精密に考えなければならない部分をうまく分けている）。平均値で近似している、子音の種類の違いを無視して、気体を構成する分子一つ一つのようにそれぞれを区別せずに扱うと言うことである。この場合はそれでもほとんど問題がない。
語彙を、含まれる子音の数によって分類（グループ1, 2, 3,...）する。 それぞれのグループについて、A, B 両国語に現れる確率を a1, a2,... b1, b2,... と計算する。分数の形で表現した方がよいかもしれない（グループに含まれる数／すべての語彙）。 この段階では、A と B を全く別々に扱う。 いまではネットと計算機があるので、文章と単語をたくさん用意してカウントできる。できるだけ多数の文章について調べるとよいのだろう（データ全体の性質を抽出する）。
言語 A で、グループ i から一つ単語を取り出す。
それと同義の単語が言語 B にあるとして、それが A と同じくグループ i に含まれ、さらに含まれる子音 (i個) も同じものである確率は
になる。この値は１よりもかなり小さいと考えられる。
さらにシノニム（synonym 同義語）がそれぞれ Si 個ずつあるとする。その分をかけ算すると、
になる。単語の長さは限られているので i はせいぜい最大で15くらいだろう。そういう長い単語は少ない。これを i = 1, 2, 3, 4,... と計算して全部足したものが、「言語 A と言語 B について多数の語彙について比較した際に、偶然に意味、子音の数と種類が一致する確率」になる。全部足さなくても、それぞれの値を並べて用いてもよいかもしれない。
ここで、言語 A と言語 B の多数の語彙から、意味、子音の数と種類が一致しているものを実際に取り出してみる。i = 1 のものでは調べた語彙 N1 個で a1 個が一致しているなら、「実測値」は a1/N1 になる。i = 1, 2, 3, 4.. と実測値を計算する。この場合、語彙を選ぶときに「食物に関する語彙」「天候に関する語彙」などのように恣意的に何らかの選択基準をつけて選択することもできる。
最初に得た、偶然に一致する確率値（またはその数列）を、実際のデータから統計的に得られた実測値（または数列）と比較することで、言語 A と言語 B の語彙に関連があるかどうかを判断する基準になる。「偶然に一致する確率」と「実際に一致した確率」を対にして比較する。
「統計的方法の長所として、関係要素の数が多くて、それら相互の交渉が複雑であればあるほど、かえってこの方法の妥当性がよくなるという点である」と、寺田先生は書かれている。マイクロアレイデータの分析に、こういうたぐいの方法が有用であることが予言されているようである。
寺田先生は語彙に含まれる子音の数と種類に注目している。現代生物学なら、遺伝子のプロモーターやタンパク質のアミノ酸配列に含まれる特定のモチーフを、語彙に含まれる子音のように解釈できるかもしれない。
PAGE 法の場合は、fold change の平均値を用い、Zscore が正規分布すると想定しているのでパラメトリックな統計であるということらしい。Zscore が正規分布すると考えてよいことは彼らの論文でも検定され確かめられている。遺伝子セットをランダムに（偶然に任せて）選んだ場合、 fold change の平均が母分布の平均と同じ値（〜０）になることが期待される。それを、何らかの基準で恣意的に選択した遺伝子セットから統計的に得られた値の平均値、分布と比較することで、データと選択基準に関連があるかどうかを判断できる。
プロモーターのシスエレメント分析に、これと似た感じの方法を用いることで優れた結果が得られている。 山本教授らによって植物遺伝子プロモーターに存在するシスエレメントの高精度な予測法が開発され、実験的な検証にも耐えるものであることが示された。すばらしい成果である。
Prediction of transcriptional regulatory elements for plant hormone responses based on microarray data. Yamamoto YY, Yoshioka Y, Hyakumachi M, Maruyama K, Yamaguchi-Shinozaki K, Tokizawa M, Koyama H. BMC Plant Biol. 2011 Feb 24;11:39. PMID: 21349196
この方法では、RAR (Relative Appearance Ratio) という値を基本にする。
8個の塩基をまとめてオクタマーとして考える。塩基は4種類あるので、4^8 = 65535 種類のオクタマーがある。単純に平均化して考えると、あるオクタマーがゲノムに出現する頻度は 1 回 / 65535 塩基になる。寺田先生の言語解析の場合は、子音の出現頻度をそのように単純化している。 あるプロモーターに、特定のオクタマーが 1 / 65535 よりも高い頻度で出現していれば、それに意味があると考えることもあるかもしれない。しかしそれではうまくいかない。
シロイヌナズナの場合、ゲノム配列がすべて決定されている。山本教授らは、ゲノムから 14960 個の遺伝子のプロモーター配列を取り出した。それらの配列で、それぞれのオクタマーがいくつずつ出現するかをカウントし、それを基本としている。実際のプロモーターではオクタマーごとに出現頻度はかなり違い、それを平均値で近似してはいけないらしい。 GCGCCCCC のような GC が並んだオクタマーが出現する頻度は低い。AAAAAAAA TTTTTTTT ATATATAT のようなオクタマーが出現する頻度は高い。
オクタマーは、65536 個の区画と考えることができる。それらの区画にそれぞれ出現頻度を数字として入れる。その結果はヒストグラム、分布として表現できる。分布が連続的でなく、いくつかの区間に分けられている場合、それらの区間に一様に分布するとエントロピーが最大になる。オクタマーごとに出現頻度が大きく違うということは、プロモーター部分の塩基配列のエントロピーはとても低いと言うことになる。塩基配列からエントロピーを計算することで有用な情報が得られることは多いらしい。「Sequence logo」という塩基配列表示法に関する説明 http://en.wikipedia.org/wiki/Sequence_logo にエントロピーが出てくる。ATGC が一様な頻度で出現するならエントロピーは高い。特定の塩基が出る頻度が高いならエントロピーは低い。エントロピーが低いときに、その原因になっている塩基を大きい文字で表示する。「エンリッチメント」というのも「分布の仕方が変化する」ということに置き換えられるので、エントロピーの計算法を用いてスコアをつけられるかもしれない。「母分布と同じ分布をしている」場合にスコアが最大または最小になればよい。しかしそれで特に何もよいことがなければ意味がない。またスコアから統計的に p-value のような値が計算できないと結果を信頼しにくくなるだろう。
山本先生の 2011 年の論文を解読してみる。
この方法でも遺伝子リストを使う。Table3 にあるように、公開されているアレイデータを元にしている。「ABA 処理による遺伝子発現を調べたデータ」等を用いている。それらのデータから、「発現が処理によって対照の3倍以上に増加した遺伝子のリスト」を作成している。それぞれ 67〜362個の遺伝子(プロモーター）が選ばれている。この場合は遺伝子セットというよりも「プロモーターセット」になる。
山本先生の論文が高く評価されているということは、「公開されているアレイデータを元に遺伝子セットを作る」ことの有用性・問題のなさを示している。
それぞれのオクタマーについて、RAR値を計算する。計算式は、
（あるオクタマーがプロモーターセットに出現する回数 / プロモーターセットに含まれるプロモーターの数）／（あるオクタマーがすべてのプロモーターに出現する回数 / 今回用意したすべてのプロモーターの数 = 14960）
になる。プロモーターセットをランダムに偶然に任せて選んだ場合、この値はどのオクタマーに対しても１になることが期待される（プロモーターセットによって濃縮されることがない）。4^8 = 65535 種類のオクタマーがあるので、RAR値も65535個計算される。「発現が処理によって対照の3倍以上に増加した遺伝子のリスト」で計算すると、「発現上昇に役割を果たしているオクタマー」が存在する可能性が高いことにより、それらのオクタマーに対するRAR値はそうでないオクタマーよりも高い値を取ることが期待される。
しかしこれだけでは不十分で問題がある。65535 種類のオクタマーには、 今回用意したすべてのプロモーターに出現する回数が非常に小さいものがある。その場合、RAR値の分母が非常に小さい値になる。そういうオクタマーがプロモーターセットに一回でも出現すると、そのオクタマーに対するRAR値が異常に大きい値になる。
そのことに対処するために「フィッシャーの Exact Test」を用いている。これはカイ二乗検定と同様に「実測値による比率が期待される比率と異なるかどうか」を検定するために使うことができる。「フィッシャーの Exact Test」は場合の数を元にしている（直接確率法）。標本数が少ないときは、この方法が正しい値を出す。計算するために、プロモーターセット、オクタマーの組み合わせごとに２ｘ２の分割表をつくる。 参考にした資料 ttp://aoki2.si.gunma-u.ac.jp/exact/exact.html
| あるオクタマーが存在する | すべて | |
| あるプロモーターセット | a | b |
| 全体 | c | d |
| 計 | a+c | b+d |
RAR値は、この表の値を使うと (a/b) / (c/d) になる。
a/b と c/d に差がないかどうかを検定する。差がないなら、全体での割合 c/d とプロモーターセットでの割合 a/b を区別することができないと言うことになる。それは、そのオクタマーについては、プロモーターセットにより濃縮されることがない (enrichment がおきない) ということと同じである。分割表ごとに、その状態が生じる確率を計算できる。それらを足し合わせて p-value が得られる。p-value が大きいときは a/b と c/d は差がないことが支持される。全体でのオクタマー存在数 c の値がとても小さい場合（極端な場合）は、足し合わせる場合の数が多くなり p-value が大きくなりやすい。p-value は LOD (log of odds) score に換算できる。LOD = -log10(p-value) p=0.05 で LOD=1.3 LODスコアが小さい＝p-value が大きい。その場合RAR値が大きくても（a/b が c/d よりもずっと大きい）、統計的には a/b と c/d に差がないので、RAR値は０ということにする。その処理をした値を RARf value と名付けた。
まとめると、「プロモーターセット、オクタマーをそれぞれ指定すると、(RAR, RARf, LOD score) という3種類の値が帰ってくる」そういう数値テーブル、関数を作成したということになる。シロイヌナズナの全遺伝子のプロモーター配列は TAIR からファイルとして誰でも入手できるものがある（このデータでは正確な転写開始点がわからないが）。 ftp://ftp.arabidopsis.org/home/tair/Sequences/blast_datasets/TAIR10_blastsets/upstream_sequences/ TAIR10_upstream_1000_20101104 というものなどがある。
つぎに、分析対象になる遺伝子のプロモーター配列を用意する。これを端から8塩基ずつ取り出しオクタマーとする。「プロモーターセット、オクタマーを指定すると、(RAR, RARf, LOD score) という3種類の値が帰ってくる」関数を呼び出して、それらの値を得る。上流から下流まで1塩基ずつずらしながら（スキャン）オクタマーを作っていき、その都度関数を呼び出す。得られる値を Figure 2 のようにプロットする。これによってシスエレメント予測ができる。プロモーターセットを変えて分析しそれらの結果を比較することもできる。
山本先生は、RAR値を用いる方法とは別に、「プロモーターで機能を果たしている重要な配列は、プロモーター内で出現する部位に偏りがある」ことに着目し、プロモーターの機能配列を網羅的に検出する方法をすでに開発されている。これもすばらしい成果である。 http://www1.gifu-u.ac.jp/~yyy/pdf/12Hieno_ppdbNewsLetter4_3.pdf Identification of plant promoter constituents by analysis of local distribution of short sequences. Yamamoto YY, Ichida H, Matsui M, Obokata J, Sakurai T, Satou M, Seki M, Shinozaki K, Abe T. BMC Genomics. 2007 Mar 8;8:67. PMID: 17346352 この成果も、分布（この場合はプロモーター内のどこに出現するか）と、その仕方の変化、偏りに着目している。重要な配列は転写開始点の近くに出現する頻度が高い（エンリッチされている）ことがグラフで示されている。エンリッチメント分析の一つと考えてもよいのかもしれない。この手法によって見いだされた配列と、RAR 値を用いて見いだされたシスエレメントが比較対照されている。二つの手法はどちらもエンリッチメント（分布の偏り）に着目しているが独立した異なるものなので、双方に共通して検出されたエレメントは信頼性が高いだろう。
山本先生の成果から、「同じような発現様式をする複数の遺伝子があっても、それら全部に共通して存在する、保存されたシスエレメント配列が見つかることはあまりない。そういうものが見つかる方が珍しい」ということが推測される。山本先生による解説 http://www1.gifu-u.ac.jp/~yyy/pdf/12Hieno_ppdbNewsLetter4_3.pdf にそういうことが書かれている。 「コンセンサス配列でも、転写制御に関わっているとは限らない」ということもあるそうである。
また、化学と生物（日本農芸化学会誌） 58(5): 267-268 (2020) では、「データマイニングによる転写因子の結合配列環境と発現応答の関連性の解析プロモーター配列から発現変動遺伝子を予測」という兒島 孝明、岡 大椰、中野 秀雄 各先生によるすばらしい解説が掲載されていた。山本先生の成果のような、プロモーターに存在する転写因子結合モチーフの配列、出現頻度、出現位置の分析は、そのプロモーターの転写様式を推定するために役立つ。しかしそれらだけでは十分ではないことがわかってきた。さらにプロモーターの塩基配列から DNA 構造を計算してその結果を結合モチーフと組み合わせることが有効であると紹介されている。 DNA 構造の計算も、兒島先生らが紹介されているパラメーター以外にもまだ工夫の余地があるだろう。
http://www.riken.jp/pr/press/2016/20160510_1/ では、理化学研究所の川上英良特別研究員らのグループによって、従来のGene Set Enrichment解析に確率的な制御関係を考慮した wPGSA法（weighted Parametric Gene Set Analysis） という方法が開発されたことが紹介されている。
http://wpgsa.org/ https://qiita.com/vanillin99/items/7adf8ce3263ede233d38 「遺伝子発現から転写因子の活性化を予測するエンリッチメント解析です」と書かれている。
「ある転写因子の配下として制御を受ける遺伝子のグループ」は、共通した性質を持つ遺伝子のグループ、セットと考えることができる。公開されている ChIP 実験のデータを大量に集めて分析しやすいフォーマットにまとめ（これだけでもものすごく大変で価値が高いことである）、 新たに開発された wPGSA法によって遺伝子発現データから転写因子の活性を極めて高い精度で予測することに成功したと紹介されている。ChIP 実験からのデータを、普通の RNAseq やマイクロアレイのデータとうまく組み合わせることで、より高度なデータ分析を実現している。複数のデータをうまく組み合わせることはとても重要で工夫の余地がたくさんあるのかもしれない。 例えば、Broad distribution spectrum from Gaussian to power law appears in stochastic variations in RNA-seq data. Awazu A, Tanabe T, Kamitani M, Tezuka A, Nagano AJ. Sci Rep. 2018 May 29;8(1):8339. doi: 10.1038/s41598-018-26735-4. PMID: 29844539 粟津先生、永野先生らのグループによるこの論文では、標準的な条件で生育しているシロイヌナズナの遺伝子発現データを時系列で、各時間ごとに約 20 サンプルを RNA-seq によって分析している。 発現量のヒストグラムは、遺伝子ごとに特徴がある。同じ条件で育成していても発現量に大きな揺らぎ・外れ値が観測される遺伝子もあれば、そうでない遺伝子もある。似た特徴を持つ遺伝子をクラスタリングでグループにすることで平均を取れるようにしてノイズを打ち消し、分布を表現する関数を推定している。どのようなメカニズムで分布関数が成立しているかについてモデルが立てられている。RNA-seq の結果を使いやすくまとめたデータが論文に添付されている。このようなデータを他のデータを組み合わせることで今までにないことがわかるかもしれない。
また、化学と生物（日本農芸化学会誌） 58(5): 267-268 (2020) では、「データマイニングによる転写因子の結合配列環境と発現応答の関連性の解析プロモーター配列から発現変動遺伝子を予測」という兒島 孝明、岡 大椰、中野 秀雄 各先生によるすばらしい解説が掲載されていた。プロモーターに存在する転写因子結合モチーフの配列、出現頻度、出現位置は、そのプロモーターの転写様式を推定するために役立つ。しかしそれらだけでは十分ではないことがわかってきた。さらにプロモーターの塩基配列から DNA 構造を計算してその結果を結合モチーフと組み合わせることが有効であると紹介されている。
「機械学習によってバラバラな細胞たちをパズルのように組み立てる〜1細胞計測データからの遺伝子発現マップの高精度予測〜」 広島大学のプレスリリース https://www.hiroshima-u.ac.jp/news/65249 １細胞RNAシーケンシング(RNA-seq)法で計測された遺伝子発現データと、位置情報を持つ ISH (in situ hybridization) データとを組み合わせ統合することで、組織を構成する細胞一つ一つの位置と、そこに存在する細胞の遺伝子発現を再構成することに成功した。
Ishimori 先生による、ChIP-seq データなどを用いた転写因子結合部位の予測に関する解説論文 Commentary が発表されている。 Transcription Factor Binding Site Prediction: Finding the Point from Many Data Motoyuki Ishimori Plant and Cell Physiology, Volume 63, Issue 10, October 2022, Pages 1324–1325, https://doi.org/10.1093/pcp/pcac128
wPGSA法
Weighted enrichment method for prediction of transcription regulators from transcriptome and global chromatin immunoprecipitation data. Kawakami E, Nakaoka S, Ohta T, Kitano H. Nucleic Acids Res. 2016 Apr 30. pii: gkw355. PMID: 27131787
この論文を解読してみる。
代謝、細胞周期、分化、ストレス応答などのほとんどの細胞活動は遺伝子発現によって駆動される。様々な環境要因の影響によって、細胞は何千もの遺伝子の発現を適切に調節し、機能を持つタンパク質を必要な量合成できるようにする。遺伝子発現を調節する主要な因子である TFs（transcription factors） は、それぞれ特異的な DNA 配列に結合し、RNA 転写の速度を決定する。 それだけでなく、クロマチンを構成するタンパク質、DNA を修飾するタンパク質、TFs と共同して働く co-regulators も、TFs が標的 DNA 配列に結合する効率に影響を与える。 これらの制御因子の組み合わせは、込み入って精巧な遺伝子発現の調節を可能にする。また内部、外部からの様々な刺激に対応した細胞の恒常性維持を可能にする。最近、動物の成熟した細胞は、いくつかの TFs を導入することで多能性を持つ状態に再プログラムできることが示された。ゆえに、TFs とその標的を同定することは、正常、または病的な状態の細胞の制御機構を理解する上で重要な第一ステップとなる。また細胞の分化を人為的に制御するためにも重要である。
TFs に関して膨大な研究があるにもかかわらず、細胞内のそれぞれの TFs の活性状態を網羅的に確かめることは未だに難しい。それは TFs の制御機構がきわめて複雑な階層構造であるためである。多種多様な翻訳後修飾、リン酸化、糖付加、アセチル化、ユビキチン化、SUMO化などが TFs の機能に大切なことが知られている。さらにその上、TFs の標的配列への結合は、クロマチンの状態、他の TFs や co-regulators との組み合わせによる相互作用によって影響される。 これらの TFs とクロマチンの状態は抗体を用いた ChIP クロマチン免疫沈降実験で分析することができるが、すべての TFs に対してそれを行うのは非現実的である。 しかし、我々は今ではマイクロアレイや RNAseq によってゲノム全体の mRNA 発現状態を知ることができる。ゆえに、それらのデータを用いて、どの TFs がどのような mRNA 発現に役割を果たしているかを分析、予測するのは最も有望な方法である。
どの TFs がある与えられた網羅的発現データに深く関与しているかを予測するために、我々は、まず TFs とそれらの標的遺伝子群の関係性を再構築する必要がある。この関係性は、GRN (gene regulatory network) と呼ばれる。最近、それぞれの TFs が結合する標的配列を元にして、GRN を推定する方法が提案された。その論文の著者らは、プロモーターの近傍に位置する結合配列の種類と数によって、それぞれの遺伝子の発現量を説明するリニアモデル、線形回帰によるモデルを構築している。 他の方法としては、クロマチン免疫沈降実験の結果を用い、TFs とその標的遺伝子を結びつける分析も行われている。結合配列（モチーフ）に基づいた方法と比較して、クロマチン免疫沈降実験に基づいたモデルは、ある細胞の状態における真の結合部位に関する情報を考慮しているという点で優れている。しかし同時に、クロマチン免疫沈降実験に基づいた方法は、実験を行った条件、細胞の種類の違いによって大きく影響を受けるという問題点もある。
我々はこの論文で、数多くの細胞種と条件で得られた、多数のマウスのクロマチン免疫沈降実験のデータを集めた。それらのデータに基づいて、TFs とそれらの標的遺伝子群の結びつきを予測する方法を提案する。この方法では、GSA (gene set analysis) において、条件に依存した TFs と標的遺伝子群の関連を考慮するために、重みを付けた t-test (weighted t-test) を用いている。この方法で予測されたマウスの TFs は、細胞の種類や実験条件にかかわらず、生物学的な実験結果ときわめてよく一致していた。具体的な事例も示す。
ChIP 実験とは、Chromatin Immunoprecipitation を略した名称で、クロマチン免疫沈降実験と呼ばれる。 この実験は、ある転写因子に対する特異抗体を用いて、その転写因子と結合している DNA 配列を同定するために行われる。分析したい TF に対して特異性が高く、高い親和性で結合する抗体が必要になる。評判がよい抗体が売られているなら買えばよいが、そうでなければ自分で用意することになり、難しい。または目的の転写因子のタンパクに、抗体と結合するタグ配列を付加した遺伝子を導入した組み替え植物や動物を作らないといけない。クロマチンの単離や抗体との結合、洗浄、溶出などの操作もうまく行わないと意味のある結果にならないだろう。わたしはやったことがない。
手順： 分析対象になる細胞を用意する。まず細胞全体をホルムアルデヒド処理して、クロマチンと転写因子とDNAの結合を固定する。細胞を溶解し、超音波処理でクロマチンを溶解性が高く分析に適した大きさに分解する。抗体を加えて結合させ、Protein G セファロースをさらに加えて固定相にトラップする。洗浄によって抗体と結合しなかったクロマチン画分を洗い流す。 抗体と特異的に結合することで固定相にトラップされたクロマチンに含まれる DNA を、固定相から溶出し、分析に用いる。 この論文では、溶出した DNA を次世代シーケンサで分析したデータ (ChIP-seq)、マイクロアレイで分析したデータ (ChIP-array) を用いている。
ChIP-seq のデータは、GEO database から得られるものを利用している。SRA toolkit （RNAseq などのデータを処理するためのプログラム群）を用いて、FASTQ フォーマットに直している。データの質のチェックには FastQC というプログラムを用いている。読まれた塩基配列をゲノム配列に align するには Bowtie2 プログラムを用いている。
ChIP-array のデータは、Ringo という Bioconductor のパッケージで処理している。ChIP によって enrich 濃縮された DNA 領域を決定するためにも、Ringo を用いている。Affimetrix のアレイでは、rMAT というパッケージを用いて、同様なことを行っている。
ChIP 実験を行うことによって、ゲノム全体から enrich 濃縮された DNA 配列（ピーク）を見いだすことができる。それらのピークを検出したデータはいくつかのフォーマットで出力されている。それらを全部 bed というフォーマットに直している。それらのピークがゲノムのどの領域に相当するかを、UCSC LiftOver Tool や ChIPpeakAnno というソフトウェアで検出している。特に TSSs transcription start sites の近傍に注目して分析している。
用いたデータの数などについて、RESULTS に記載されている。マウスでは ChIP のデータが豊富に入手できるらしい。ChIP では転写因子に対する特異性と結合能が高い性質のよい抗体が必要になる。 マウスは研究が進んでいてそういう抗体が多いのだろう。また次世代シーケンサが普及したことによってデータが得られやすくなっている。454 個の転写因子について、8,578,738 857万8738 個の結合データを得ることができた。 実験数で数えると、273 の ChIP-chip 実験と、2924 の ChIP-seq 実験になる。足し合わせると 3197 で、8578738 を 3197 で割り算すると答えは 2683.37 になる。 単に割り算しただけでは、「一つの転写因子がゲノム上に平均的にもっている結合部位の数」はわからない。
まず、ChIP の結合データに基づいて、遺伝子発現の調節を確率の観点からとらえることを試みた。「発現調節確率 probability of regulation 」という値を考えている。これは、対象にしている転写因子が、核内の標的配列に対して結合している確率である。これを p とする。 発現調節確率 p は、最尤推定法 maximum-likelihood estimation によって推定した。
最尤推定法では、尤度関数 likelihood function を考えないといけない。この場合、二項分布 binomial distribution を基本にして考えている。 転写因子は、標的配列に対して「結合する」「結合しない」の二つの状態を取り得る。だから二項分布の考え方を使う。 この場合の尤度関数 likelihood function は、発現調節確率 p と、結合していると判定された回数 k の、二つの引数を持つ関数として定義されている。さらにデータの総数 n も式に入っている。0 <= k <= n ということになる。
この式を解読してみる。この式は、二項分布 binomial distribution における確率を求める式に基づいている。 「確率 p で起きることがわかっている事象が、n 回試行（コインを投げる、実験を行う）を行った際に、ちょうど k 回起きる確率を表す式」
の方は、n 回の試行（n 回実験を行う・この論文の場合は n 個の実験データを用意して）で、ちょうど k 回事象が起きて（結合して）、(n - k) 回起きない（結合しない）という結果になる確率を表している。
具体的な例を作ると、3 回の実験で 2 回結合する場合、結合することを 1、しないことを 0 として、0, 1, 1 という場合、1, 0, 1 という場合、1, 1, 0 という場合の 3 通りがある。
それぞれの場合について、その通りになる確率は で、共通した同じ値になる。
この値をすべての場合について足し合わせると、求める確率になる。都合のよいことにすべての場合で同じ なので、この値に場合の数をかけ算するとよいことになる。
だから、場合の数を求める式が必要になる。
の部分は、高校数学で勉強する、n 個の粒子から k 個の、お互いに区別が付かない粒子を選び出す場合の数である。
n 個の粒子が、すべて区別が付くとして、順番に並べる場合の数は分子にある n! になる。
しかしそこから k 個を選ぶ場合、選ばれずに残る粒子には順番が付かないので、その分の場合の数を割り算で減らさないといけない。
それが分母の (n - k)! になる。さらに選ばれる k 個の粒子もお互いに区別が付かないなら、その分の場合の数も割り算で減らさないといけない。それが分母の k! になる。
これで、場合の数と、それぞれの確率が求められたので、両者をかけ算（＝それぞれの確率を場合の数だけ足し合わせる）すると、「確率 p で起きることがわかっている事象（転写因子が標的配列に結合する）が、n 個の実験データにおいて、ちょうど k 回結合している確率を表す式」になる。
ここで、「n 個の実験データにおいて、何回結合しているか」は、実験データから直接求めることができる。 だから n と k はデータから求められた値として決まっていることになる。 変化させることができるのは p の値だけで、この値を動かして、一番上の式の値が大きくなるところを選ぶ。それによって確率 p を実験データから推定できる。
上の式のままでは階乗が入っていて、計算がやりにくい。 そこで対数をとった「対数尤度関数」を作ると、次の式のようになる。log(AB) = log(A) + log(B) なので、3 個の項の和にできて式がわかりやすくなる。
一番目の項の値は n と k だけで決まるので、定数として K で表している。
対数尤度関数を p について微分（p 以外は固定するので偏微分）して、その値（元の尤度関数の傾き）が 0 になる p を求めることで、確率 p を推定できる。
基本的な公式 から、log p を微分すると 1 / p になるので、書いてあるように
ということになる。
でよいことになる。単純な式で p を決めることができた。
この値を
という記号で書いている。この値は、遺伝子発現の変化（対数で表している Log Fold Change = LFC）に、転写因子それぞれの結合能による重みを付けて補正するために用いられる。
つぎにこの論文では、ある転写因子の配下にある遺伝子のセット (gene set) の遺伝子発現の変化（対数で表している Log Fold Change = LFC）と標準偏差を、遺伝子全体の発現量の平均と標準偏差に対して t スコアを計算することで比較している。
この論文での t スコア（二群のデータから計算する t スコア）は、以下の式に則って計算される。
は グループの発現変化 (LFC) の平均
は転写因子 T の支配下にある遺伝子のセットの発現変化 (LFC) の平均
は、すべての遺伝子の発現変化 (LFC) の平均
n は グループのメンバーの数
は グループの発現量の標準偏差 SD
の部分は、分子の方は二つの群を合わせた分散の推定値の平方根＝SD に相当する。データが一群の際の SD に対応している。
分母は単に
で、
に対しては
をかけ算することになる。これはデータが一群の場合に t スコアを計算する式に出てくるものに相当する。
一番簡単な、一群のデータに対して t スコアを計算する場合の式は
t スコア＝ N はセットに含まれるデータの数
となる。この論文の場合データは二群だが、片方はデータ全体でデータの数がとても大きいということから、それを普通の t スコアの計算法で式に取り入れても、結局遺伝子セットが含む遺伝子の数 n しか意味を持たなくなる。
で、n2 は 2 万で n1 は 30 なら、
とほとんど変わらない。
だからデータの数としては、遺伝子セットが含む遺伝子の数 n がこの場合ふさわしい。それで一番簡単な一群のデータの場合の式と似たものになっている。
普通に考えると、; 転写因子 T の支配下にある遺伝子のセットの発現変化 (LFC) の平均 は、例えばマイクロアレイの結果を RMA法などで処理した値をそのまま使う。PAGE 法ではそうするようになっている。
しかしこの論文ではそこを工夫して、上の方で求めた発現調節確率 p を使った重み付けを行い、優れた分析力を発揮できるようになっている。この分野の論文では、新しい工夫、考え方が論文のどこかに提案されていることが、高い評価を得るために必要なのだろう。
; 転写因子 T の支配下にある遺伝子のセットの発現変化 (LFC) の平均を求めるために、以下のような式で計算している。
この式の分子は、データとして得られた、遺伝子セットに含まれるそれぞれの遺伝子の Log Fold Change (LFC) = に、その遺伝子に対する転写因子 T が結合し、発現を調節している確率
を重みとしてかけ算している。
これによって、発現調節確率が高い場合は、その配下遺伝子の LFC の値が大きく寄与することになり、発現調節確率が低い場合は、いくらデータから求められた LFC が大きい値でも、寄与が小さくなることになる。
分母ではすべての を足し合わせている。これは規格化のために行っていることのようで、これがないと
の値に大きいものが多い場合と、0 に近い値が多い場合で
の値の違いがとても大きくなってしまう。
すべての
を足し合わせた値で割り算することで、発現調節確率それぞれの値の大きさよりも、それぞれの発現調節確率の値の比率が重要になるようになっている。
他の分野でも、a と b の差を考えるときに、単に a - b とせず、a + b で割って と規格化することがある。
このようにも考えられる。
Z はすべての場合の確率を足し合わせた値で、それぞれの確率をこの値で割り算すると、相対的な頻度・確率になる。それに、その場合に実現される値（この場合 Ei）をかけ算すると、期待値になる。平均値は期待値の特別な場合と考えてもよい。すべての期待値を全部足して平均値としている。
とても簡単な例を構成すると、
p1 = a, p2 = b, E1 = X, E2 = Y とすると、求める重み付け平均値 Et = (aX + bY) / (a + b)
ここで a, b がどちらも C 倍になったとすると、分子の方は C 倍に変化する。分母の a + b で割り算しておくと、こちらも C 倍になるので、重み付け平均値は変化しないことになる。 a, b がどちらも C 倍になっても a と b の比率は変化しないので、重み（発現調節確率）の比率が重要になるようになっている。
発現調節確率は k 結合が観測された回数 / n 実験データの数 で推定されるので、n = 1 や 2 では一応値は出るが、分析には使いにくい。Fig. 3A では、例として n = 20 になっている。
以前はマウスや人以外の生物ではデータが少なかったので、発現調節確率を求めることはできなかった。しかし最近多くのの生物種でも大量のデータが取得され公開されるようになった。
PlantPAN 3.0 http://plantpan.itps.ncku.edu.tw/
Transcription Factor Binding Site Prediction: Finding the Point from Many Data Motoyuki Ishimori Plant and Cell Physiology, Volume 63, Issue 10, October 2022, Pages 1324–1325, https://doi.org/10.1093/pcp/pcac128
Ishimori 先生の解説論文 Commentary において PlantPAN 3.0 が紹介されていた。
PlantPAN3.0: a new and updated resource for reconstructing transcriptional regulatory networks from ChIP-seq experiments in plants. Chow CN, Lee TY, Hung YC, Li GZ, Tseng KC, Liu YH, Kuo PL, Zheng HQ, Chang WC. Nucleic Acids Res. 2019 Jan 8;47(D1):D1155-D1163. doi: 10.1093/nar/gky1081. PMID: 30395277 ChIP-seq 実験のデータを元にして転写制御ネットワーク、 transcription factor binding sites (TFBSs) のデータベースを構築している。position weight matrix (PWM) が計算されている。
発現調節確率や相関係数は、遺伝子と遺伝子の間の結合を切るのに必要なエネルギーの大きさのようなものだと解釈することもできる。平均値を出す際に、エネルギーを重み付けとしてうまく取り入れる方法になっている。
「エンリッチメント」を指標とする分析法は、ごく普通の実験生物学の様々な分野に使えるような気がする。しかしそのためにはデータの分布を見るために、データを大量に取ることができないといけない。手間を掛けて精密な値を取るよりも、少々誤差が出ても大量に処理ができる実験法が望まれるのだろう。
遺伝子のマッピングはエンリッチメントそのものである。変異体と野生型（それぞれは異なる accession、例えば C型 と L型。二つの型は染色体全体に存在する多数のDNA マーカーで区別できるようになっている）を交配する。F2 世代の種子を収穫し、多数を播種育成して変異体の性質を示すものとそうでないものに分ける。変異体の性質を示すものでは、変異遺伝子の近傍にあるマーカーが C 型になっている割合が高くなる。変異遺伝子の部分は全部 C 型であるはずである。すなわちエンリッチされている。変異遺伝子と離れた部位では、C 型、ヘテロ接合、L 型が１：２：１で出てくる（偶然により期待される比率）。
mRNA のコドン使用頻度は「エンリッチメント」で分析できるかもしれない。生物時計に関わる重要な遺伝子では珍しいコドンが使われていて、それが様々な環境に適応するために有利に働いているという論文が紹介されていた。 不適合コドンと時計機能Nature 495, 7439 2013年3月7日 http://www.natureasia.com/ja-jp/nature/highlights/42351
Global transcriptomic profiling of aspen trees under elevated [CO(2)] to identify potential molecular mechanisms responsible for enhanced radial growth. Wei H, Gou J, Yordanov Y, Zhang H, Thakur R, Jones W, Burton A. J Plant Res. 2013 Mar;126(2):305-20. doi: 10.1007/s10265-012-0524-4. Epub 2012 Oct 13. PMID: 23065025
と言う論文では、「Protein domain enrichment analysis」という分析が行われていた。
A compendium of RNA-binding motifs for decoding gene regulation. Ray D ら Nature. 2013 Jul 11;499(7457):172-7. doi: 10.1038/nature12311. PMID: 23846655 と言う論文では、モチーフの保存度を表現するのに Z score が使われていた （Fig. 1)。
http://www.riken.jp/pr/press/2013/20130912_2/ 「キャッサバの系統間におけるDNA配列の違いを網羅的解析により同定」 理化学研究所 環境資源科学研究センターと同バイオマス工学のすばらしい成果
全遺伝子を14種類の機能グループに分類する。それぞれの割合が求まる。区画が14個あるヒストグラムのようにも見ることができる。次に「読み過ごし変異が起きている遺伝子（たぶん、こういう変異が起きる確率が高すぎず低すぎず、解析するのにちょうどよいのだろう）」を全遺伝子から取り出す。これも一つの「Gene set」と見ることができる。この Gene set についても、同様に機能グループに分類して割合を求める。「生物的・非生物的刺激応答（response to abiotic or biotic stimulus）」と「ストレス応答（response to stress）」の2つの機能グループには、「読み過ごし変異が起きている遺伝子」の割合が全遺伝子の割合よりも高く、濃縮されていた。こういう機能を持つ遺伝子は、変異・進化速度が速いということなのだろう。
In silico discovery of transcription regulatory elements in Plasmodium falciparum. Young JA, Johnson JR, Benner C, Yan SF, Chen K, Le Roch KG, Zhou Y, Winzeler EA. BMC Genomics. 2008 Feb 7;9:70. doi: 10.1186/1471-2164-9-70. PMID: 18257930 という論文では、Gene Enrichment Motif Searching (GEMS) という方法が用いられていた。
Plant-derived antifungal agent poacic acid targets β-1,3-glucan では、「Pathway Z score」というスコアが用いられていた。
Causal analysis （因果解析）approaches in Ingenuity Pathway Analysis では、「activation Z-score」とか何種類かの Z-score が用いられていた。
Context-specific transcriptional regulatory network inference from global gene expression maps using double two-way t-tests. Qi J, Michoel T. Bioinformatics. 2012 Sep 15;28(18):2325-32. doi: 10.1093/bioinformatics/bts434. PMID: 22962443 という論文では、Twixtrix という、転写因子と標的のネットワークを推定する方法が発表されていた。そこでは、ある遺伝子の発現が 2 つのセット間で最も大きな差が生じるように、複数のサンプル条件を二つのセットに分けることをまず行う。 その際に、 t スコアが用いられていた。
この論文の基本になる仮定として、「転写因子（制御因子）とその標的は、どちらも、その制御因子の遺伝子に特異的になるように作成したサンプルコントラスト（ = critical contrast, ある遺伝子の発現が 2 つのセット間で最も大きな差が生じるように、複数のサンプル条件を二つのセット C1, C2 に分けること） において、差次的に発現している（セット間の発現量に大きな差がある）」 ということを考える。 それ以外には仮定は必要ない。「Minimal Statistical model」と呼んでいる。仮定が少なければ少ないほど、余計なパラメーターを減らすことができる。 Overfitting などの問題も軽減される。と書かれている。
しかし試してみると、結局相関係数を計算するのとあまり変わらない結果になった。能書きが本当に正しいかはどうかと思う。しかし、せっかく試してみたので、そのことについて書いておく。 Twixtrixについて
この方法も、「発現量に偏りが生じるように設定した二つのデータセットを作成する」ということで、分布の偏りを利用していることになる。 そういう「アレイデータ分析に使える基本的な仮定」は、ほかにどんなものがあるだろうか。これは考える価値があるだろう。できるだけたくさん思いつかなければならない。 私はいまのところそういうことを文章で書いてみることしかできないが、様々な研究例を参考にして数学の形に書き直さないといけない。
しかし、マイクロアレイでの遺伝子発現の変化などは「処理区／対照」の対数として扱うことが多い。その場合は 0 を中心とした値になり、負の数も取る。
RNA-seq のカウント数はマイナスになることはないので、非負性を元にして分析ができる。"nmf rna-seq" で検索するとたくさん出てくる。
生物種もパラメーターと考えることができる。いくつかの生物種からのデータを統合することで信頼性の高い解析ができることが既に実証されている。 ATTED-II in 2016: A Plant Coexpression Database Towards Lineage-Specific Coexpression. Aoki Y, Okamura Y, Tadaka S, Kinoshita K, Obayashi T. Plant Cell Physiol. 2016 Jan;57(1):e5. doi: 10.1093/pcp/pcv165. Epub 2015 Nov 6. PMID: 26546318 パラメータを変化させるだけでなく、データにその性質を保つような何らかの変換・縮約を施して、それでも同様に検出されるかどうかという方法も考えられる。うまく縮約してから計算しても問題ないなら、計算が速くなって都合がよいだろう。 これもすでに調べられている。
Multi-dimensional correlations for gene coexpression and application to the large-scale data of Arabidopsis. Kinoshita K, Obayashi T. Bioinformatics. 2009 Oct 15;25(20):2677-84. doi: 10.1093/bioinformatics/btp442. Epub 2009 Jul 20. PMID: 19620096
SVD-based anatomy of gene expressions for correlation analysis in Arabidopsis thaliana. Fukushima A, Wada M, Kanaya S, Arita M. DNA Res. 2008 Dec;15(6):367-74. doi: 10.1093/dnares/dsn025. Epub 2008 Oct 17. PMID: 18931094
上の論文に出てくる SVD = 特異値分解について勉強したことを、忘れないようにここSVDに書いておく。
http://prime.psc.riken.jp/?action=drop_index 理化学研究所 PRIME では、様々なメタボロームデータが公開されている。普通体重や身長のような一般的な特徴量を多数測定すると平均値を中心として左右に対称に近い分布をしたヒストグラムが得られる。 しかしメタボロームのデータをみるとぜんぜんそうなっていない、指数分布または期待される発生回数が少ないポアソン分布のようにも見えるヒストグラムになっていることがある。 これはどういうことなのか。データの質が悪いのだろうか。 それとも対照の値や基準となる分子の量で割り算したりするとよいのだろうか。左右対称に近い、正規分布に近い分布をしたヒストグラムになっている方が回帰などの分析がやりやすい・正しい予測をしやすいというようなことが書かれた文章を読んだことがある。 ポアソン分布で分布が表されると言うことは、その物事が起きるかどうかの確率が低い・確率的なことが大きく関わっていると言うことだろう。そうなら決定論的な回帰式などが当てはまりにくくなる・予測が当たりにくくなるのだろう。だからそうならないように前処理を何か工夫しないといけないのだろう。
そのままの値ではなく、対数を取ると左右対称に近くなることがあることがわかった。対数正規分布に近いのかもしれない。 http://iss.ndl.go.jp/books/R000000024-I003953411-00 佐々木,小林,森山,松下先生による発表 「対数正規分布への再訪」で検索すると PDF ファイルがみつかる。 「過去の履歴に影響を受ける様々な現象が対数正規分布によって記述できる可能性がある」と書かれている。咀嚼によって食品を細かくする過程は、噛む動作一回ごとの履歴によって影響を受け、噛むごとに細かく破壊されていく。研究費獲得金額は、過去の研究成果、研究実績と獲得実績によってきわめて大きな影響を受ける。 メタボロームとは関係ないような気もするが、それが考え違いかもしれない。
と思ったが、ポアソン分布はパラメーター λ が大きければ正規分布とほとんど同じになるということを勉強した。単に測定の感度が足りなくて単位時間にその分子を検出するという事象がめったに起きていないというデータがあるのかもしれない。 ポアソン分布のように見えるデータはかなり以前に測定されたものだった。最近は技術が進んでいてそういうことは少ないらしい。
マイクロアレイのデータでは相関係数がとても有用に用いられ、すぐれた成果が多数公開されている。 Twixtrix も、ｔスコアを絶対値にせずに使うなら結局相関係数とほとんど同じであることがわかった（絶対値にしてもあまり意味がないような気がする）。 なぜマイクロアレイのデータでは相関係数がとても有用なのだろうか。
相関係数は負の値も取る。しかし遺伝子発現で負の相関があっても、それに意味があって生物学的な発見につながった例はあまり聞いたことがない。
負の相関が生じる原因は、いくつか考えられる。
このことは、物事を考える上であまり都合がよくない。同じ生物に関わることでも、バイオマスと関連がある細胞壁多糖類やデンプンなどの場合、ある成分が増加すると別の成分は減少すると言うことが起きることがある。 バイオマスというものはエネルギーの塊ととらえることができる。だからエネルギー保存則に従って、量にトレードオフが生じる。それによってある程度バイオマス量に関する定性的な予測が可能になる。 また、メタボロームのデータでは負の相関に意味があることが多い。それはバイオマスの場合と同じように前駆体が一定量存在してそれを取り合うような関係にあるときに生じる。 遺伝子発現のほうでそういうトレードオフが見つけられないだろうか。転写因子のような特別なカテゴリーの遺伝子に限れば、トレードオフが成り立っている例があってもよいのではないか。 川上博士らの wPGSA法のように、転写因子（制御因子）に特に注目するのはよい考えなのかもしれない。
遺伝子から転写された転写産物はさらに翻訳され、生じたタンパク質は化学反応を触媒したり構造を形成したりする。間接的になる。
GSEA 法では、遺伝子セットを決める必要がある。特にこうしなければならないと決まっているわけではない。Gene ontology データベースに頼るのは代表的な方法である。しかしそれでは満足できないこともある。公開されているマイクロアレイデータベースの結果などを用いて自分で決めた基準でセットを作ることもできる。それでうまくいくように思えるが、主観的、恣意的になりやすいということが問題になる可能性もある。Gene ontology 自体も、主観的に作られた物であることが問題点と言われているらしい。https://sites.google.com/site/scriptofbioinformatics/maikuroarei-guan-xi/go-jie-xi-r https://sites.google.com/site/scriptofbioinformatics/home 露崎弘毅博士によるサイト
https://kaken.nii.ac.jp/d/p/25117727/2013/3/ja.ja.html 野村 真樹先生の科学研究費報告書 遺伝子セット、遺伝子モジュールを得ることが、研究の一つの手段、目標として取り上げられている。
http://www.scls.riken.jp/scruise/software/eem.html 新井田 厚司先生によって、EEM というプログラムが開発されている。「遺伝子セット情報に基づいてmRNA発現データ中で共発現している遺伝子群、発現モジュールを抽出します」と書かれている。遺伝子セットを作る基準として「cisエレメントの解析」「ChIP-chip解析」などが考慮されていると書かれている。入力となる遺伝子セットは（Seed Gene Sets）と書かれている。まず種になるセットを作っておいてそれを解析したいデータの発現プロファイルに応じて変えていくように思える。
http://ci.nii.ac.jp/naid/110007226105 「関連遺伝子セットの多重解の存在 プリチャード真理 江口真透」 という文献があった。この場合「患者の状態を判別するために最適な遺伝子セットを作成する」ことが目標になっている。
転写産物の量を、「バネがついたおもりの、自然長からの伸び」と考えてみる。 バネのおもりと反対の側は、固定されている。バネの動きではフックの法則が成り立つ。 転写産物の種類ごとに、バネとおもりのセットになる。それがたくさん平行に並んでいる。 相関が高い遺伝子の集団は、お互いに近くにまとまっている。お互いに発現の仕方が似ていること、プロモーターの構造が似ていることに対応する。 それによってお互いに動きが同期するようになっている。 こういうことにすれば転写因子を考えなくても済む。
バネがついたおもりが振動している。バネの自然長からの伸びをランダムなタイミングで多数回観測して、それらの値をヒストグラムにしてみる。 転写産物の場合経時的に小さな時間間隔で何回も量を測定するのは難しいので、一度に多種類の転写産物を観測して、得られた値をヒストグラムにする（PAGE 法のやり方）ことで代用できないだろうか。 バネやそれに類するものは調和振動子と呼ばれている。細胞や組織も叩けば振動するので調和振動子として観測することもできる。組織の性質（硬さや水分含量）を観測する手法として使われている。
マイクロアレイのデータ一つは、細胞の一つの状態を表していると考えることができる。二つ（対照と何か刺激を与えたもの）を組み合わせると状態変化になる。
関数のように書いてみると、 細胞（ある状態変化）＝ 関数（遺伝子１の発現量変化，２の発現量変化，･･･、遺伝子 N の発現量変化） 最大で N = 数万の遺伝子が、パラメーターになることになる。以前の研究スタイルでは、数万の遺伝子から生物学的知識に基づいて恣意的にいくつかの遺伝子を選んでパラメーターとしていた（測定、議論の対象にしていた）。しかしそれではうまく解明できない難しい研究対象が見いだされてきた。そこで網羅的解析が導入された。
この場合 N = 数万を全部使うことが可能になる。バイオインフォマティクスの分野では、それによってきわめて有用な発見、優れたデータベースの構築がなされている。しかし、生物と関係ない分野の先生が書かれていたことだが、一般的にパラメーターは少ない方が物事を考えるにはよいらしい。また遺伝子発現量変化に含まれているノイズはそのまま入ってきてしまう。
「遺伝子セット」をうまく作れば、理屈上多数の種類の遺伝子からの情報を生かしながらパラメーターを減らすことができる。そのほうが分析しやすい場合があるかもしれない。また、遺伝子セットにすることで、遺伝子発現量に含まれているランダムなノイズをある程度打ち消すことができるかもしれない。 細胞（ある状態変化）＝ 関数（遺伝子セット１の Zscore，２の Zscore，･･･、遺伝子セットN の Zscore） N を数万よりは小さくできる。Zscore は標準正規分布をする数なので、扱いやすいかもしれない。 こういうことが役に立つかどうかわからないが、何かの足しになるかもしれない。この場合、遺伝子セットをどのように決めればよいかが問題になる。また、スコアの付け方も Zscore だけでなく、何らかの力学とか幾何学や熱力学のモデルに基づいたスコアの付け方を考えた方がよいかもしれない。そういうことができれば、他の分野のすぐれた成果にうまく乗っかることができるようになるかもしれない。
遺伝子セットと、それらから計算されるスコアをうまく作ることで、データから「温度変化」「体積変化」「圧力変化」「エントロピー変化」「エネルギー変化」に相当する状態量を取り出せるかもしれない？ 農水省の研究所では、野外で生育しているイネの生活環全体にわたって、マイクロアレイのデータを網羅的、組織的に採取することが行われている。同時に温度変化、日照時間、降水量の変化、花が咲く時期、収穫日、収量などのデータも得られる。それらによって得られた膨大なデータをうまく分析すると、イネの遺伝子発現と様々なデータの関係を見いだし予測することが可能になる。だからすでにそういう分析も行われている。 永野惇先生が、植物の生長調節2014年第2号（植物化学調節学会会誌）に、「野外でのトランスクリプトーム解析から見えること」というすばらしい解説を書かれている。イネのトランスクリプトームの野外環境応答の研究成果について解説されている。
Deciphering and prediction of transcriptome dynamics under fluctuating field conditions. Nagano AJ, Sato Y, Mihara M, Antonio BA, Motoyama R, Itoh H, Nagamura Y, Izawa T. Cell. 2012 Dec 7;151(6):1358-69. doi: 10.1016/j.cell.2012.10.048. PMID: 23217716
データベース FiT-DB: http:////fitdb.dna.affrc.go.jp/ によってデータも公開されている。
永野先生の解説では、トランスクリプトームデータと窒素栄養条件の相関を用いて、窒素栄養条件を良く表現する遺伝子セット指標を作成した研究が紹介されていた。
Gene expression biomarkers provide sensitive indicators of in planta nitrogen status in maize. Yang XS, Wu J, Ziegler TE, Yang X, Zayed A, Rajani MS, Zhou D, Basra AS, Schachtman DP, Peng M, Armstrong CL, Caldo RA, Morrell JA, Lacy M, Staub JM. Plant Physiol. 2011 Dec;157(4):1841-52. doi: 10.1104/pp.111.187898. Epub 2011 Oct 6. PMID: 21980173
遺伝子セットの組み合わせ（場合の数）はとても多いので、なにかの重要な指標となる数に相当する情報を含むものがその中には必ず存在することが理論的に示せるかもしれない？経済学では「市場均衡解の存在証明」とか「Nash 均衡の存在証明」とか「斉一成長経路の存在証明」とか「一般均衡価格体系の存在証明」とか、存在を証明することが多くてそのためのテクニックも発達しているらしい。これは数学の方法を使うらしい。「存在 証明 数学」で検索すると、有用な解説が出てくる。
http://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/39717/2/kubo2009IEICE.pdf 最近のベイズ理論の進展と応用 久保拓弥先生 「主観事前分布」を使ったベイズモデルというものが以前はよく使われたと紹介されている。主観的に決めた事前分布によって観測値を補正し、より正しい値を推定するような感じに思える。この方法の恣意的になりやすいという欠点を改善するために、階層ベイズモデルというものが考案された。それについて解説されている。遺伝子セットを決めることにそのまま使えるわけではないが、解析者の主観にできるだけ左右されにくい方法（データ自身に遺伝子セットを決めさせる）を考える参考になるかもしれない。
主観的に自分で決めた基準で遺伝子セットを作るのは、「主観事前分布」を使ったベイズ分析に相当するような気がする。そういうことを言っているらしい論文もすでにあった。
Metabolomics data exploration guided by prior knowledge. van den Berg RA, Rubingh CM, Westerhuis JA, van der Werf MJ, Smilde AK. Anal Chim Acta. 2009 Oct 5;651(2):173-81. doi: 10.1016/j.aca.2009.08.029. Epub 2009 Aug 25. PMID: 19782808
遺伝的アルゴリズムの研究が盛んに行われている。 http://mikilab.doshisha.ac.jp/dia/research/pdga/research.html 同志社大学 知的システムデザイン研究室 遺伝的アルゴリズム研究グループ こういうのもよりよい遺伝子セットを作るのに使えるのではないか。「GSEA にとって、よりよい遺伝子セットとはどんなものか？ 良さのスコアをどう付けるのか？」ということが問題になる。 「セットに含まれる遺伝子はできる限り性質が似ていることが望ましい」という基準もありうる。逆に「注目している性質以外は、できるだけ多様な遺伝子を含む方がよい」という基準もありうる。 熱力学では一種類の粒子に注目した場合、同じ種類の粒子はお互いに全く区別できないことを前提とする。それによって様々な状態量とそれらの間の法則が成り立つようになる。熱力学はすべての物事の基本である。だからこの場合も「セットに含まれる遺伝子はできる限り性質が似ていることが望ましい」という方針が考えられる。
ここにメンバーとして例えば30個の遺伝子を含むリストがあるとする。これらの遺伝子の「区別できなさ」を表す指標となる値を計算できるようにする必要がある。それはいろいろなものがありうるので、試してみないといけない。
「化学」2010年10月号に山本先生が「その実験データ、ちょっとどーなん！？」という記事を書かれている。エラーバー（ある信頼度の元で期待される誤差）について書かれている。t 分布の精神を元にするという方針に基づいている。書かれていることに少し補足をして、考えてみる。
遺伝子セットのような少数のメンバーを含むセットについてスコアを考える場合、
z スコア＝
t スコア＝
N はセットに含まれるデータの数
のうちで、t スコアの方を考える。
二つのデータセットを比較する場合、Reference value に「対照となるデータ群の平均」を当てはめる。それによって二つのデータセットの平均に差があるかないかを判断する基準にする。
今回の遺伝子セットの場合 Reference value は「真の値：その遺伝子セットに含まれるすべての遺伝子の発現量の平均値が本来取るべき値」で、その値と「セットの平均値」との差は「誤差」になる。
だから 誤差 ＝ t スコア * セットの SD / SQRT(N) この式が36ページの (1) として書かれている。
この場合 t スコアとして、t 表の値 t(自由度, (1 - 信頼度)) を使う。どうしてそれでよいのかよく理解できていない。
自由度は「平均値がわかっている」ということからその分 1 をサンプル数 N から引くので N - 1 になる。 エラーバーの大きさは信頼度 95% とすると t((N - 1), 0.05) * (セットの SD) / SQRT(N) になる。「真の値」がわからなくてもよいことになる。
遺伝子セットの場合、セットに含まれる遺伝子の数を増やしすぎるとセットの SD は大きくなる。 しかし遺伝子の数が少ないと t スコアの方が大きくなる。だから誤差を最も小さくする数があると推測できる。誤差が小さい＝真の値のあたりに固まっている＝区別できにくい というように考えてもよいように思える。しかしこういう考え方が最適かどうかはわからない。
「遺伝子セットを作る」ということは、「少数の遺伝子を選んでそれらに特別なラベルを付ける」ということである。これは判別分析のような方法と共通したところがある。そういう多変量解析の手法でうまく使えるものがあるかもしれない。一つの遺伝子がもつ性質・プロパティで正確に値が求められるものを、できるだけ多くそろえてみる（多変量のデータを用意する）のがよいのかもしれない。
「データを〜でソートして最上位50個」などの方法は、そのデータの特徴になる部分を取り出す方法の一つと考えられる。その場合、データ一セットのみから特徴になる部分を取り出そうとしている。「いくつかのデータをまとめて（複数セット）、それぞれの特徴がもっともはっきり区別できる／分散が最大になる」ような量を作り出すのは、多変量解析で主成分分析のような方法がある。そういう量を、限定された数の遺伝子、要素で作り出せるような方法を作るとよいかもしれない。
vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12:
