RMA(Robust Multichip Array)
R 言語 + BioConductorパッケージを使うと、自分の実験で得た、またはデータベース ArrayExpress で公開されている *.cel ファイルを簡単に RMA 法で処理できる。
バイオインフォマティクスの専門家である門田先生が書かれたかなり古いテキストを参考にしている。いまでは RMA 法よりも改良された方法が多数ある。それらの解説が門田先生のページにある。私は古い方法を使っている。
#BioConductorのインストール
#source("http://www.bioconductor.org/biocLite.R")
#biocLite() １回行えばOK
#処理したいcelファイルが置いてある場所を指定
setwd("D:/〜〜/〜〜〜/〜〜〜〜")
#The *.cel files were treated with R/BioConductor package as follows.
library(affy)
data <- ReadAffy()
data.out <- rma(data)
これだけで Background correcting, Normalizing, Calculating Expression が行われる。結果は対数で出る。私が行った実験の場合、これで出た結果は生物学者の視点から見て納得できるものになっていた。どういうことかというと、
ということがある。
ディファレンシャルハイブリダイゼーションはいまでは使われることがない方法である。当時のそういう方法では発現量が低い場合差が検出できないことが知られている（バックグラウンドのシグナルによる）。また普段の発現量はとても低く、何か刺激を受けた際に極めて強く発現誘導される遺伝子（ストレス応答性遺伝子に多い）が優先して見いだされてくる方法だった。
Google で "M-A plot" の画像を検索すると多数のプロットが見つかる。中には「こんなデータで単純に倍率を求めて遺伝子を選択したら大間違いの元だろう」と思えるような画像もある。倍率以外の方法で遺伝子を選択する手法がすでに開発されている。
結果は、タブ区切りのテキストファイルとしてセーブしてエクセルに読み込んで、注釈を付けたりすると便利である。これも公開されている資料の通りにしている。
write.table(data.frame(exprs(data.out) ,check.names=FALSE), "data_rma.txt", sep="\t", col.names=NA,quote=FALSE) #結果をdata_rma.txtに保存
rma 以外にも、expresso という関数もある。これを使うように説明されている資料も多い。しかしこの関数はパラメータが多数ある。それらを変えることによって結果がそれぞれ少しずつ変化する。どれが正しいと言うことはできない。生物学的な解釈から判断するしかない。私がパラメータを変えて比べてみた限りでは、生物学者にとって面倒なだけで特にメリットはない。RMA関数を何もパラメータをつけずにそのまま使って問題を感じることはなかった。
RMA法についてよく知らないのでは間違いの元なので、しくみについてもさらに知らなければならない。
http://topepo.github.io/caret/other.html というページに記事がある。
Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Biostatistics. 2003 Apr;4(2):249-64. PMID: 1292552
GeneChip のデータではPM（パーフェクトマッチ） 値とMM （ミスマッチ）値が得られる。RMA 法は PM 値のみを使う。 Irizarry らの論文を見ると、(PM-MM) の値を使うとかえってデータの質が悪くなるというデータが示されているように見える。MM はミスマッチだから PM よりもずっとシグナルが下がるはずだが、そうなっていないことが多い。
単に PM の値を使うと、PM の値が大きくなるにつれて標準偏差 SD がとても大きくなる。それは望ましいことではない。PM の値を対数に変換すると、PM の値が大きくなっても SD の値の上昇は 1.5 倍に収まると書かれている (p254 の上)。これは、「ノイズはシグナルの〜％（かけ算）」という仕組みで観測される値が決まるかららしい。
どんな測定でも観測される値は、真の値に系統的な誤差、偶然による誤差が影響した値になっている。影響の仕方は足し算で影響することもあればかけ算で影響することもある。マイクロアレイの結果の場合、Scanner effect というものがあると書かれている（254ページ）。これは系統的な誤差をもたらす。
RMA 法では、background adjustment、Quantile Normalization、summarization という３つの段階がある。http://www.biomatrix.co.jp/product/download/data/FAQ_GeneSpring.pdf （株）バイオマトリックス研究所から、わかりやすい説明が公開されている。８ページからRMA 法の説明がある。これを読めばよい。 と思ったが、やはり Irizarry らの論文をよく読まなければいけなかった。しかし単に RMA 法を値を入れれば結果が得られるブラックボックスのように使うなら、そこまですることはない。
background adjustment について自分で書いてみる。
１つのプローブの、観測されたシグナル値（バックグラウンド値を含む）＝ シグナル値 ＋ バックグラウンド値 バックグラウンド値は、プローブが基盤に非特異的に吸着するとかの原因で生じる。
i = 1..I （アレイの枚数）, j = 1..J （遺伝子ｎに対してプローブがＪ個存在する）, n = 1..N （遺伝子１からN）
ここで、「シグナル値の内で、ノイズの影響が少ない、最大値に近い部分は指数分布する」「シグナル値の内で、最小値に近い部分の値＝バックグラウンド値は正規分布する」と、想定する。これは実験的に得られている結果から導かれたものである。シグナル値全体の分布ではなく、最大値へ向かう方向の分布と最小値へ向かう方向への分布をうまく使い分けている。他のデータの分析にも応用できるかもしれない。
Irizarry らの論文では、Fig. 5 にそのヒストグラムがある。左端に近いところに山の頂上があり、そこから左端へは急に落ちて、右端に向けては指数関数的になだらかに下がっている。こういうヒストグラムは対数正規分布と呼ばれるものと似ている。しかしそれとは異なってシグナル値の部分はもっと頻度が高いことも、Fig. 5 には示されている。このことから単一の分布ではなく、二つの分布が足しあわされたものと想定することができる。文章では、p260 の中央あたりに書いてある。
（勉強中）
Quantile Normalization について自分で書いてみる。
Quantile は「分位」と訳される。データをいくつかにできるだけ均等に分割することを考える。２つに均等に分割するポイント、数を二分位数という。中央値 median がそれに相当する。４つに均等に分割するなら四分位数が第一、第二、第三と生じる。
しかし、この場合は quantile を「順位」と考えた方がよいかもしれない。
それぞれのアレイ（一枚ずつ）について、シグナル値のヒストグラムを書いてみる。補正していない状態では最小値、最大値がアレイごとに異なりずれている。最小値（発現量ランキングの順位は最下位）から最大値（順位は第一位）まで、順位が同じ遺伝子では発現量の値が一致するように補正する（ように思える。最終的に得られる値では、そこからまた補正されるので順位が同じでも発現量の値は少しずつ異なってくる）。
補正はどう行われるか。複数のアレイデータで、同じ順位の遺伝子の発現量をそれぞれのアレイから取り出し、それらの中央値に値をそろえる？ （というように想像する。勉強中）。
Irizarry らの論文では、p254 の下の方に quantile normalization について書かれている。
（勉強中）
結果的に、すべてのアレイデータで、発現量ランキングの順位が同じならシグナル値は同じにされる（ように思える。最終的に得られる値では、そこからまた補正されるので順位が同じでも発現量の値は少しずつ異なってくる）。シグナル値の分布は同じ形になる。 その分布を、「横軸は発現量ランキングの順位」「縦軸はシグナル値」にしたグラフに書いてみると、それによってできるグラフの曲線はシグモイド関数に似た形を取る（そう考えるとつじつまが合う。また実際に書いてみると確かに最上位、最下位では傾きが小さくなっている。特に最上位のほうは著しい）。 シグモイド関数の傾きは、ランキングの順位が全体の中央で最大になり、最上位、最下位のあたりでは傾きが小さくなる。
どうしてそうなるか。どんな測定でも得られる測定値、シグナルは測定対象物の量と完全に比例するわけではない。あるレベルを超すと、シグナル値が飽和してそれ以上上がらなくなることがほとんどである。 マイクロアレイの場合も発現量が高い部分では、シグナル値が飽和する。そのため発現量ランキングの順位トップのあたりでは、順位がかなり変動してもシグナル値の変化は小さい。 発現量が一番少ないランキングの順位最下位のあたりでは、ノイズレベルよりも値が下がることはない。逆向きに飽和していることになる。そのため最下位のあたりでも、順位がかなり変動してもシグナル値の変化は小さい。
RMA 法で処理したデータを M-A plot すると、発現強度の順位が中央の領域で点のばらつきが最大になり、両端（最大値、最小値）では０に収束したような形になる。 これは、上に書いたことによる効果による。シグモイド関数の傾きをグラフにすると、中央の部分で値が最大になる。両側は０に収束する。
これは倍率を指標として分析するときに都合がよい性質である。
実際に発現量ランキングで最下位に近い（ほとんど転写されない）遺伝子は、どんな条件でも転写されないままで最下位に近いことが多いのかもしれない。発現量ランキングで最高位に近い遺伝子は、転写のレベルでも飽和していて変化しにくいということもあるのかもしれない。
RMA 法で処理したデータは、複数のアレイデータのシグナル値がお互いに近くなるように補正される（ヒストグラムの形が同じになることで、自然にそうなるだろう）。それにも関わらずシグナル値がずれている遺伝子を選ぶことで、本当にシグナル値に違いがあるものを選べる可能性が高くなっているのかもしれない。 「分布がシグモイド関数の形を取る」というような情報から、その分布が形成される機構に関する情報を得られる可能性がある。天文学では星や星間分子を観測して位置、密度の分布を調べる。そこから理論、法則につなげることができるらしい。
Science 11 April 2014: Vol. 344 no. 6180 pp. 183-185 DOI: 10.1126/science.1248724 Unfolding the Laws of Star Formation: The Density Distribution of Molecular Clouds Jouni Kainulainen, Christoph Federrath, Thomas Henning
http://www.sciencemag.org/content/344/6180/183.abstract
シグモイド関数は、微生物の放射線に対する生存曲線でも出てくる。それは「マルチヒットモデル（ある程度強い放射能で細胞内の複数の標的が同時に損傷することで、効果が出る）」とつなげられる。標的が一つでよいなら指数関数になる。 http://www.agr.nagoya-u.ac.jp/~food/Dr.Namiki%20Review.pdf 並木満夫先生 プロモーターから転写が開始される際も、プロモーターに存在する複数のエレメントに同時にそれぞれに対応する転写因子や RNA ポリメラーゼが結合する必要がある。そう考えると、転写も「マルチヒットモデル」があてはまる。そのためにシグモイド関数が出てくるという考えもあり得る。
シグモイド関数の形は、正規分布などの分布の累積分布関数として出てくる。
群馬大学の青木先生の解説で、プロビットモデル Probit model というものが解説されている。「許容量」は正規分布で、用量 - 反応曲線はシグモイド関数の形になっている。
（勉強中）
Summarization について自分で書いてみる。
１つのプローブのシグナル値＝ RMA値 ＋ プローブによる偏り ＋ ノイズ
i = 1..I （アレイの枚数）, j = 1..J （遺伝子ｎに対してプローブがＪ個存在する）, n = 1..N （遺伝子１からN）
上のように「シグナル値がどのようにして決定されるか」をモデル化している。ここで出てくる RMA 値等は対数なので、足し算は元の発現量にとっては乗算に相当する。logAB = logA + logB このことは、「発現量が大きい遺伝子は、プローブのアフィニティーの違いによるシグナルの変化、ノイズによるシグナルの変化も大きくなる。」と書き直すことができる。「ノイズはシグナルの〜％」というように割合で考えるというモデルになっている。
たくさんのパラメーターが生じるがそれらの値を決めるために、いくつかのこと（拘束、制約する条件）を仮定する。
ノイズは平均値０でシグナルと独立に分布する（どんな分野でもよく使われている仮定：プローブ一つずつに対してノイズ値を割り当てる。それらを全部平均すると０になり、ノイズ値の分布はある分散をもつ正規分布になる）。
とする（遺伝子一つについて複数（Ｊ個）のプローブがあるが、Ｊ個のプローブでアフィニティ、偏りを平均すると０になる。これらは複数のアレイで同じ値であるとする）。
RMA値は0以上の値だろう。RMA値は１枚のアレイ上の遺伝子一つについて値が1つだが、シグナルの方は１枚のアレイ上の遺伝子一つについて値がＪ個存在する。そのため、「ノイズは全部0」「アフィニティ値も全部0」では左右の値が等しくなることはない。条件をさらによく満たすようなパラメーター値が存在し、その値が決められる。
その方法に、median polish と呼ばれる手法が使われる。"median polish microarray" で検索すると、Methods in Microarray Normalization という本の一部分が出てくる。"メディアンポリッシュ" で検索すると、「ＲとＢｉｏｃｏｎｄｕｃｔｏｒを用いたバイオインフォマティクス」という本の一部が出てくる。
行列として並んでいる数値に median polish を行うと、（行、列）ごとに一つづつ割り当てられた数＋それぞれの行の値＋それぞれの列の値＋ 切片に相当する一つの数 という具合に、それぞれの要素を４つの因子に分解して表現できる。数値を行列として並べたものを、この方法で分解する。分解して生じた因子を、アフィニティー（偏り）に対応させる。
メディアンポリッシュは、行列になっているデータに対して、「行の効果」「列の効果」を取り出すために使われる。何かレースがあったとする。行がスタート位置、列が第１レースから第１２レースの結果とする。もし特定のスタート位置から出た選手が勝つ頻度が高い（偏りがある）なら、それは「行の効果」として検出される。
RMA 法では、プローブを行に、繰り返し実験を列に配置して行列を作る。これにメディアンポリッシュを適用すると、「プローブの効果（プローブによる偏り）」「実験の効果（それぞれの実験固有の偏り）」が検出される。
「プローブの効果（プローブによる偏り）」の値を除外し、残りの値からそれぞれのプローブに対応するシグナル値を計算する。一つの遺伝子について複数のプローブがあるので、それらから得られるシグナル値を平均することでランダムな誤差を打ち消すことができる。こうして RMA 値が計算される。
メディアンポリッシュという操作については、もう少し考えてみる必要がある。
（勉強中）
RMA 値は、生の数値からプローブの効果による系統的な誤差（偏り）、ノイズによるランダムな誤差（平均すると０になる）を取り除いた値になっている。それによって信頼性が高い、もっともらしい結果が得られる。
RMA 法はうまく働いて生物学的な知見をデータから取り出すのに有効であることが実証されている。しかし私は値を入れたらそれに対応する値が出てくるブラックボックスのようにしか使えていない。この処理で得られる値と細胞内の仕組みをつなげること、また化学反応の基盤にある熱力学とつなげることができるのか、いまのところわからない。
（勉強中）
RMA 法の問題点： この方法で出てくる値は生物学的に見るともっともらしいものになる。しかし複雑な処理を何段にも行っているのでそれらは単純な数式だけで表現しにくい。何かもっと単純な方法で RMA 法と相関が高い結果を出してくれる回帰式のようなものが出来たら、その係数から出てきた値の解釈がしやすくなるかもしれない。
適当に事前知識に基づいて様々なパラメーターを恣意的に考える。野球のデータ解析では出塁率とか様々なものが次々と考案されている。それらを用い、RMA 法によって計算された値と相関が高い結果を出して、パラメーター数ができるだけ少ない回帰式を見つける。 その回帰式に選び出されたパラメーター（記述子）が、それぞれどんな意味を持っているかを考え解釈する。これは人間が行う必要がある。
vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12:
