QTL 解析の統計モデルと検定の多重性調整
栗木 哲∗

概 要
個体の形質を規定する遺伝子 (量的形質遺伝子座，QTL) を探索するための統計的手法を QTL 解析という．QTL 解析においては，ロッドスコアとよばれる尤度関数の極大点を探索することによって，QTL の位置が推測される．本稿では，実験交配生物
に対する QTL 解析について，その統計モデルと，QTL 探索においてしばしば問題となる検定の多重性調整について解説する．1 章では，QTL 解析が，多重比較・多重検定の一種であることを説明し，問題の枠組みを与える．2 章では QTL 解析のための代表的な統計モデルを，背景となる遺伝知識と併せて可能な限り簡潔に説明する．また各モデルに対して定まるロッドスコアの確率構造を導く．3 章では 2 つの数学的手法 (非線形再生理論およびオイラー標数法) によって，QTL の有無判定のためのロッドスコ
アの閾値が合理的に設定されることを見る．

1

はじめに

1.1

QTL 解析と変化点問題

個体のある形質 (形態や性質) が遺伝的な効果によってもたらされると考えられる場合，その原因となる遺伝子を探し出すことは，遺伝学研究の重要な目的の一つとなる．その形質が主として連続量で記述され，また一般には複数の遺伝子と環境要因によって規定されるものである場合，量的形質とよばれる．たとえばマウスの脂肪体重比 (肥満度) は，典型的な量的形質である．量的形質の原因となる遺伝子が，QTL (量的形質遺伝子座，quantitative

trait loci) である．
QTL を探索するための統計手法を QTL 解析という．QTL 解析は連鎖とよばれる遺伝現象を積極的に利用した，代表的な連鎖解析である．QTL 解析においては，連鎖と形質発現の双方を確率的な現象と捉えて統計モデルを設定することによって，目的遺伝子の探索が行われる (鵜飼 (2000) , Wu, et al. (2007))．
最初にマウスの肥満の原因となる遺伝子を探索するためのデータ解析の例を示す．ここでは肥満度の代用特性である血中アディポネクチン濃度 (単位 log10 [ng/m]) に着目し，それを量的形質としている．また解析対象のマウスは，標準的マウス近交系である B6 と，日本産亜種由来の MSM 系統の雑種 206 個体である．この 2 系統は形質が多くの点で対照的であるため，QTL 解析に適したものである．解析結果は，図 1 のロッドスコア (LODscore) に要約されている．図の横軸はマウスの 20 対の染色体における遺伝子座の位置であり，その点に位置する遺伝子の遺伝子型と，血中アディポネクチン濃度とのある種の連∗ 統計数理研究所，総合研究大学院大学，〒 106-8569 東京都港区南麻布 4–6–7，Email:1

kuriki@ism.ac.jp

関の尺度 (ロッド) がプロットされている．この図によると，第 3 および第 16 染色体上に

QTL が存在することが示唆される．

5

lod.Log_ADI

4

3

2

1

0
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

X

図 1: ロッドスコア
ところでこの解析例のように，観測値として与えられている系列や関数のデータがある
時点において変化を示すと考えられるときに，その変化時点を統計的に推測する問題を変化点問題という．一般に変化点問題では (i) 変化点の有無の判断のための統計量の閾値の設定，(ii) 変化点の位置の区間推定，(iii) 複数の変化点が存在する可能性があるときはその個数の推測，などが問題となる．本稿では QTL 解析および関連する連鎖解析において，比較的研究が進んでいる (i) の閾値の設定という問題に焦点を絞り，現時点で知られていることやその背景となる数理について概観する．残念ながら (ii), (iii) については現時点では不十分な結果しか知られておらず，本稿ではほとんど触れることはしない．例えば汎用的手法とされるブートストラップもこれらの解決の役には立たない (Manichaikul, et al.

(2006))．変化点問題は，より広い立場では特異モデルとよばれるクラスの統計モデルであるが，特異モデルにおいては，正則なモデルで成り立つ種々の漸近的性質が成り立たない(福水，他 (2004))．変化点問題において，ブートストラップやモデル選択規準を，少なくとも正則なモデルと同じ形で用いることはできないのは，このことによる．
ところで，連鎖解析，QTL 解析は，ヒトを対象にするものと実験交配が可能な生物を対象にするものに大別され，その両者では解析の対象とするデータの形が大きく異なる．本稿では後者に対するものを想定する．

1.2

多重性調整 (有意水準の調整)

図 1 のようにロッドスコアが明確なピークを持つ場合，その位置付近に QTL が存在すると推測される．しかしながら連鎖や量的形質の発現は確率的な事象であり，それゆえロッ
ドスコアもランダムなグラフである．現われたピークがランダムなゆらぎによる見せかけ2

のものでないかどうかを判定するためには，そのための基準，すなわち閾値を合理的に決める必要がある．
ロッドスコアの定義域を Γ とし，ロッドスコアを LOD(γ), γ ∈ Γ で表わすとする．LOD(γ)
は，“QTL が γ 付近に存在しない” という帰無仮説に対する検定統計量であり，その値がある閾値 c を越えたときに位置 γ の付近に QTL が存在すると判定することができる．ロッ
ドスコアのピークを探索し，その付近に QTL が存在するかどうかを判定することは，全ての γ ∈ Γ について仮説を同時に検定していると考えられる．そのことから多重検定の考え方によって，閾値 c を決めることができる．帰無仮説を

H0 : QTL がどの位置にも存在しないとし，その仮説が真であるにも関わらずロッドが閾値 c を越え QTL がどこかで発見される (正確には，発見されたと判断される) 事象を偽陽性 (false positive) と定義する．偽陽性確率を α (0.05 あるいは 0.01 など) 以下に調整するためには，閾値 c を


(1.1)
P ∃γ, LOD(γ) ≥ cα | H0 = αとなる c = cα とすればよい．この偽陽性確率は，弱い意味での FWER (family-wise error
rate in the weak sense) ともよばれる (Hochberg & Tamhane (1987))．上式は


(1.2)
P max LOD(γ) ≥ cα | H0 = αγ∈Γと書き換えられる．すなわち閾値 cα は，確率過程 LOD(·) の最大値の上側 α 点である．
ところで γ を固定すれば，統計量 LOD(γ) に対する有意点として

P LOD(γ) ≥ c̃α (γ) | H0 = α
をみたす点 c̃α (γ) が定義される．LOD(·) が確率 1 で一定値をとることがない限り，

maxγ∈Γ LOD(γ) は LOD(γ) よりも確率的に大きな値をとり，cα > c̃α (γ) (0 < α < 1) である．逆に，多重検定であることを考慮しないでロッドスコアの棄却点として c̃α (γ) を用いると，偽陽性確率は α をこえてしまう．この現象を検定の多重性という．また水準 α の有意点として，c̃α (γ) の代わりにより値の大きな閾値である cα を用いることを，多重性調
整，あるいは有意水準の調整という．
多重検定の多重性調整の方法としてボンフェロニ法が良く知られている．これは，検定の回数 (Γ の要素数) を |Γ| とおくとき，全ての γ ∈ Γ について LOD(γ) を水準 α/|Γ| で検定する，すなわち棄却点として c̃α/|Γ| (γ) を用いる方法である．このとき

(1.3)



≤
P ∃γ, LOD(γ) ≥ c̃α/|Γ| (γ) | H0 P LOD(γ) ≥ c̃α/|Γ| (γ) | H0
γ∈Γ

=



α/|Γ| = αγ∈Γであるので，偽陽性確率は α 以下に調整される．しかし後で詳しく見るように，LOD(·)
は強い相関を持った確率過程であるため，(1.3) の左辺の偽陽性確率は α よりも非常に小さな値となり，それにともなって QTL の検出確率 (検出力) も小さなものとなる．とくに

3

マーカー数が多いとき，あるいはエピスタシスとよばれる複数の QTL による交互作用を検定するときには検定の回数が莫大となりこの傾向が顕著となる．また，2.6 節で説明する
区間マッピング法では，ロッドスコアはマーカー間で連続的に補間されるため，|Γ| = ∞,c̃α/|Γ| (γ) = ∞ となり，ボンフェロニ法は意味をなさなくなる．以上の理由から，QTL 解析においてはボンフェロニ法を用いることはできない．
次の 2 章では，実験交配における QTL 解析と関連する連鎖解析の統計モデルをいくつか紹介する．さらにそれらのモデルにおいて現われるロッドスコア LOD(·) の確率過程としての構造を調べる．

3 章では，2 章で与えたロッドスコア LOD(·) の構造から，その最大値 maxγ∈Γ LOD(γ)の上側 α 点 cα を求める方法を説明する．経験則やシミュレーションに基づく方法に触れた
後，理論的な近似法について解説する．確率過程，確率場の最大値の分布については長い研究の歴史がある一方で，近年においても本質的な進展が見られている (Siegmund (1985),

Piterbarg (1996), Adler & Taylor (2007), 栗木・竹村 (2007))．本稿では，逐次解析，非線形再生理論を用いる方法と，オイラー標数法とよばれる積分幾何的な手法の 2 通りによって，確率過程としてのロッドスコアの最大値分布の近似を与え，そのことを通して多重性調整が可能であることを見る．

2

QTL 解析の統計モデルとロッドスコア

2.1

データの形

ここでは QTL 解析が対象とするデータの形と，その背後に想定される確率構造を説明する．交配の実験計画として，BC (戻し交配, backcross) と F2 の 2 種類を考える．(これらの実験交配については，次節で説明する．)
個体数を n とする．またマーカー遺伝子座の数を m とする．マーカー遺伝子座 (しばしばマーカーと略す) とは，何らかの方法でその遺伝子型が観測できる遺伝子座をいう．個体のそれぞれ t = 1, . . . , n について，着目している量的形質 (表現型) の測定値 y (t) (スカ
(t)

(t)

ラ) と m 個のマーカー遺伝子に対する遺伝子型のベクトル z (t) = (z1 , . . . , zm ) が得られ(t)ている．遺伝子型 zi

は戻し交配の場合は 2 値，F2 の場合は 3 値をとる．取り扱いの容

易さから，
(t)

zi

= 1, −1 (戻し交配の場合),

1, 0, −1 (F2 の場合)と表わすことにする (表 1)．
表 1: 個体データ
個体番号

表現型

1
..
.
n

(1)

y

..
.
y (n)

遺伝子型

z

(1)

(1)

(1)

= (z1 , . . . , zm )
..
.
(n)
(n)
z (n) = (z1 , . . . , zm )4

これらの個体データとは別に，m 個のマーカー遺伝子 i = 1, . . . , m のそれぞれについて，それが属する染色体の番号 ci と，その染色体上での位置 di の情報が与えられている．

di は基準となる点からの遺伝的距離 (単位はモルガン M，またはセンチモルガン cM，これらの意味は次節で説明する) で記述され，同一染色体の中では昇順 (i < j ならば di < dj )とする (表 2)．
表 2: マーカーデータ
マーカー番号

1

2

···

m

マーカー名

∗

∗

···

∗

染色体番号

1

1

···

c

座の位置 (M)

d1

d2

···

dm

QTL 解析では，個体データである遺伝子型と表現型の組 (z (t) , y (t) ) を確率変数と考え，モデル化を行う．ただし通常の多変量解析と同様に，個体間では独立と考える．
以下では，m 次元の遺伝子型ベクトルデータ z (t) の周辺分布について説明する．これは連鎖によって引き起こされるものである．遺伝子型 z (t) が与えられたときの表現型 y (t) の分布を表現するための統計モデルについては後の節で説明する．マーカーの位置 di の単位はモルガンとする．
マーカー遺伝子の添字 i = 1, . . . , m に対応させる形で，±1 に値をとる確率変数 i (i =1, . . . , m) を考える．ただしこの列はマルコフ系列で，⎧ ⎨ 1 1 ± e−2(di+1 −di )1
P (1 = ±1) = , P (i+1 = ±i | i ) = 2⎩1
2(i, i + 1 は同じ染色体上),(i, i + 1 は異なる染色体上)2で定義されるものとする．同時確率分布は


1−2(di+1 −di )
1
+


e
i i+1
2m i=1
m−1(2.1)

P (1 , . . . , m ) =(ただし座 i と座 i + 1 が同じ染色体上にないならば di+1 − di = ∞ とおく) である．任意の i, j についてP (i = ±1) =1,
2

P (j = ±i | i ) =⎧ ⎨ 1 1 ± e−2|dj −di |2(i, j は同じ染色体上),⎩12(i, j は異なる染色体上)となることに注意する．次に (δ1 , . . . , δm ) ∈ {−1, 1}m を (1 , . . . , m ) と独立に同じ分布に
(t)

(t)

従うランダムベクトルとする．このとき，遺伝子型 z (t) = (z1 , . . . , zm ) に仮定される確率モデルで最も基本的なものは，

(2.2)

(t)

d

(t)
(z1 , . . . , zm
) =
d

=

(1 , . . . , m )(戻し交配の場合),
1(1 + δ1 , . . . , m + δm ) (F2 の場合),
2(各 t について独立に)5dと表わされる．ここで = は両辺の分布が等しいことを意味する．
遺伝子型のこのような確率構造は，連鎖により引き起こされるものである．次節ではそのことについて説明する．

2.2

実験交配と連鎖

一対の染色体 (相同染色体) の一本は母親由来，一本は父親由来である．各個体の遺伝
子型を，記法 A1 B1 · · · /A2 B2 · · · によって表わす．これは，一方の親に由来する染色体の遺伝子型が A1 B1 · · · ，もう一方の親に由来する染色体の遺伝子型が A2 B2 · · · であることを意味するものとする．ここで A1 と A2 のペア，あるいは B1 と B2 のペアは同じ座に位置する一対の遺伝子 (対立遺伝子) である．全ての遺伝子座についてその遺伝子型がホモ(すなわち A1 = A2 , B1 = B2 , . . .) であるとき，近交系という．
近交系は植物の場合は自殖，動物の場合は兄妹 (けいまい) 交配を繰り返すことによって容易に作り出すことができる．実際，自殖あるいは兄妹交配を，遺伝子型を状態とする
マルコフチェインによってモデル化すると，その吸収状態が近交系に対応することが確認できる．(遷移行列の固有値の絶対値で 1 でないものの最大は，自殖あるいは兄妹交配のそれぞれについて 0.5, 0.809 であり，前者の方が収束速度が速い．)
ある近交系の個体 P1 と別の近交系の個体 P2 を掛け合わせた雑種第 1 代を F1 世代という．F1 個体とその親 P1 (P2 ) との掛け合わせを戻し交配 BC1 (BC2 )，また F1 個体同士の自殖，あるいは兄妹交配によって得られる雑種第 2 代を F2 世代という．
以下では 2 つの遺伝子座に着目する．近交系においては，どの遺伝子座においても遺伝
子型はホモであるので，異なる近交系の個体の遺伝子型は

P1 = AB/AB,

P2 = ab/abと書くことができる．また P1 と P2 の配偶子 (生殖細胞，すなわち卵，精子) の遺伝子型
はそれぞれ AB, ab であるので，雑種第 1 代の遺伝子型は

F1 = AB/abとなる．

F1 個体の配偶子の遺伝子型としては，次の 4 種類
(2.3)

F1 の配偶子 = AB, ab (それぞれ確率 1−r
2 で),Ab, aB (それぞれ確率 2r で)が現われる．その理由を述べるために，図 2 に沿って配偶子の生成過程 (減数分裂) を説明しよう．
減数分裂では最初に相同染色体の対が分離し，4 本の染色分体となる．さらに互いに別の親に由来する 2 本の染色分体 (図中で色の異なるもの) は，図のように交差を起こす場
合がある．(これは確率的な事象である．) その後 4 本の染色分体は互いに分かれて F1 個体の配偶子が生成される．いま着目している 2 座の間で，奇数回の交差が起きたとすると，生成される生殖細胞の遺伝子は Ab または aB となる．この現象を組換えという．この結果として，F1 個体の配偶子の遺伝子型は 4 種類 AB, ab, Ab, aB でその頻度は (2.3) の通

6

りとなる．なお r は 2 座間で組換えが起きる確率であり，組換え価とよばれる．P1 , P2 の各個体においても同じ過程で配偶子が生成されるが，4 つの染色分体の遺伝子型は同じであるので結果として組換えは観察されないことに注意する．

$

%

$

%

D

D

E

$

%
E

D

E

$

%

ίᄩᴾྙẅ(1-r ) /2ὸ

D%ίᄩᴾྙẅᴾᴾᴾr/2

ὸ

$

E

ίᄩᴾྙẅᴾᴾᴾr/2

ὸ

D

E

ίᄩᴾྙẅ(1-r ) /2ὸ

図 2: 減数分裂と交差
さらに P1 = AB/AB と F1 = AB/ab の戻し交配 BC1 を考えると，その遺伝子型は，

BC1 = AB/AB, AB/ab (それぞれ確率 1−r
2 で),AB/Ab, AB/aB (それぞれ確率 2r で)の 4 通りとなる．1 つの遺伝子座に着目した場合，遺伝子型は A/A または A/a の 2 通り
であり，一般性を失うことなく 1, −1 で表わすことができる．
また F1 個体同士の掛け合わせにより得られる F2 個体では，4 × 4 = 16 通りの遺伝子型が現われる．ただし通常の方法では，2 種類のヘテロ A/a と a/A は識別されない．A/A,

A/a, a/a の 3 通りについては，そのマーカー遺伝子が共優性であれば識別することができる．それらを一般性を失うことなく 1, 0, −1 で表わすことにする．2 つの座に着目した場合は，識別できる遺伝子型は 3 × 3 = 9 通りとなる．
交差の確率モデルとして最も基本的なものは，交差の生起をポアソン事象と考えるものである (Haldane (1919))．交差をポアソン事象と考えることにより，染色体上の 2 点間でおこる交差の平均値をもって，その 2 点間の距離と定義することができる．この距離が遺伝的距離で，その単位はモルガン (M) である．この遺伝的距離を遺伝子座の距離と定義
することによって，交差の生起事象は強度関数が 1 の定常ポアソン点過程と考えることができる．x (M) 離れた 2 点で，i 回交差がおきる確率は xi e−x /i! であるので，2 点間の組
換え価 r = r(x) は(2.4)

r(x) = P (奇数回交差が起きる) = xi
i: odd

i!

e−x =

1
(1 − e−2x )2となる．関数 r(x) をホールデンの地図関数という．

F1 個体の配偶子の遺伝子型の確率構造を考える．マーカー遺伝子座 i の遺伝子型が P1由来であるとき i = 1，P2 由来であるとき i = −1 とする．ポアソン性の仮定の下で，(1 , . . . , i ) と i+1 − i は独立であるので 1 , . . . , m はマルコフ性を持ち，またP (i+1 = i | i ) = P (座 i, i + 1 間で組換えは起らない) = 1 − r(|di+1 − di |),すなわち 1 , . . . , m は確率分布 (2.1) に従う ±1 列である．

7

遺伝子型 zi の定義の仕方より，戻し交配 BC においては F1 個体の配偶子の i = 1, −1の値が，BC の遺伝子型に一致する．また F2 個体においては， 12 (i + δi ) = 1, 0, −1 の値
が遺伝子型となる．以上で (2.2) が導出された．
組換え価 r(x) は，0 ≤ r(x) ≤ 1/2 に値をとる単調増加関数である．2 座が同じ遺伝子
座 (x = 0) のとき組換えは起らず (r = 0)，また 2 つの遺伝子座が非常に離れている，あるいは別の染色体上にあるとき (x = ∞) は確率 r = 1/2 で組換えが起きる．2 つの遺伝子座が近くにあり，組換え価が小さな値をとっている状態を，2 座が連鎖するという．染色体の長さは典型的には 100cM (=1M) 程度であり，交差は平均 1 回しか起きない．そのため実験交配の個体の遺伝子型は強い正の相関を持ち，多重性の調整においてはそのことの配慮が必要となる．
注 2.1 ここでは交差をポアソン点過程によりモデル化しているが，より一般には再生過程を用いてモデル化することができる (Karlin & Liberman (1983))．

2.3

単一マーカー分析

本節では QTL 解析の統計モデルとして基本となる単一マーカー分析 (single marker

analysis) を説明する．とくに断らない限りは F2 集団を扱う．
QTL 解析が対象とするデータは，表 1，表 2 の形であった．単一マーカー分析とは，マーカー遺伝子座 i の遺伝子型で集団を 3 群に分け，その 3 群で形質 (表現型) の平均が等しいという仮説に対する分散分析統計量，あるいは尤度比検定統計量を T (i) とおき，その統計量を最も大きくする遺伝子座 î = argmax T (i) の付近に QTL が存在すると判断するという手順である．
(t)

この手順に対応する統計モデルは次のようなものである．遺伝子型 zi= 1, 0, −1 に対応した 3 つの群の形質の平均を 3 つのパラメータによってμ − δ,

μ + α + δ,

μ−α+δと表現する．QTL は 1 つだけどこかの座 (i とおく) に存在すると仮定すると，モデルは以下のようになる．

(2.5)

1 つの i ∈ {1, . . . , m} が存在し，
(t)

(t)

y (t) = μ + αzi + δw(zi ) + ε(t) ,ただし

w(z) =⎧⎨1

ε(t) ∼ N (0, σ 2 ) (t = 1, . . . , n).

(z = ±1),⎩−1 (z = 0)

とおいた．
このモデル (2.5) に含まれる未知パラメータは (i, μ, α, δ, σ 2 ) である．パラメータ α, δ はそれぞれ QTL の加法効果，優性効果と解釈されている．QTL の効果がない (α = δ = 0)
ときは，QTL の位置パラメータ i は推測不能な量となる．この例のようにパラメータが特別な値をとるときにモデルの識別性が崩れるモデルを特異モデルという (福水，他 (2004))．

8

QTL の位置 i が既知であるという仮定をおくと，モデルは正則となる．その仮定の下で，パラメータ (μ, α, δ, σ 2 ) の最尤推定量を (μ̂(i), α̂(i), δ̂(i), σ̂ 2 (i))，また QTL が存在しな
い (QTL の効果がない) という帰無仮説 α = 0, δ = 0 の下での最尤推定量を (μ̃, 0, 0, σ̃ 2 )とおく．サンプルサイズ n を明示する形で，尤度関数を Ln と書くとき，座 i が QTL であるという仮定の下でのロッドスコア LODn (i) および尤度比検定統計量 LRTn (i) はLODn (i) = log10

Ln (μ̂(i), α̂(i), δ̂(i), σ̂ 2 (i))= 0.217 LRTn (i)Ln (μ̃, 0, 0, σ̃ 2 )である．(ロッドスコアは尤度比の常用対数として定義される．すなわち尤度比検定統計量の (2 log 10)−1 = 0.217 倍である．) QTL の位置 i も未知とするモデル (2.5) の下で，最尤推定量は (μ̂(î), α̂(î), δ̂(î), σ̂ 2 (î))，ただしî = argmax LODn (i)であり，単一マーカー分析における QTL の推測手順がモデル (2.5) の下での i の最尤推定量を与えることが分かる．
以降では，多重性調整のための max1≤i≤m LRTn (i) の分布計算の準備として，QTL が存在しないという帰無仮説の下での LRTn (i) (i = 1, . . . , m) の同時漸近分布を与える．尤
度比検定の一般論より，帰無仮説の下では個体数 n についての漸近的性質として，各 i に対して LRTn (i) は漸近的に自由度 2 のカイ 2 乗分布に従う．しかしそれらは独立ではない．カイ 2 乗確率変数の間の相関構造は以下のように表わされる．
命題 2.1 座間 i, j の組み換え価を 12 (1 − ρij ) とおく．QTL が存在しないという帰無仮説

H0 の下で，i = 1, . . . , m の同時分布の意味で分布収束
LRTn (i) ⇒ Ti = Ui2 + Vi2(2.6)(n → ∞)

が成り立つ．ただし (U1 , V1 , . . . , Um , Vm ) は平均 0 の 2m 次元正規分布ベクトルでCov(Ui , Uj ) = ρij ,Cov(Vi , Vj ) = ρ2ij ,Cov(Ui , Vj ) = 0

をみたすものである．とくに各 i について Ti は自由度 2 のカイ 2 乗分布に従う．
(t)= i = ±1 を母由来の相同染色体= δi = ±1 を父由来のそれとする．このとき証明 しばしば個体を識別する添字 (t) を省略する．i(t)の第 i 座の遺伝子型，δi

(t)

zi

= zi =

1(i + δi ),2(t)w(zi ) = wi = i δiと表わすことができる．

2 つの m 次元ベクトル (1 , . . . , m ), (δ1 , . . . , δm ) ∈ {−1, 1}m は独立で，それぞれの各成分は平均 0，分散 1，また組換えに由来する相関構造
E[i j ]= E[δi δj ]==1 × P (i = j ) + (−1) × P (i = j )1
1
(1 + ρij ) − (1 − ρij ) = ρij22
9を持っていた．zi , wi の 1, 2 次モーメントは E[zi ] = 0, E[wi ] = 0,

Cov(zi , zj ) = (E[i j ] + E[δi δj ])/4 = ρij /2,
Cov(wi , wj )= E[i δi j δj ] = E[i j ]E[δi δj ] = ρ2ij ,Cov(zi , wj )=(E[i j δj ] + E[δi j δj ])/2 = 0である．
一般性を失うことなくパラメータの真値を μ = 0, σ 2 = 1 とおく．
⎛
⎛
⎞
⎞
..
..
..
. ⎟
⎜ . ⎟
⎜ .
⎜
⎜
⎟
⎟
y = ⎜y (t) ⎟Xi = ⎜zi(t) w(zi(t) )⎟
⎝
⎝
⎠
⎠
..
..
..
.
.
.
1≤t≤n

1≤t≤n

とおく．テイラー展開より

LRTn (i) ≈ y T QXiT (XiT QXi )−1 XiT Qy,

Q = In − 1n 1Tn /n, 1n = (1, . . . , 1)T .

ここで ≈ は両辺の差が op (1) であること意味する．さらに
 (t)


n
1
1 T
1  zi − z̄i(t)(t)(zi − z̄i , wi − w̄i ) ≈ 2
Xi QXi =(t)nn0wi − w̄i
t=1


0
1

に注意すると

2 1 
22  (t)(t)(t)(y − ȳ)zi+(y (t) − ȳ)wi≈ u2i + vi2 ,
n t=1
n t=1
n

LRTn (i) ≈
ただし


ui =

n

2  (t) (t) 
≈
y zi
n t=1
n



n

1  (t) (t)  √≈ n δ̂(i)y wivi = √n t=1

nα̂(i) ,
2である．ここで ui , vi の平均は 0, 共分散関数はCov(ui , uj ) = 2 Var(y) Cov(zi , zj ) = ρij , Cov(vi , vj ) = Var(y) Cov(wi , wj ) = ρ2ij ,√Cov(ui , vj ) = 2 Var(y) Cov(zi , wj ) = 0.
あとは中心極限定理による．

注 2.2 戻し交配の場合は，i = 1, . . . , m の同時分布の意味でLRTn (i) ⇒ Ui2(n → ∞).

注 2.3 命題 2.1 は，組換え価がホールデンの地図関数 (2.4) であることは仮定していない．

(2.4) の仮定の下では，i, j 座間の遺伝的距離を rij とするとき，ρij = e−2rij であるので
ρij = ρi,i+1 ρi+1,i+2 · · · ρj−1,j .
すなわち (Ui )i≥1 および (Vi )i≥1 はマルコフ性を持った正規分布変数列となる．

10

2.4

分離比の検定

ここで述べる分離比の検定は，QTL 解析ではないが，数理的には前節の単一マーカー
分析と非常に似た構造を持つ連鎖解析である．

F2 個体の 1 つの遺伝子座に着目する．対立遺伝子を A, a とする．(2.2) によると，それが共優性である限り，遺伝子型は分離比の期待比率

A/A : A/a : a/a = 1 : 2 : 1を持つはずである．(A/a と a/A は区別されていない．) しかし実際の観測値は，理論値である 1:2:1 の比率を持った母集団からのサンプルとはみなされないことがある．このような現象は分離のゆがみとよばれる．そのような現象が起こる理由として，その近くに致死遺伝子 (生殖隔離障壁) が存在することが考えられる．生殖隔離障壁とは，それが特定の遺伝子型をとったときに生殖率，稔性が下がるような遺伝子をいう．Harushima, et al.

(2001) は，イネの F2 集団を用いて，そのような分離のゆがみを伴う生殖隔離障壁を検出している．
いま，分離のゆがみの検出のために，表 1，表 2 のデータが利用可能であるとする．た
だし表 1 の表現型のデータ y (t) はここでは不要である．
分離比が理論値に従っているかどうかを検定するためには，多項分布のカイ 2 乗適合度検定を用いることができる．各遺伝子型の個体数の観測度数を nA/A , nA/a , na/a とおく．分離比の検定統計量

Tn =(nA/A − n/4)2(nA/a − n/2)2(na/a − n/4)2
+
+
n/4
n/2
n/4

(n = nA/A + nA/a + na/a )は，分離比が 1 : 2 : 1 であるという帰無仮説の下で，自由度 2 のカイ 2 乗分布を漸近分布に持つ．
分離のゆがみを検出するためには，m 個のマーカー遺伝子 (i = 1, . . . , m) について，このカイ 2 乗検定を同時に行う必要がある．そのために，ここでも多重性の調整が必要となる．第 i 座における分離比の検定統計量を Tn,i とする．最大値 max1≤i≤m Tn,i の分布を近似する準備として，Tn,i (i = 1, . . . , m) の同時漸近分布を求めよう．
次のカイ 2 乗統計量の分解に注意する．

Tn =(nA/A − n/4)2(nA/a − n/2)2(na/a − n/4)2++= Un2 + Vn2 ,
n/4
n/2
n/4

ただし


Un =

2
(nA/A − na/a ),n1
Vn = √ (nA/A + na/a − nA/a ).
n

ここで，もし A/a と a/A が識別可能 (相が既知) とすると，頻度のデータは表 3 の形に集計される．また Vn = √1n (nA/A + na/a − nA/a − na/A ) である．1 つの個体 (t 番目とする)の，この表に対する寄与を考える．4 つのセルのうちの 1 か所で 1 回カウントされている(t)

はずである．相同染色体の母方染色体の第 i 座の遺伝子型 i= ±1 と，父方の対応する(t)遺伝子型 δi = ±1 に立ち返ると，この個体の，表 3 への寄与は表 4 となる．表 4 において，1 つのセルは 1，他の 3 つのセルは 0 である．

11

表 3: 相が既知の場合

A
a

A

a

nA/A
na/A

nA/a
na/a
n

表 4: 個体 t の表 3 への寄与

A
a

A

a

(t)
(t)
1
4 (1 + i )(1 + δi )
(t)
(t)
1
4 (1 − i )(1 + δi )(t)1
2 (1 + δi )

(t)
(t)
1
4 (1 + i )(1 − δi )
(t)
(t)
1
4 (1 − i )(1 − δi )(t)1
2 (1 − δi )(t)1
2 (1 + i )(t)1
2 (1 − i )122
このことから，第 i 座の分離比の統計量は Tn,i = Un,i
+ Vn,i
，ただし

Un,i

=

  n
n


2 1
1
(t)
(t)
(t)
(t)
(1 + i )(1 + δi ) −(1 − i )(1 − δi )n t=1 44
t=1=

1  (t)√u ,n t=1 i
n

Vn,i

=

n
n

1
1  1
(t)
(t)
(t)
(t)
√
(1 + i )(1 + δi ) +(1 − i )(1 − δi )4
n t=1 4
t=1

−

n

1
t=1

4

1  (t)√v ,n t=1 i
n

=

1 (t)
(t)
(t)
ui = √ (i + δi ),
2

(t)
(t)
(1 + i )(1 − δi ) −

n

1
t=1

(t)

vi

4


(t)
(t)
(1 − i )(1 + δi )(t) (t)= i δiと分解できる．

2 つの座 i, j に着目する．座間の組み換え価を 12 (1 − ρij ) とすると，前節ですでに計算(t) (t)(t) (t)(t)(t)
したように，E[i j ] = E[δi δj ] = ρij . したがって，E[ui ] = 0, E[vi ] = 0,
(t)

(t)

Cov(ui , uj ) = ρij ,

(t)

(t)

Cov(vi , vj ) = ρ2ij ,

(t)

(t)

Cov(ui , vj ) = 0が成り立つ．以上から中心極限定理によって，次が従う．
命題 2.2 座間 i, j の組換え価を 12 (1 − ρij ) とおく．i = 1, . . . , m の同時分布の意味で，分布収束

Tn,i ⇒ Ti(n → ∞)

が成り立つ．ただし Ti は命題 2.1 の (2.6) で定義したものである．

12つまり単一マーカー分析におけるロッドスコアの同時漸近分布と全く同じものが現わ
れる．

2.5

エピスタシス，遺伝子座相互作用の検出

今までは，QTL はたかだか 1 つ存在するというモデルを扱ってきた．しかし量的形質
は，複数の QTL によって引き起こされるものと考えられているため，そのようなモデル
では不十分である．とくにエピスタシスとよばれる QTL 間の交互作用を検出するためには，複数の QTL の存在を仮定した統計モデルを使う必要がある．
たとえば 2 つの QTL の存在を仮定した場合，単一マーカー分析のモデル (2.5) に対応
するものとして，次のモデルが考えられる．

i, j ∈ {1, . . . , m} が存在し，
(t)

(t)

(t)

(t)

y (t) = μ + α1 zi + δ1 wi + α2 zj + δ2 wj(t) (t)(t)(t)(t) (t)(t)(t)+β1 zi zj + β2 zi wj + β3 wi zj + β4 wi wj+ε(t) ,

ε(t) ∼ N (0, σ 2 )

ただし
(t)
(t)
wi = w(zi ) =⎧⎨1(t = 1, . . . , n),(t)= ±1),⎩−1 (z (t) = 0).
i(zi

座 i と座 j が QTL であるという仮定の下で，エピスタシスが存在しないという帰無仮説 β1 = · · · = β4 = 0 の尤度比検定を考える．サンプルサイズ n に対するロッドスコア (尤度比検定統計量) を，LRTn (i, j) と書く．多重性調整のため，エピスタシスが存在しないという帰無仮説の下での LRTn (i, j) (i, j = 1, . . . , m) の同時漸近分布を与えたい．尤度比
検定の一般論より，帰無仮説の下では個体数 n についての漸近的性質として，各 i, j に対して LRTn (i, j) は自由度 4 のカイ 2 乗分布に従う．しかしそれらは独立ではない．その相関構造は複雑なものとなるが，2 座 i, j が同じ染色体の上にない場合は，次のような直積型構造であることが示される．
命題 2.3 座間 i, j の組換え価を 12 (1 − ρij ) とおく．座 i と座 j が同じ染色体上にある (ない) ことを i ∼ j (i  j) と書く．エピスタシスが存在しないという帰無仮説の下で，全ての(i, j) ∈ {(i, j) | 1 ≤ i, j ≤ m, i  j}についての同時分布の意味で，分布収束
22LRTn (i, j) ⇒ Tij = U1,ij+ · · · + U4,ij(n → ∞)

が成り立つ．ただし (Uk,ij ) は k が異なると互いに独立な平均 0 の正規分布の配列でCov(Uk,ij , Uk,i j  ) = ρii ρjj  (k = 1), ρ2ii ρjj  (k = 2), ρii ρ2jj  (k = 3), ρ2ii ρ2jj  (k = 4)
をみたすものである．

13

春島，他 (2006) は，遺伝子座相互作用によって引き起こされる生殖隔離障壁を検出するために，2 座 i, j の遺伝子型の組合せとして得られる 3 × 3 表の独立性検定を行い，相互作用を独立性の乖離として検出することを試みた．その検定統計量の同時分布は帰無仮説の下で命題 2.3 と同じであることを示すことができる．

2.6

区間マッピング法と Haley-Knott の回帰分析

単一マーカー分析では，各マーカー遺伝子座についてロッドスコアが計算された．ここで説明する 2 つの方法は，ロッドスコアを補完によってマーカー遺伝子座の位置以外でも定義するものであり，マーカー間隔が密でない場合に有効である．
ここでも F2 集団で考える．QTL がある位置に 1 つ存在してそれが形質に影響を与えるというモデルを考える．このような仮想 QTL のことを putative QTL という．個体 t の(t)

QTL の遺伝子型を z∗ とおく．これは 1, 0, −1 の値をとる潜在変数である．さらにこれは，(t) (t)(t)(t)(t)
相同染色体上の母方，父方の遺伝子型 ∗ , δ∗ (= ±1) によって，z∗ = 12 (∗ + δ∗ ) と表される．
以下では混乱のない限り個体の添字 (t) を省略する．仮想 QTL の位置を γ とする．いま di ≤ γ ≤ di+1 ，つまり QTL はマーカー遺伝子座 i と i + 1 の間に存在するとする．

 = (1 , . . . , m ) と ∗ の同時分布，δ = (δ1 , . . . , δm ) と δ∗ の同時分布は，(2.1) と同様にマ
ルコフ性によって γ の関数として陽に書き下すことができる．このことから，z = (z1 , . . . , zm ) =1
1( + δ) = (1 + δ1 , . . . , m + δm )22が与えられたときの，QTL の遺伝子型 z∗ の条件付き分布が以下のように与えられる．

P (z∗ | z ; γ) =(2.7)

P (z, z∗ ; γ)
,
P (z)

ただし



P (z, z∗ ; γ) =z=(+δ)/2, z∗ =(∗ +δ∗ )/2

1

m−1



1 + j j+1 e−2(dj+1 −dj )

1 + δj δj+1 e−2(dj+1 −dj )2m+1
j=1, j=i
× 1 + i ∗ e−2(γ−di ) 1 + ∗ i+1 e−2(di+1 −γ)
× 1 + δi δ∗ e−2(γ−di ) 1 + δ∗ δi+1 e−2(di+1 −γ) ,(2.8)

P (z) =

z=(+δ)/2

di ≤ γ ≤ di+1 ,

m−1


1
−2(dj+1 −dj )−2(dj+1 −dj )
1
+

1
+
δ
.

e
δ
e
j
j+1
j
j+1
2m j=1(2.7) の P (z∗ | z ; γ) は，γ → di または di+1 のとき zi または zi+1 に確率 1 で値をとる
一点分布となる．P (z∗ | z ; γ) は γ の連続関数であるが，マーカー点で滑らかではない．
仮想 QTL の遺伝子型の情報は，m 個のマーカーのなかで，とくに QTL に隣接するマーカー (ﬂanking marker) i と i + 1 が多く持っていると考えられる．そのため QTL の遺伝

14

子型の予測のために (2.7) の代わりに P (z∗ | zi , zi+1 ; γ) を考えることもできる．これは簡便法であるが，戻し交配の場合のように，マーカー遺伝子型自体にマルコフ性が成り立つ
ばあいには， (2.7) と正確に一致する．

Lander & Botstein (1989) の区間マッピング法 (interval mapping) とは，次の統計モデ
ルを仮定した解析法である．
(t)

(t)

(2.9)

z∗

∼ P (z∗ | z (t) ; γ),(2.10)y (t)= μ + αz∗ + δw(z∗ ) + ε(t) ,

(t)

(t)

(t = 1, . . . , n)，ただしw(k) =⎧⎨1

ε(t) ∼ N (0, σ 2 )(k = ±1),⎩−1 (k = 0).

また Haley & Knott (1992) は，区間マッピング法の簡便法として，潜在変数をその期待値で置きかえた次の回帰分析を提案した．

(2.11)

y (t)

(t)

(t)

= μ + αE[z∗ | z (t) ; γ] + δE[w(z∗ ) | z (t) ; γ] + ε(t) ,ε(t) ∼ N (0, σ 2 )(t = 1, . . . , n).

両者において，仮想 QTL の位置 γ の関数として尤度関数が定義されるので，連続な曲線としてロッドスコアが定義される．
また両者ともに論文で提案されたオリジナルの形は，QTL の遺伝子型の分布として隣接マーカーの情報だけを用いた P (z∗ | zi , zi+1 ; γ) を仮定するものであるが，本稿では全マーカーの情報を用いた P (z∗ | z ; γ) で考えることにする．
以下では，区間マッピング法の尤度関数の形を書き下し，ロッドスコアの確率過程としての構造を調べていくことにする．Haley-Knott の回帰分析の場合は最後に触れる．
(t)

マーカーの遺伝子型 z (t) と QTL の位置 γ が与えられたとき z∗ は 3 値の離散分布に従う．記法を簡単にするために，その確率を (2.7) を使って
(t)

(t)

πk (γ) = P (z∗ = k | z (t) ; γ),k = −1, 0, 1とおくと，z (t) が所与のときの y (t) の分布は，コンポーネント数が 3 の正規分布の有限混
合分布
1


(t)

πk (γ)fk,θ (y (t) )k=−1となる．ただし正規分布 N (μ + αk + δw(k), σ 2 ) の密度関数を fk,θ (·), θ = (α, δ, μ, σ 2 ) とおいた．したがって，(z t , y (t) ) (t = 1, . . . , n) の同時密度関数は，(2.8) を用いて
n  
1


(t)
πk (γ)fk,θ (y (t) )P (z (t) )t=1 k=−1と書ける．
ロッドスコアを描くためには，各 γ についてこの尤度を θ について最大化する必要がある．そのために EM アルゴリズム (Wu (1983) など) を用いることができる．ここで扱う

15

モデルは，混合確率は個体 t に依存すること，また (γ を所与とすると) 混合確率には推測
対象の未知パラメータが含まれないこと，の 2 点で通常の有限混合モデルとはやや異なっ
ている．
(t)

z∗ と一対一に対応するダミー変数ベクトル
(t)

(t)

(t)

e(t) = (e−1 , e0 , e1 ),⎧⎨1 (z (t) = k)(t)ek =⎩0 (z (t) = k)を導入する．マーカー遺伝子の遺伝子型，表現型，ならびに潜在変数 (z (t) , y (t) , e(t) ) (t =1, . . . , n) の同時分布は
n 


e(t)

k
(t)
πk (γ)fk,θ (y (t) )P (z (t) )1

t=1 k=−1と書けるので(t)π (γ)fk,θ (y (t) )(t)Eθ ek | (z (t) , y (t) )t=1,...,n = 1 k (t)(t)
k=−1 πk (γ)fk,θ (y )となる．これを用いて，EM アルゴリズムは次のようにまとめられる．

1. 初期値(t)(t)

êk := πk (γ)(k = −1, 0, 1; t = 1, . . . , n).

(t)

2. 以下を θ̂ と êk が収束するまで繰りかえす．
n1θ̂ := argmaxθê(t)k(t)
πk (γ)fk,θ (y (t) ),t=1 k=−1(t)πk (γ)fk,θ̂ (y (t) )(t)êk := 1(t)(t)
k=−1 πk (γ)fk,θ̂ (y )(k = −1, 0, 1; t = 1, . . . , n).

ステップ 2 において θ̂ は，Qn (θ ; γ) =n 1
t=1 k=−1

(t)

êk

1
σ


(t)
2
2
(y
−
μ
−
αk
−
δw(k))
+
log
σ
+
(
定数
)
2

を θ = (α, δ, μ, σ 2 ) について最小化する操作 (重み付き最小 2 乗法) で陽に求めることができる．

1 章の図 1 は区間マッピング法により描いたものである．同じものを第 3 染色体について拡大して描いたものが図 3 である．ロッドスコア曲線が，マーカー間を補完していることが分かる．
注 2.4 ここで区間マッピング法におけるモデリングを振りかえってみる．統計モデルは，連鎖を記述する部分 (2.9) と，形質発現を記述する部分 (2.10) で構成されていた．後者の形質発現モデルを工夫することによって，QTL が複数ある場合や，表現型が非正規分布
や離散分布，あるいは多変量分布に従う場合などを扱うことができる．このような統一的な扱いは，Sen & Churchill (2001) により提案され，R/qtl (Broman, et al. (2003)) として実装されている．

16

5lod.Log_ADI

4

3

2

1

0
0

20

40

60

80

Map position (cM)

図 3: 区間マッピング法によるロッドスコア (第 3 染色体)2.7

区間マッピング法のロッドスコア

区間マッピング法のロッドスコア (の定数倍である) LRTn (γ) は，各 γ について尤度比
検定統計量であり，QTL が存在しないという帰無仮説の下で漸近的に自由度 2 のカイ 2 乗
分布に従う．ここでは，QTL が存在しないという帰無仮説の下で，LRTn (·) を確率過程(漸近的カイ 2 乗確率過程) とみなしたときの相関構造を確定する．これは多重性調整のために必要である．
対数尤度は

L(γ)
n (θ) =

n

t=1

log

1


(t)

πk (γ)fk,θ (y (t) ) + (θ を含まない項),

θ = (α, δ, μ, σ 2 )

k=−1

であった．
パラメータ θ を θ = (θ1 , θ2 )，ただし θ1 = (α, δ), θ2 = (μ, σ 2 ) と分割表現する．帰無仮説は H0 : θ1 = 0 である．一般性を失わずに，真値を θ0 = (θ10 , θ20 ) = (0, 0, 0, 1) とする．

γ を固定し，θ に関するフィッシャー情報行列を ∂L(γ)  
I11 I12n = Varθ0(2.12)
I(θ0 ; γ) =∂θ θ0
I21 I22
とおく．尤度比検定の一般論から，真値 θ0 の下で (γ)(γ)
 ∂L(γ)nI∂Ln ∂Ln −1−1∂θ1
I11·2 I, −I12 I22
LRTn (γ) ≈
,
−1
∂L(γ)
n
∂θ1
∂θ2
I21
−I22
(θ10 ,θ20 )∂θ2−1
ただし I11·2 = I11 − I12 I22
I21 である．

17

帰無仮説 α = δ = 0 のとき，fk,θ は k によらない．このことから，スコア関数を真値で
評価すると
(γ)
∂Ln 
∂θ θ0(2.13)

=

n




(t)
(t) ∂
(t) 
k πk (γ)fk,θ (y ) ∂θ log fk,θ (y ) 

 (t)
(t)
θ0
k πk (γ)fk,θ (y )

t=1

=

=


∂

log fk,θ (y (t) )∂θ
θ0
t=1 k=−1⎞⎛(t)y (t) k kπk (γ)
⎟
⎜

n

⎜y (t) k w(k)πk(t) (γ)⎟
⎟
⎜
⎟
⎜
(t)
y
⎠
⎝
t=1
(t)2
(y
− 1)/2
1
n 


(t)

πk (γ)(γ)(γ)

である．Ln (θ0 ) は γ に依存しない (特異モデルであるため) が，∂Ln /∂θ|θ0 は γ に依存することに注意する．
スコアベクトルの分散共分散関数をR(γ, γ̃) = ∂L(γ)  ∂L(γ̃)  
1
n 
n 
Covθ0
 ,

n
∂θ θ0 ∂θ θ0

とおく．
1


1


 (t)

(t)
kπk (γ) = E z∗ | z (t) ; γ ,

k=−1



(t)
(t)
w(k)πk (γ) = E w(z∗ ) | z (t) ; γk=−1に注意すると，簡単な計算により⎛⎜
R(γ, γ̃) = ⎝R11 (γ, γ̃)
O

O

⎞

⎟
1 0 ⎠
0 1/2 4×4

ただし(2.14)

R11 (γ, γ̃)


 
!
(t)


(t)
(t)
k kπk (γ)= E (t)k kπk (γ̃),
k w(k)πk (γ̃)
k w(k)πk (γ)
!
E[z∗ | z ; γ]E[z∗ | z ; γ̃], E[w(z∗ ) | z ; γ̃]= E
E[w(z∗ ) | z ; γ]が分かる．ここで外側の期待値はマーカーの遺伝子型 z についてとる．フィッシャー情報
T
行列 (2.12) は I(θ0 ; γ) = nR(γ, γ) であり，ブロック対角行列 (I12 = I21= 0) となる．

C(γ) を各成分が γ について滑らかな 2 × 2 行列で C(γ)R11 (γ, γ)C(γ)T = I2 をみたすものとする．中心極限定理より以下が示される．
命題 2.4 有限個の γ の，有限次元周辺分布の分布収束の意味で

LRTn (γ) ⇒ T (γ) = U (γ)2 + V (γ)218(n → ∞).

ただし，U (·), V (·) は平均 0 の正規過程で，その共分散関数は (2.14) の R11 (γ, γ̃) を用いて
 


U (γ)
U (γ̃)
UV(2.15)R (γ, γ̃) = Cov,= C(γ)R11 (γ, γ̃)C(γ̃)TV (γ)
V (γ̃)
として記述される．とくに固定した γ について LRTn (γ) は漸近的に自由度 2 のカイ 2 乗
分布に従う．
注 2.5 上の命題は，確率過程としての収束は述べていない．確率過程としての収束は，尤度比確率場の弱収束に関する Ibragimov & Has’minskii (1981) の方法で示すことができる(Yoshida (2006), Introduction を参照)．
注 2.6 Haley-Knott の回帰分析 (2.11) のスコア関数を真値 θ0 で評価すると，(2.13) と同じ形となる．ロッドスコアは区間マッピング法と同じ極限分布を持つ．

3

多重性調整のための方法

3.1

経験的方法とシミュレーション

前章ではいろいろな QTL 解析のモデルについて，QTL が存在しないという帰無仮説のもとで，ロッドスコアの確率構造を確定した．その結果を出発点として，本章では多重性調整のために必要な，ロッドスコアの最大値の分布の近似法について説明する．そのために利用可能な数学的な方法として，非線形再生理論やオイラー標数法がある．それらについては，後の 3.2, 3.3 節で説明することとし，本節では QTL 解析の多重性調整のために行われている経験則とシミュレーションによる方法を紹介する．

Lander & Kruglyak (1995) は，ロッドスコアの多重性調整の目安として，多重性未調
整の p 値を表 5 に従って解釈することを提唱している．しかしゲノムワイドの多重性調整
は，染色体長のみならずマーカーの密度に大きく依存する．(3.2 節，表 4 の数値実験を参照．) そのため，機械的にこの表を用いることは行われてはいない (石川 (2006) など)．
表 5: 多重性未調整 p 値の解釈
実験交配の手法

suggestive

signiﬁcant

BC (1 d.f.)
F2 (1 d.f., 加法効果)−3

3.4 × 10
3.4 × 10−3

1.0 × 10−4
1.0 × 10−4

F2 (1 d.f., 優性効果)
F2 (2 d.f)

2.4 × 10−3
1.6 × 10−3

7.2 × 10−4
5.2 × 10−5

多重性調整の方法として現在広く用いられている方法は，Churchill & Doerge (1994) の提案による並べ替え検定である．その手順は次のようなものである．

1. 個体数を n とし，集合 {1, . . . , n} の置換の全体を Πn とおく．
19

2. N を十分大きな数とし，以下の手順を N 回繰りかえす．その繰り返しを k = 1, . . . , Nとする．

(i) ランダムに π ∈ Πn を選ぶ．
(ii) 表現型 (y (t) )t=1,...,n を，置換 π によって入れ替えたデータセット (表 6) について
QTL 解析を行い，ロッドスコアの最大値を数値的に探索する．それを MaxLODk
とおく．

3. {MaxLODk }k=1,...,N の経験分布をロッドスコアの最大値の経験分布とみなして p 値の推定値を#{k | MaxLODk ≥ MaxLOD (実現値)}値=p"Nと計算する．

表 6: 並べ替え検定のためのデータセット
個体番号

表現型

1
..
.

y

(π(1))ny (π(n))

遺伝子型

z

..
.

(1)

(1)

(1)

(n)

(n)

= (z1 , . . . , zm )..
.

z (n) = (z1 , . . . , zm )この並べ替え検定は直感的に分かりやすい方法であり，R/qtl などの多くのプログラムに実装されている．しかし次のような問題点がある．
最初の問題点は，並べ替え検定の一般論に関わる問題である．並べ替え検定が生成する{MaxLODk }k=1,...,N の経験分布は，帰無仮説が成り立つ場合とそうでない場合とでは当然異なるものである．そのために，帰無仮説が正しい場合には閾値を正しく推定できたとしても，帰無仮説が正しくない場合には閾値を過大評価する可能性がある．そのときは，ピークの検出確率 (検出力) の低下を招くことになる．
もうひとつの問題点は，計算量の問題である．説明した並べ替え検定の手順では，ロッドスコアの最大値を数値的に探索する必要があるため，計算時間がかかる．QTL を 1 つしか想定しない場合は問題ないが，エピスタシスの検定などで複数の QTL を想定する場
合は，探索すべき組合せ数が莫大となり，計算が実行可能でなくなる場合がある．
並べ替え検定は汎用的でしばしば用いられる手法であるが，これらの理由のためその解釈や利用に注意を払う必要がある．可能な限り，他の多重性調整の方法の結果と併用するのがよいと思われる．問題を単一マーカー分析あるいは分離比の検定における多重性調整に限定すると，AR モデルを用いたモンテカルロシミュレーションも利用可能である．

1. 1 , . . . , m , δ1 , . . . , δm を標準正規分布 N (0, 1) に従う i.i.d. 列とする．U1 = 1 , V1 = δ1
とおく．

20

2. i = 2, . . . , m について以下を計算する．
#
Ui = αi Ui−1 + 1 − αi2 i (αi = e−2(di −di−1 ) ),#Vi = βi Vi−1 + 1 − βi2 δi (βi = e−4(di −di−1 ) ),Ti

= Ui2 + Vi2 .

3. max1≤i≤m Ti を計算する．
上の手続きを十分な回数繰り返すことによって，命題 2.1, 2.2 の max1≤i≤m Ti の経験分布を求めることができる．しかしながら，並べ替え検定のときと同様，複数の QTL を想定した場合には計算量が膨大となる．

3.2

非線形再生理論による近似

本節以降では，ロッドスコアの最大値の分布を理論的に求める (近似する) ための方法
を 2 つ紹介する．最初に述べる方法は，逐次解析，非線形再生理論を用いるものであり，
単一マーカー分析や分離のゆがみの検定において，マーカーが密で間隔が等間隔に近い場
合 (|di+1 − di | ≈ Δ  1) に有効な方法である (Dupuis & Siegmund (1999))．
まず一般的な形で問題設定を行う．Zk (t), t ∈ I ⊂ R (k = 1, 2) を互いに独立な正規過程で，平均 0，共分散関数がCov(Zk (t), Zk (t̃)) = e−ρk |t−t̃|(ρk > 0)であるもの (Ornstein-Uhlenbeck 過程) とする．自由度 2 のカイ 2 乗確率過程の平方根をY (t) =$Z1 (t)2 + Z2 (t)2(t ∈ I)で定義する．
命題 3.1 格子間隔を Δ > 0 とし，J = {j ∈ Z | jΔ ∈ I} とおく．b → ∞, Δ → 0,√
b Δ → c (> 0) のとき

(3.1)

% 2π


$
2
dθ
P max Y (jΔ) ≥ b ∼ |I| b2 e−b /2ρ̄(θ) ν(c 2ρ̄(θ)).
j∈J
2π
0

ここで | · | は集合のルベーグ測度，ρ̄(θ) = ρ1 cos2 θ + ρ2 sin2 θ．また Φ(·) を標準正規分布の分布関数とするとき⎧
 ∞
 1 √ ⎪⎪⎨2x−2 exp −2
n−1 Φ − x n(x > 0),2
ν(x) =
n=1⎪⎪⎩1(x = 0).
注 3.1 x ≤ 2 くらいでは ν(x) ≈ exp(−0.583x) と近似できる．実用的にはこの範囲で足りることが多い．

21

注 3.2 形式的に Δ = 0, c = 0 とおいて得られる式

2
P sup Y (t) ≥ b ∼ |I|Kb2 e−b /2 ,

% 2π



t∈I

K=
0

dθ
1
ρ̄(θ) = (ρ1 + ρ2 ) (b → ∞)2π2も成り立つ (Piterbarg (1996), Corollary 7.1 より)．

I = [0, 1], Δ = 0.01 (染色体長 1M, マーカー数 m = 100 に相当) の場合の数値例
を図 4 に示す．一点鎖線 (− · −) と実線 (—) はそれぞれ格子点上の最大値の上側確率
P maxj∈J Y (jΔ) ≥ b の命題 3.1 による近似値とシミュレーションによる推定値を表している．すなわち図の縦軸は p 値に対応する．ただしここでは最大値分布の近似として (3.1)の右辺をそのまま用いるのではなく，それと漸近的に同等な 1 − exp{−(3.1) 右辺 } を用いている．この両者は p 値が 0.5 程度より小さい範囲で非常に精度良く一致していることが見て取れる．


また破線 (– –) は Δ = 0，すなわち連続集合 I 上の最大値分布 P supt∈I Y (t) ≥ b

である．Lander & Botstein (1989) は区間マッピング法を提案するとともに，Ornstein-

Uhlenbeck 過程の最大値の分布を用いてゲノムワイドな多重性調整を行うことを提唱している．それがこの場合に対応するが，マーカー間隔 10cM (Δ = 0.01) という密な場合であってもその近似はあまり良くない．
点線 (· · · ) は自由度 2 のカイ 2 乗分布の上側確率である．これは多重調整を行わない場
合の p 値に対応する．
命題 3.1 の証明は，章末で与える．なおここで与えた近似は，マーカー間距離が常に一定値 Δ であるというモデルに基づくものである．実際問題への適用の際には，Δ にマー
カー間距離の平均値を代入して用いることになる．Dupuis & Siegmund (1999) は類似の問題設定において，そのような場合であってもよい近似を与えることを数値計算によって確認している．
また命題 2.3 で与えた 2 次元格子上のカイ 2 乗確率場の最大値についても，同様の近似を与えることができる (栗木 (2007))．

3.3

ランダム関数の零点の個数の期待値

次に，滑らかなサンプルパスを持つ確率場の最大値の上側裾確率を近似するためのオイ
ラー標数法 (Euler characteristic heuristic) を，添字集合が 1 次元という特殊な場合について説明する．この方法は，信号処理の分野でライスの公式 (Rice’s formula) として知られているものと同等である．区間マッピング法のロッドスコアは，漸近的には区分的に滑らかなサンプルパスを持つカイ 2 乗確率過程であったので，オイラー標数法によってその最大値の分布を近似することができる．オイラー標数法の一般論については，栗木・竹村(2007) を参照のこと．
Z(t), t ∈ I ⊂ R は，実数に値をとるランダムな C 1 関数で，各 t について (Z(t), Ż(t))(Ż(t) = dZ(t)/dt) が縮退しない分布を持つとする．Z(t) = u となる t (方程式 Z(t) − u = 0の零点) の個数をNu = #{t ∈ I | Z(t) = u}
22

1.0
0.6
0.4
0.0

0.2

upper probab.

0.8

simulation (10000 rep)
asymptotic prob.
discrete correction
chi2 (df=2)

0

5

10

15

20

25

30

max chi2

図 4: 格子点上の最大値の分布( −·−：命題 3.1 による近似, —：シミュレーション, – –：連続近似, · · · ：χ22 分布)

とおく．
ところで，任意の連続関数 f に対して
%
%
(3.2)
|Ż(t)| f (Z(t)) dt = Nu f (u) duが成り立つことが容易に分かる (図 5)．

(3.2) の両辺について，Z(·) の期待値をとると%%

E[Nu ]f (u)du =E |Ż(t)| f (Z(t)) dt
% %


=
E |Ż(t)| | Z(t) = u f (u)pZ(t) (u) du dt
を得る．最右辺において，pZ(t) は Z(t) の周辺密度である．これが任意の連続関数 f で成り立つので，%(3.3)

E[Nu ] =
E |Ż(t)| | Z(t) = u pZ(t) (u) dta.s.

である (Azaı̈s & Wschebor (2005))．
いま閾値 u が大きいとする．このとき，関数 Z(t) は閾値 u を超えることは稀であり，もし一度超えてもまたすぐに閾値を下回ることが予想される．つまり Nu が 0 または 2 以外23

f (u)

Nu=0

Nu=2

Nu=1

u
t=

t=-

Z (t )

図 5: (3.2) の図による説明の値をとる確率が小さいことが期待される．そのような状況では

1
E[Nu ] ≈ P (Nu = 2)2
≈ P ∃t ∈ I, Z(t) ≥ u= P sup Z(t) ≥ ut∈Iである．オイラー標数近似とは，最大値分布の近似式として


1
(3.4)
P sup Z(t) ≥ u ≈ E[Nu ] (u が大きいとき)
2t∈Iとおき，(3.4) の右辺として，(3.3) を用いるというものである．ここではこのような直感的な議論にとどめるが，正規確率場，カイ 2 乗確率場に対するオイラー標数法は，正則条
件の下で非常に良い近似を与えることが知られている．(栗木・竹村 (2007) の参考文献を参照．)
このオイラー標数法によって，F2 集団の区間マッピング法によるロッドスコアの最大値の分布を近似しよう．仮想 QTL の位置を γ とする．区間マッピング法で補間されたロッドスコア (尤度比検定統計量) LRTn (γ), γ ∈ Γ の，QTL が存在しないという帰無仮説のもとでの漸近分布は，連続かつ隣り合うマーカーの間では滑らかなパスを持つカイ 2 乗確
率過程 T (γ), γ ∈ Γ であった．全ての隣り合うマーカー間について E[Nu ]/2 を計算し，それを全て足しあわせることによって最大値分布の近似が可能である (Rebaı̈, et al. (1994))．
命題 3.2 RU V (γ, γ̃) を，命題 2.4 で与えた共分散関数 (2.15) とする．

A(γ) =

∂ UV
,R (γ, γ̃)
∂γ̃
γ̃=γ

B(γ) =∂2

RU V (γ, γ̃)
∂γ∂γ̃
γ̃=γ

とおく．h = (h1 , h2 )T を円周上の一様分布 Unif(S1 ) とする．ロッドスコアの極限過程 T (·)の，Γ = [d1 , dm ] 上最大値のオイラー標数近似は

(3.5)



P sup T (γ) ≥ u
γ∈Γ



m−1 %
#

1 1/2 −u/2  di+1 ≈ √ u eEhT (B(γ) − A(γ)T A(γ))h dγ
2πi=1 diで与えられる．

24

この命題の証明も，本章の章末で与える．その導出は Davies (1987), Theorem A.1 と本質的に同じである．

3.4

命題の証明

3.4.1

命題 3.1 の証明

ここで与える証明は，Kim & Siegmund (1989) と Siegmund (1992) を組合せたものである．t を固定し，j ∈ Z を添字とみなした確率変数列
Y'j = b(Y (t + jΔ) − y)

(Z1 (t),Z2 (t))=(y cos θ,y sin θ)$を定義する．Y (t) =Z1 (t)2 + Z2 (t)2 = y に注意する．b, y → ∞, b ∼ y, Δ → 0,√
b Δ → c (> 0) の極限を考える．Y (t + jΔ) を t のまわりでテイラー展開し，また正規分
布の条件付分布の公式より，(Yj )j∈Z の有限次元の周辺分布は正規分布で，その平均，分
散はE[Y'j ]= −ρ̄(θ)c2 |j|,

Cov(Y'j , Y'j̃ ) = ρ̄(θ)c2 (|j| + |j̃| − |j − j̃|)⎧⎨2ρ̄(θ)c2 min(|j|, |j̃|) (j と j̃ は同符号),
=
⎩0
(異符号)
であることを確認できる．これは Xi (i ∈ Z) を独立な正規分布 N (−ρ̄(θ)c2 , 2ρ̄(θ)c2 ) の列としたときの，両側ランダムウォーク⎧⎪⎪
X + · · · + Xj⎪
⎨ 1
Sj = 0⎪⎪⎪⎩X + · · · + X−1(j > 0),(j = 0),−|j|(j < 0)の分布に等しい．すなわち(3.6)(. . . , Y'−1 , Y'0 , Y'1 , . . .) ⇒ (. . . , S−1 , S0 , S1 , . . .)である．
条件 (Z1 (j 0 Δ), Z2 (j 0 Δ)) = (y cos θ, y sin θ) を与えた条件付確率測度を P' で表わす．

−x = b(b − y) (すなわち y = b + x/b) とおく．まず (3.6) よりP' max0 Y (jΔ) < b= P' max0 b(Y (jΔ) − y) < −x
j>j
j>j




= P' max Y'j < −x → P max Sj < −xj>0に注意する．

25

j>0これからが証明の本体である．Y (jΔ) ≥ b をみたすような点 j の最大値を j 0 とおく．事
象 {maxj∈J Y (jΔ) ≥ b} は j 0 の値によって排反に分割される．



 
(3.7)
P max Y (jΔ) ≥ b =P max Y (jΔ) < b, Y (j 0 Δ) ≥ b .
j∈J

j>j 0

j 0 ∈JZ (j 0 Δ)

さらに事象を tan−1 Z21 (j 0 Δ) の値によって排反に分割する．

(3.7) 右辺 =% 2π 
0j 0 ∈J
Z2 (j 0 Δ)
P max0 Y (jΔ) < b, Y (j 0 Δ) ≥ b, tan−1∈dθZ1 (j 0 Δ)
j>j

% 2π 

%
=

y≥b

0

j 0 ∈J
Z2 (j 0 Δ)
P max0 Y (jΔ) < b, Y (j 0 Δ) ∈ dy, tan−1∈dθ,Z1 (j 0 Δ)
j>j

ただし dθ = (θ, θ + dθ), dy = (y, y + dy) である．

tan−1

Y (j 0 Δ) = y,

Z2 (j 0 Δ)
=θ
Z1 (j 0 Δ)と値を与えることは (Z1 (j 0 Δ), Z2 (j 0 Δ)) = (y cos θ, y sin θ) と値を与えることと等価なので，


Z2 (j 0 Δ)
P max0 Y (jΔ) < b, Y (j 0 Δ) ∈ dy, tan−1∈dθZ1 (j 0 Δ)
j>j




Z2 (j 0 Δ)∈dθ
= P' max0 Y (jΔ) < b × P Y (j 0 Δ) ∈ dy, tan−1
Z1 (j 0 Δ)
j>j


dθ
→ P max Sj < −x × P (Y (j 0 Δ) ∈ dy) ×
j>0
2π
(
(
である．ここで y = b + x/b と変数変換する． y>b = x>0 ならびに
2

2

P (Y (j 0 Δ) ∈ dy) = ye−y /2 dy ∼ e−b /2 e−x dx(b → ∞)に注意する．
さらに逐次解析で知られた関係式
% ∞


$
e−x P max Sj < −x dx = ρ̄(θ)c2 ν(c 2ρ̄(θ))
j>0

0

および

Δ

−1 2

c ∼b ,2Δ
j 0 ∈J%∼

dt(Δ → 0)Iを組合せると (3.1) を得る．

3.4.2

命題 3.2 の証明

仮想 QTL の位置を γ ∈ [di , di+1 ] とする．命題 2.4 の共分散関数 RU V (γ, γ̃) (2.15) は引数について十分に滑らかである．とくに ∂  3  ∂ 3

∂γ

∂γ̃
26

RU V (γ, γ̃)が γ = γ̃ の近傍で存在し有界なので，U (·), V (·) は確率 1 で連続微分可能で，また微分過
程 U̇ (·), V̇ (·) も連続な正規過程となる (Adler & Taylor (2007), Theorem 1.4.2)．
微分過程 U̇ (·), V̇ (·) も平均 0 で，その共分散関数は
⎞
⎛
U (γ)

⎟ 
⎜
⎜V (γ)⎟
I
A(γ)
2
⎟
Cov ⎜
⎜U̇ (γ)⎟ = A(γ)T B(γ)⎠⎝
V̇ (γ)
である．ここで

0

=

0

=

d
Cov(U (γ), U (γ)) = 2Cov(U̇ (γ), U (γ)),
dγd
Cov(U (γ), V (γ)) = Cov(U̇ (γ), V (γ)) + Cov(U (γ), V̇ (γ))
dγであるので A(γ)T = −A(γ) (交代行列) である．
以下では引数の γ は省略する．(U̇ , V̇ ) の，(U, V ) を与えた条件付き分布は

であり，また

 U

U̇ 
, B − AT A
∼ N AT
V̇ (U,V )V
U̇Ṫ = 2(U, V )
V̇


U
,
T = (U, V )
V

より，



Ṫ 




U
U 
T
, 4(U, V )(B − A A)
∼ N 2(U, V )A
V
V



U
= N 0, 4(U, V )(B − AT A)V(U,V )

Tである．したがって#= ξ ∼ N (0, 1)
U 
4(U, V )(B − AT A) V (U,V )Ṫ

となり，ξ の条件付き分布は条件 (U, V ) によらないことが分かるので，)
UT(ξ ∼ N (0, 1) は (U, V ) と独立)Ṫ = ξ 4(U, V )(B − A A)
Vと表わすことができる．
さらに

h=


h1
h2


1 U
=√
T V

√ $
とおくと，Ṫ = 2ξ T hT (B − AT A)h であり，h は ξ, T と独立に，円周上の一様分布
Unif(S1 ) に従う．
27

以上より



E |Ṫ | | T = u pT (u) ===E |Ṫ | 1{T ∈(u,u+du)} /du

#

√
hT (B − AT A)h
2E[|ξ|] E T 1{T ∈(u,u+du)} /du × E#
2 1/2 1 −u/2 EhT (B − AT A)h .
2
u
e
π
2

これを γ で積分して，

1
1
E[Nu ] = √ u1/2 e−u/222π% di+1E#hT (B(γ) − A(γ)T A(γ))h dγ,diただし右辺の期待値は h ∼ Unif(S1 ) についてとる．
これを全ての隣り合うマーカーについて足しあわせることによって (3.5) が得られる．

謝辞
1 章で紹介したマウスの QTL 解析例 (図 1, 3) は，国立遺伝学研究所の城石研究室 (前野哲輝氏，城石俊彦氏) によるものです．また同研究所の春島嘉章氏，倉田のり氏からは，
本稿の内容について有益なコメントをいただきました．本稿は，数理科学総合セミナー II(2006 年度前期，東京大学数理科学研究科）の講義録に加筆したものです．同大学の吉田朋広氏には，注 2.5 をご指摘いただいたきました．これらの方々に感謝いたします．

参考文献[1] Adler, R. J. and Taylor, J. E. (2007). Random Fields and their Geometry, Springer.
[2] Azaı̈s, J-M. and Wschebor, M. (2005). On the distribution of the maximum of a
Gaussian ﬁeld with d parameters, Ann. Appl. Probab., 15 (1A), 254–278.
[3] Broman, K. W., Wu, H., Sen, Ś. and Churchill, G. A. (2003). R/qtl: QTL mappingin experimental crosses, Bioinformatics, 19 (7), 889–890.
[4] Churchill, G. A. and Doerge, R. W. (1994). Empirical threshold values for quantitative triat mapping, Genetics, 138 (3), 963–971.
[5] Davies, R. B. (1987). Hypothesis testing when a nuisance parameter is present only
under the alternatives, Biometrika, 74 (1), 33–43.
[6] Dupuis, J. and Siegmund, D. (1999). Statistical methods for mapping quantitative
trait loci from a dense set of markers, Genetics, 151 (1), 373–386.
[7] 福水健次, 栗木哲, 竹内啓, 赤平昌文 (2004). 「特異モデルの統計学 — 未解決問題への新しい視点」, 統計科学のフロンティア 7, 岩波書店.

[8] Haldane, J. B. S. (1919). The combination of linkage values and the calculation of
distances between the loci of linked factors, J. Genetics, 8, 299-309.
28

[9] Haley, C. S. and Knott, S. A. (1992). A simple regression method for mapping quantitative trait loci in line crosses using ﬂanking markers, Heredity, 69, 315–324.
[10] 春島嘉章, 栗木哲, 水多陽子, 藤澤洋徳, 倉田のり (2006). 配偶体内または接合体内の異なる遺伝子座間の相互作用による生殖的隔離障壁の検出, 育種学研究, 8, 別冊 2 号,

157.
[11] Harushima, Y., Nakagahra, M., Yano, M., Sasaki, T. and Kurata, N. (2001). A
genome-wide survey of reproductive barriers in an intraspeciﬁc hybrid, Genetics,
159 (2), 883–892.
[12] Hochberg, Y. and Tamhane, A. C. (1987). Multiple Comparison Procedures, Wiley.
[13] Ibragimov, I. A. and Has’minskii, R. Z. (1981). Statistical Estimation: Asymptotic
Theory, Springer.
[14] 石川明 (2006). 動物モデルによる多因子性疾患の QTL 解析：その基礎的理論と解析
方法 (改訂版), NAGOYA Repository, http://hdl.handle.net/2237/6779[15] Karlin, S. and Liberman, U. (1983). Measuring interference in the chiasma renewal
formation process, Adv. Appl. Probab., 15 (3), 471–487.
[16] Kim, H.-J. and Siegmund, D. (1989). The likelihood ratio test for a change-point in
simple linear regression, Biometrika, 76 (3), 409–423.
[17] 栗木哲 (2007). 直積型の相関構造を持つカイ 2 乗確率場の最大値の分布, 日本数学会
2007 年度年会 統計数学分科会講演アブストラクト, 89–90.
[18] 栗木哲, 竹村彰通 (2008). チューブの体積と正規確率場の最大値の分布, 数学, 60 (2),
134–155.
[19] Lander, E. S. and Botstein, D. (1989). Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps, Genetics, 121 (1), 185–199.
[20] Lander, E. S. and Kruglyak, L. (1995). Genetic dissection of complex traits: guidelines for interpreting and reporting linkage results, Nature Genetics, 11 (3), 241–247.
[21] Manichaikul, A., Dupuis, J., Sen, S. and Broman, K. W. (2006). Poor performanceof bootstrap conﬁdence intervals for the location of a quantitative trait locus, Genetics, 174 (1), 481–489.
[22] Piterbarg, V. I. (1996). Asymptotic Methods in the Theory of Gaussian Processes
and Fields, Translations of Mathematical Monographs, 148, AMS.
[23] Rebaı̈, A., Goﬃnet, B. and Mangin, B. (1994). Approximate thresholds of interval
mapping tests for QTL detection, Genetics, 138 (1), 235–240.
[24] Sen, S. and Churchill, G. A. (2001). A statistical framework for quantitative trait
mapping, Genetics, 159 (1), 371–387.
[25] Siegmund, D. (1985). Sequential Analysis, Springer.

29[26] Siegmund, D. O. (1992). Tail approximations for maxima of random ﬁelds, in Probability Theory, L. H. Y. Chen, K. P. Choi, K. Hu and J-H. Lou (eds.), Walter deGruyter, 147–158.
[27] 鵜飼保雄 (2000). 「ゲノムレベルの遺伝解析 — MAP と QTL」, 東京大学出版会.
[28] Wu, C. F. J. (1983). On the convergence properties of the EM algorithm, Ann.
Statist., 11 (1), 95–103.
[29] Wu, R., Ma, C.-X., and Casella, G. (2007). Statistical Genetics of Quantitative
Traits, Springer.
[30] Yoshida, N. (2006). Polynomial type large deviation inequalities and convergence ofstatistical random ﬁelds, ISM Research Memorandum No. 1021, to appear in Ann.
Inst. Statist. Math.

30