Rで連鎖解析の基礎をなぞる - ryamadaの遺伝学・遺伝統計学メモ

家系図を個人とその親子関係のグラフとする
辺の数は、伝達の数。これをとする
- 伝達に関して、それが、父由来か母由来かを表すとすると、長さnのベクトルで、値は0か1かで表せる
  - $n=4$ の場合はこんな感じ( $2^n \times n$ の行列

TreeVector<-function(n=4){
 N<-2^n
 mat<-matrix(0,N,n)
 for(i in 2:N){
  mat[i,]<-mat[i-1,]
  mat[i,1]<-mat[i,1]+1
  for(j in 1:(n-1)){
   if(mat[i,j]==2){
    mat[i,j]<-0
    mat[i,j+1]<-mat[i,j+1]+1
   }else{
    j<-n
   }
  }
 }
 mat
}
n<-4 #伝達数
TreeVector(n=n)
      [,1] [,2] [,3] [,4]
 [1,]    0    0    0    0
 [2,]    1    0    0    0
 [3,]    0    1    0    0
 [4,]    1    1    0    0
 [5,]    0    0    1    0
 [6,]    1    0    1    0
 [7,]    0    1    1    0
 [8,]    1    1    1    0
 [9,]    0    0    0    1
[10,]    1    0    0    1
[11,]    0    1    0    1
[12,]    1    1    0    1
[13,]    0    0    1    1
[14,]    1    0    1    1
[15,]    0    1    1    1
[16,]    1    1    1    1

DNA配列上のすべての座位は、この $2^n$ のどれかのパターンである。それを $P_i;i=1,2,...,2^n$ とする
ある座位のパターンがPiで、別の座位のパターンがPjとすると、PiとPjのマンハッタン距離がその２座位間で、組み換えが起きた伝達の数である
- $2^n \times 2^n$ の組み合わせについて、それを行列化する

RecNumberMat<-as.matrix(dist(m,method="manhattan",diag=TRUE,upper=TRUE))
RecNumberMat
   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1  0 1 1 2 1 2 2 3 1  2  2  3  2  3  3  4
2  1 0 2 1 2 1 3 2 2  1  3  2  3  2  4  3
3  1 2 0 1 2 3 1 2 2  3  1  2  3  4  2  3
4  2 1 1 0 3 2 2 1 3  2  2  1  4  3  3  2
5  1 2 2 3 0 1 1 2 2  3  3  4  1  2  2  3
6  2 1 3 2 1 0 2 1 3  2  4  3  2  1  3  2
7  2 3 1 2 1 2 0 1 3  4  2  3  2  3  1  2
8  3 2 2 1 2 1 1 0 4  3  3  2  3  2  2  1
9  1 2 2 3 2 3 3 4 0  1  1  2  1  2  2  3
10 2 1 3 2 3 2 4 3 1  0  2  1  2  1  3  2
11 2 3 1 2 3 4 2 3 1  2  0  1  2  3  1  2
12 3 2 2 1 4 3 3 2 2  1  1  0  3  2  2  1
13 2 3 3 4 1 2 2 3 1  2  2  3  0  1  1  2
14 3 2 4 3 2 1 3 2 2  1  3  2  1  0  2  1
15 3 4 2 3 2 3 1 2 2  3  1  2  1  2  0  1
16 4 3 3 2 3 2 2 1 3  2  2  1  2  1  1  0

- 同様に、組み換えの起きなかった伝達の数は次のように表せる

NonRecNumberMat<-n-RecNumberMat #隣接ツリー間の非組み換え回数

座位間の組み換え割合。座位間は組み替えるか組み換えないかで、に分けられる
- 異なる染色体上の２座位であるとか、非常に遠い場合には、 $\theta=(1-\theta)=0.5$ 、非常に近くて、組み換えが起きないときに0となる
- ２座位間で、 $n$ 個の伝達があって、そのうち、 $k$ 個で組み換えがあり、 $n-k$ 個で組み換えなかったとすると、パターンの組み合わせには、尤度の違いができる
- それは、組み換え伝達数と組み換えなし伝達数の行列と、 $\theta$ を使って、次のように計算できる

theta<-0.3
theta^RecNumberMat*(1-theta)^NonRecNumberMat
...16x16行列なので、表示は省略...

２座位間に原因座位があるかどうかを探索する
- ２つのマーカーM1,M2がある
- それぞれのジェノタイプデータから、それぞれ、伝達パターン個の取り方をm1,m2という長さのベクトルで表す
  - ジェノタイプデータを与えて、このパターンを決めることになるが、今は、たった４伝達(夫婦と子２人とか)で考えているので、以降の計算がつまらなくなるので、以下のように適当にベクトルを作る。
- 表現型のデータがある
  - 表現型についても、優性・劣性の形式を仮定し、浸透率を仮定すると、 $2^n$ のどれをどれくらいの比率でもつべきかがわかる。Gとする。
  - Ｇがm1に近い場合をあえて作りたいので、以下のように勝手に作る。本当は、表現型のパターンと遺伝形式と浸透率(など)で地道に計算するべきものである

m1<-sample(c(0,1),n^2,0.2)
m1<-m1/sum(m1)
m2<-sample(c(0,1),n^2,0.2)
m2<-m2/sum(m2)
G<-0.9*m1+0.1*m2
G<-G/sum(G)

> m1
 [1] 0.0000000 0.1428571 0.1428571 0.0000000 0.0000000 0.0000000 0.0000000
 [8] 0.0000000 0.0000000 0.1428571 0.0000000 0.1428571 0.1428571 0.1428571
[15] 0.1428571 0.0000000
> m2
 [1] 0.1111111 0.1111111 0.1111111 0.0000000 0.1111111 0.0000000 0.0000000
 [8] 0.1111111 0.0000000 0.0000000 0.0000000 0.1111111 0.1111111 0.0000000
[15] 0.1111111 0.1111111
> G
 [1] 0.01111111 0.13968254 0.13968254 0.00000000 0.01111111 0.00000000
 [7] 0.00000000 0.01111111 0.00000000 0.12857143 0.00000000 0.13968254
[13] 0.13968254 0.12857143 0.13968254 0.01111111

M1とM2の間は組み換え割合が $\theta=0.3$ であるとわかっている(マッピングに用いるマーカーはそれがわかっているからマーカーとして用いることが便利)

theta<-0.4

M1とM2の間を細かく刻んで尤度を計算する位置をとる。それぞれの場所にGを置いたときに、は
- $\theta=\theta_{1G}(1-\theta_{2G})+\theta_{2G}(1-\theta_{1G})$ の関係を満たすが、そのようなときに、次のようにパラメタ化できる(昨日の記事を参照)

x<-seq(from=0,to=1,by=0.01)
t<-x*pi/2

theta1G<--sqrt((0.5-theta)/sin(2*t))*cos(t)+0.5
theta2G<--sqrt((0.5-theta)/sin(2*t))*sin(t)+0.5
range<-which(theta1G>=0 & theta2G>=0)
x<-x[range]
theta1G<-theta1G[range]
theta2G<-theta2G[range]

のそれぞれについて、m1,G,m2の伝達のパターンをすべて網羅して、そのときの尤度が計算したい
- $(2^n)^3$ の組み合わせパターンがある
- それを全部やるのも１つの方法

Like<-rep(0,length(x))
LikeNull<-rep(0,length(x))
for(i in 1:length(x)){ #Gの位置をずらすループ
 for(j in 1:length(m1)){ #m1のパターンのループ
  for(k in 1:length(G)){ #Gのパターンのループ
   for(l in 1:length(m2)){ #m2のパターンのループ
    tmp1<-m1[j]*G[k]*m2[l] #組み換えによる起こりやすさを考慮しない確率
    numrec1G<-RecNumberMat[j,k] #M1とGとの間の組み換えあり伝達数
    numNonrec1G<-NonRecNumberMat[j,k] #M1とGとの間の組み換えなし伝達数
    numrecG2<-RecNumberMat[k,l] #M2とGとの間の組み換えあり伝達数
    numNonrecG2<-NonRecNumberMat[k,l] #M2とGとの間の組み換えなし伝達数
    tmp2<-tmp1*theta1G[i]^numrec1G*(1-theta1G[i])^numNonrec1G*theta2G[i]^numrecG2*(1-theta2G[i])^numNonrecG2
    Like[i]<-tmp2+Like[i] #すべての伝達パターン組み合わせについて足し合わせ。Gの位置ごとに集計
    numrec12<-RecNumberMat[j,l] #帰無仮説では、GはよそにあるのでM1とM2の間の組み換えあり伝達数を使う
    numNonrec12<-NonRecNumberMat[j,l] #同様にM1とM2の間の組み換え無し伝達数
    tmp3<-tmp1*theta^numrec12*(1-theta)^numNonrec12*(1/2)^n #M1とM2の間は組み換え有り無しを考慮し、GについてはtM1-M2のつながりパターンに対して、全伝達で組み換え割合が1/2なので、(1/2)^nをかけている
    LikeNull[i]<-LikeNull[i]+tmp3
   }
  }
 }
}
ylim<-c(0,max(Like,LikeNull))
plot(x,Like,ylim=ylim,type="l")
par(new=T)
plot(x,LikeNull,ylim=ylim,type="l",col="red")

尤度の計算をするのに、ループを４重にしてぐるぐるするのも気が利かないので、方法を変える
- M1からスタートする。このとき $2^n$ パターンの比率が決まる
- 次にGと組み合わせる

m1G<-m1%*%t(G) #16x16行列

- $(2^n)^2$ のパターンについて、考慮して、Gの位置に応じて、組み換えありなしを考慮する
- $(2^n)^2$ のパターンごとに尤度が出る

 m1Gx<-m1G*theta1G[i]^RecNumberMat*(1-theta1G[i])^NonRecNumberMat

- 次にM2との関係を考える
- ただし、M2とこの $(2^n)^2$ の関係は、組み換えについてだけ、考慮するのであれば、この時点で、M2に面しているもの(それはGの $2^n$ パターン)のみによって決まるので、それに関して尤度を集計する

m1Gx<-apply(m1Gx,2,sum) #長さ16のベクトルに戻る

- GとM2の間の組み換えありなし伝達により尤度を変化させる

 m2G<-m1Gx%*%t(m2)
 m2G<-m2G*theta2G[i]^RecNumberMat*(1-theta2G[i])^NonRecNumberMat

- 全パターンについて足し合わせる

Like2[i]<-sum(m2G)

２つの方法の尤度を帰無仮説のそれと較べてみる

plot(x,Like,type="l",ylim=ylim)
par(new=T)
plot(x,LikeNull,type="l",col="red",ylim=ylim)
par(new=T)
plot(x,Like2,type="h",col="blue",ylim=ylim)

尤度比検定をすると

LOD<-log(Like/LikeNull,10)

対立仮説と帰無仮説の尤度を計算する

CalcLike<-function(n,m1,m2,G,theta1,theta2){
 trvec<-TreeVector(n=n)
 RecNumberMat<-as.matrix(dist(trvec,method="manhattan",diag=TRUE,upper=TRUE)) #隣接ツリー間の組み換え回数
 NonRecNumberMat<-n-RecNumberMat #隣接ツリー間の非組み換え回数
 LL<-0
 for(j in 1:length(m1)){
  for(k in 1:length(G)){
   for(l in 1:length(m2)){
    tmp1<-m1[j]*G[k]*m2[l]
    numrec1G<-RecNumberMat[j,k]
    numNonrec1G<-NonRecNumberMat[j,k]
    numrecG2<-RecNumberMat[k,l]
    numNonrecG2<-NonRecNumberMat[k,l]
    tmp2<-tmp1*theta1^numrec1G*(1-theta1)^numNonrec1G*theta2^numrecG2*(1-theta2)^numNonrecG2
    LL<-tmp2+LL
   }
  }
 }
 LL
}

CalcLike2<-function(n,m1,m2,G,theta1,theta2){
 m1G<-m1%*%t(G)
 m1Gx<-m1G*theta1^RecNumberMat*(1-theta1)^NonRecNumberMat
 m1Gx<-apply(m1Gx,2,sum)
 m2G<-m1Gx%*%t(m2)
 m2G<-m2G*theta2^RecNumberMat*(1-theta2)^NonRecNumberMat
 sum(m2G)
}

CalcLikeNull<-function(n,m1,m2,G,theta){
 trvec<-TreeVector(n=n)
 RecNumberMat<-as.matrix(dist(trvec,method="manhattan",diag=TRUE,upper=TRUE)) #隣接ツリー間の組み換え回数
 NonRecNumberMat<-n-RecNumberMat #隣接ツリー間の非組み換え回数
 LL<-0
 for(j in 1:length(m1)){
  for(k in 1:length(G)){
   for(l in 1:length(m2)){
    tmp1<-m1[j]*G[k]*m2[l]
    numrec1G<-RecNumberMat[j,l]
    numNonrec1G<-NonRecNumberMat[j,l]
    tmp2<-tmp1*theta^numrec1G*(1-theta)^numNonrec1G*(1/2)^n
    LL<-tmp2+LL
   }
  }
 }
 LL


}