多点遺伝的関連研究に対する主成分回帰アプローチ

背景

関連解析は家族ベースの連鎖解析より検出力が高い(Risch and Merikangas 1996)
密なSNP情報が利用可能であるため、連鎖不平衡構造に注目が集まっている。
- Clark 2004によるハプロタイプ解析のメリット
  - 特定のアミノ酸の組み合わせが、タンパク質折りたたみに影響している可能性がある
  - ハプロタイプは、祖先から受け継がれる染色体上の断片そのものである
  - ハプロタイプ解析は、単点解析よりも検出力が高い（いくつかのSNPを統合し、自由度を減らすため）
しかしハプロタイプ解析には問題がある
- マーカー数が増えると指数関数的に次元が増大し、多点検定の問題が悪化する
- また、まれなハプロタイプをどう扱うかという問題もある
もうひとつの方法は遺伝子型ベースの単点解析である
- ハプロタイプベースの方法より結果が理解しやすく、多く行われてきた
- 場合によっては単点解析のほうがハプロタイプ解析より検出力が高いという報告さえある
- この場合、有意差はpermutation法、Bonferroniの方法、Sidakの方法などで推定する
さらにもうひとつの方法は、候補SNPにおける遺伝子型スコアを独立変数とし、重回帰分析を行うという方法（Chapman et al. 2003)
- Fan and Knapp(2003)はHotellingのT²検定法を用いて関連解析を行った。→ロジスティック回帰分析と等価の方法
- 密に分布するSNPはたいてい相関しており、遺伝子型スコアの共線性が重大な問題である。
- 共線性のインパクトを減少させ、検出力を増大させるため対象領域の代表SNP (たとえばtagSNPs）を選択するという方法が、H-clustなどのアルゴリズムでは使用されている(Rinaldo et al, 2005)
これらの伝統的な方法のほか、多点遺伝子型データに特化した新しい方法も開発されてきている。
- ひとつの方法は、対象領域でのSNPの関連パターンをあてはめるというものである（意味不明）。
  - Lazzeroni(1998)は連鎖不平衡係数 $\delta$ を用いた方法。

$\delta=\frac{P(A|D)-P(A|N)}{1-P(A|N)}$

ここで、Dはdisease,Nはnormalの状態、AはマーカーアレルAをcarryする状態。

- - Cordell and Elston (1999)は、Fiellerの理論をもちいたもの。
  - Conti and Witte(2003)は、連鎖不平衡構造の階層化モデリング
  - Zhang et al.(2003)はベイズ適応回帰分析（？）を用いた方法。
- Schaid(2005)らは、個体ペア同士の遺伝子型の類似性をスコア化したprespecified kernelを利用したノンパラメトリックな方法を提案した。ペアのスコアの平均をケースとコントロールで比較した。
- Wessel and Schork(2006)は、kernelのかわりにゲノム類似性を測定して類似性ベースの回帰分析を提案した。
  - とても柔軟性に富んだ方法である
  - ゲノム類似性のたしかな測定法はわからない
- Wang and Elston(2007)では遺伝子型データのフーリエ変換を行い変換された成分の重み付けを行った。
今回著者らは、最低限の説明変数で最大限に必要なSNP遺伝子型多様性を捕らえて行う回帰モデルとして、主成分回帰 (principal component regression, PCReg)アプローチを採用した。
- この方法では、SNP遺伝子型スコアの共分散行列から計算された小数の主成分を説明変数として重回帰分析を行った。
- 主成分は、オリジナルのすべてのSNPの多様性のほとんどを説明できるように選択される。
- SNP間の連鎖不平衡が強ければ強いほど、必要な主成分は少数で済む。

方法

n individuals
遺伝子型データは0,1,2、対象領域にm個のSNP
g_ijは個体iのj番目SNPの遺伝子型
G=(g_ij)は $n \times m$ の遺伝子型行列
Gの各列の平均が0となるように中心化する
m SNPsの遺伝子型スコアの共分散行列を $\mathbf{V}$ とするが、これは $\mathbf{G}^t\mathbf{G} \div (n-1)$ に等しい。

(n-1)で割るのはなぜだろう？

$\lambda_1 \ge \lambda_2 \ge ... \ge \lambda_m$ は共分散行列 $\mathbf{V}$ のcharacteristic root（固有値？）で、 $\mathbf{a}_k$ は $\lambda_k$ と関連するcharacteristic vectorである（？？固有ベクトルのこと？）
定義から、 $\mathbf{Va}_k=\lambda_k\mathbf{a}_k$ である。
個体iのk番目のPC(主成分)は $(g_{i1},g_{i2},...,g_{im})\mathbf{a}_k$ である。
k番目のPCの分散は $\lambda_k$ である。
集合的に言うと、これらのPCは、SNP遺伝子型スコアの総分散をわけたものである。 $\sum^m_{k=1}\lambda_k$ は全SNPの遺伝子型スコアの分散と等しい（あるいは、 $\mathbf{V}$ の対角行列の和と等しい）。
PCをならびかえる方法から言って、特定のPCは、その他すべてのPCと比較して、オリジナルのSNP遺伝子型スコアの総分散を説明するに当たってより重要である。
詳しい説明は、どんな本にだってのってるよ

二つの例

PCのかたちは扱いにくいけど、簡単な例ならなんとかなる

一例目

すべての二組のSNPの相関係数は定数である。 $r \gt 0$
すべての遺伝子型スコアの分散は定数である。
- すべてのSNPのMAFが同じなら、この仮定は妥当であろう。
最大の固有値は $\lambda_1=s^2[1+(m-1)r]$ で、その固有ベクトルは $(m^{-\frac{1}{2}},...,m^{-\frac{1}{2}})^t$ である。
その他の固有値は全て $(1-r)s^2$
したがって最大の固有値による総分散の説明の割合はで、これは常により大きい。
- これはSNP数mが増えると減少するが、相関係数rが増えると増加する。r=1なら、mにかかわらず1である。
- m=30, r=0.5なら割合は0.52である。これは、たったひとつのPCで30SNPの遺伝子型データの52%の分散を説明できることを意味する。

二例目

SNPを二つの分画に区切り、それぞれ $m_1$ 、 $m_2$ の数とする。 $m_1+m_2=m$ である。
ブロック1の相関係数は定数 $r$ とする。
ブロック2のSNPは互いに相関しておらず、ブロック1のSNPとも相関しないものとする。
この場合、最大の固有値はで、固有ベクトルはで、このベクトルのうち非0要素が個となる。
- 最大のPCが説明する割合は $\frac{\{1+(m_1-1)r\}s^2}{ms^2}=\frac{rm_1}{m}+\frac{1-r}{m}$ である。
- これは最初の例よりも小さい値である。
つまり、SNP同士の相関が小さくなると、最初のPCで説明できる割合は減少する。
$m_2$ 個のPCが存在し、それぞれ分散 $s^2$ 、また $m_1-1$ 個のPCは分散 $(1-r)s^2$ である。
相関のあるブロック1のSNPのみが分解されたという結果である。
いろいろm、m₁、rを動かしてみると、最大のPCで説明できる割合はrに依存するがmにはあまり関係がないということがわかる。

PCSと回帰分析

最初の数個のPCでほとんどの多様性は説明できる。後述する実データでは、26SNPsを3つのPCで90%説明できる。
オリジナルのSNPに表現型との関連がないなら、PCの回帰係数は全て0になるはずである。検定はF検定を行う。
F統計量は自由度 $df_1$ 、 $df_2$ のF分布に従うが、 $df_1$ は回帰に用いたPCsの数であり、 $df_2$ は $n-(1+df_1)$ である。
$\mathbf{y}$ はn個体の表現型ベクトルとする。中心化し、平均を0とする。表現型は連続型になっているので回帰モデルを次のように書ける。

$\mathbf{y}=\mathbf{Gb}+\epsilon (1)$

yとGがともに中心化されているので切片は必要ない。
$\mathbf{A}=[\mathbf{a}_1,...,\mathbf{a}_m]$ は $n \times m$ 行列で、j列に $\mathbf{V}$ のj番目の固有ベクトルを配している。Aは直交系である。
$\mathbf{A}=[\mathbf{A}_1\vdots \mathbf{A}_2]$ と分割し、A₁はPCRegの説明変数として用いるPC、A₂は除外するPCである。
対応して $\mathbf{b}=[\mathbf{b}_1^t \vdots \mathbf{b}_2^t]^t$ と分割する。
回帰式(1)は次のように書ける。

$\mathbf{y}=\mathbf{GAA}^t\mathbf{b}+\epsilon=\mathbf{\tilde{G}}_1\mathbf{\tilde{b}}_1+\mathbf{\tilde{G}}_2\mathbf{\tilde{b}}_2+\epsilon$

ここで $\mathbf{\tilde{G}}_1=\mathbf{GA}_1$ 、 $\mathbf{\tilde{G}}_2=\mathbf{GA}_2$ 、 $\mathbf{\tilde{b}}_1=\mathbf{A}_1^t\mathbf{b}_1$ 、 $\mathbf{\tilde{b}}_2=\mathbf{A}_2^t\mathbf{b}_2$ である。
PCRegでは、回帰モデルをつぎのようにとる。

$\mathbf{y}=\mathbf{\tilde{G}}_1\mathbf{\tilde{b}}_1+\tilde(\epsilon)$

(1)式と同様、誤差 $\tilde{\epsilon}$ は平均0の正規分布に従う。
$\mathbf{\tilde{G}}_1=\mathbf{GA}_1$ であるから $\mathbf{\tilde{G}}_1$ の全項の和は0である。やはりこの式も切片をもたない。
ひきつづく解析でいくつのPCを使用するかにはたくさんのルールがある。
- ひとつの一般的なルールは、80-90%の総分散を説明できるPCを選択するというもの
- もうひとつは、固有値の平均を超える固有値を持つすべてのPCを選択するというもの。
  - scree plotが有用: $\{\lambda_1,...,\lambda_m\}$ を $\{1,...,m\}$ に対してプロットする。
- 統計学的検定を行うやりかたもある。
著者らの検討では「80-90%ルール」がうまくいった。
PCReg法は、これまでのいくつかの方法とは異なり、検出力はアレルのスコアリング法によらない。
より多くのPCをPCRegに含めれば、より検出力が高くなる、わけではない。
- PCRegのregression sum of squares (SSR)は

$SSR=\mathbf{y}^t\mathbf{\tilde{G}}_1(\mathbf{\tilde{G}}_1^t\mathbf{\tilde{G}}_1)^{-1}\mathbf{\tilde{G}}_1^t\mathbf{y}=\sum_k \lambda_k \tilde{b}^2_k$