数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

Rのarrayでハプロタイプのディプロタイプを表す

R
  • こちらを見た
    • 自分と異なる発想に触れると考えが広がるのでありがたい
  • ２座位の場合
    • 座位１のアレル数n_1,座位２のアレル数n_2
    • ハプロタイプ種類数n_h=n_1\times n_2
    • ディプロタイプ種類数(父由来(paternal),母由来(maternal)を区別する)n_h^2
  • vector,matrix,array。「平衡」状態は掛け算でできて、しかも、要素ごとに行わなくてもできる(matrixの掛け算がarrayにも拡張できているところがRらしくてよい)
N1<-3    #ローカス1のアレル数
N2<-4    #ローカス2のアレル数
# アレル頻度
f1<-runif(N1)
f1<-f1/sum(f1)
f2<-runif(N2)
f2<-f2/sum(f2)

# 連鎖平衡のときのハプロタイプ頻度
H<-f1%*%t(f2)

# 確認
apply(H,1,sum)
apply(H,2,sum)
f1
f2

# HWEのときのディプロタイプ頻度
# ベクトル化して、掛け合わせて、アレイの枠に納める
B2<-array(c(H)%*%t(c(H)),c(N1,N2,N1,N2))
# 要素ごとに計算するのと同じこと
B<-array(0,c(N1,N2,N1,N2))     #二倍体の存在頻度B:これはarray   (N1×N2×N1×N2の4次元)
for (i in 1:N1){
	for(j in 1:N2){
		for (k in 1:N1){
			for (l in 1:N2){
				B[i,j,k,l]<-H[i,j]*H[k,l]
				}
			}
		}
	}


#確認
B-B2
apply(B2,1,sum)
apply(B2,2,sum)
apply(B2,3,sum)
apply(B2,4,sum)
f1
f2
  • 実行歴
> N1<-3    #ローカス1のアレル数
> N2<-4    #ローカス2のアレル数
> # アレル頻度
> f1<-runif(N1)
> f1<-f1/sum(f1)
> f2<-runif(N2)
> f2<-f2/sum(f2)
> 
> # 連鎖平衡のときのハプロタイプ頻度
> H<-f1%*%t(f2)
> 
> # 確認
> apply(H,1,sum)
[1] 0.3662530 0.4191157 0.2146312
> apply(H,2,sum)
[1] 0.007674897 0.331545544 0.062147874 0.598631686
> f1
[1] 0.3662530 0.4191157 0.2146312
> f2
[1] 0.007674897 0.331545544 0.062147874 0.598631686
> 
> # HWEのときのディプロタイプ頻度
> 
> B<-array(0,c(N1,N2,N1,N2))     #二倍体の存在頻度B:これはarray   (N1×N2×N1×N2の4次元)
> for (i in 1:N1){
+ for(j in 1:N2){
+ for (k in 1:N1){
+ for (l in 1:N2){
+ B[i,j,k,l]<-H[i,j]*H[k,l]
+ }
+ }
+ }
+ }
> B2<-array(c(H)%*%t(c(H)),c(N1,N2,N1,N2))
> 
> #確認
> B-B2
, , 1, 1

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2, 1

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 3, 1

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 1, 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2, 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 3, 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 1, 3

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2, 3

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 3, 3

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 1, 4

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2, 4

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 3, 4

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

> apply(B2,1,sum)
[1] 0.3662530 0.4191157 0.2146312
> apply(B2,2,sum)
[1] 0.007674897 0.331545544 0.062147874 0.598631686
> apply(B2,3,sum)
[1] 0.3662530 0.4191157 0.2146312
> apply(B2,4,sum)
[1] 0.007674897 0.331545544 0.062147874 0.598631686
> f1
[1] 0.3662530 0.4191157 0.2146312
> f2
[1] 0.007674897 0.331545544 0.062147874 0.598631686
>
  • この先、データハンドリング上、方向は大きく分けて２つあるか(やらなくてもよいけれど…)
    • 座位数を任意にする
    • 「平衡(LE(連鎖平衡)、HWE)」からの逸脱をさせる
      • LEからの逸脱は交叉・組み換え
      • HWEからの逸脱はfixation indexなど
      • 『ryamada本』なら、こんな図のあたり

http://www.genome.med.kyoto-u.ac.jp/func-gen-photo/albums/StatGenetTextbook/1-16.jpeg

  • 座位数を任意にする
    • ２k座位の場合
      • 座位１のアレル数n_1,座位２のアレル数n_2,...,座位kのアレル数n_k
      • ハプロタイプ種類数n_h=n_1\times n_2\times ... \times n_k=\prod_{i=1}^k n_k
      • ディプロタイプ種類数(父由来(paternal),母由来(maternal)を区別する)n_h^2
        • この２乗は２倍体の２なので不変
# 座位数
Nm<-4
Na<-rpois(Nm,3)+1
Na
#最大座位数
maxNa<-max(Na)
# アレル頻度
F<-matrix(0,Nm,maxNa)
F
for(i in 1:Nm){
 F[i,1:Na[i]]<-runif(Na[i])
 F[i,]<-F[i,]/sum(F[i,])
}
F
# 確認
apply(F,1,sum)


# 連鎖平衡のときのハプロタイプ頻度
H<-array(c(F[1,]),rep(maxNa,1))
for(i in 2:Nm){
 H<-array(c(H)%*%t(F[i,]),rep(maxNa,i))
}
#H<-array(F[1,]%*%t(F[2,]),rep(maxNa,2))
#H<-array(c(H)%*%t(F[3,]),rep(maxNa,3))

# 確認
for(i in 1:Nm){
 print(apply(H,i,sum))
}
F

apply(H,1,sum)
apply(H,2,sum)
f1
f2

# HWEのときのディプロタイプ頻度


B2<-array(c(H)%*%t(c(H)),rep(maxNa,Nm*2))
#確認
sum(B2)-1
for(i in 1:(Nm*2)){
 print(apply(B2,i,sum))
}
F
  • 実行歴
> # 座位数
> Nm<-4
> Na<-rpois(Nm,3)+1
> Na
[1] 4 3 5 5
> #最大座位数
> maxNa<-max(Na)
> # アレル頻度
> F<-matrix(0,Nm,maxNa)
> F
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
> for(i in 1:Nm){
+  F[i,1:Na[i]]<-runif(Na[i])
+  F[i,]<-F[i,]/sum(F[i,])
+ }
> F
          [,1]      [,2]       [,3]       [,4]       [,5]
[1,] 0.4004534 0.2042133 0.32550353 0.06982986 0.00000000
[2,] 0.3457009 0.1557042 0.49859489 0.00000000 0.00000000
[3,] 0.3472646 0.1618986 0.04297426 0.18060354 0.26725901
[4,] 0.2898551 0.2479129 0.14146643 0.27727000 0.04349553
> # 確認
> apply(F,1,sum)
[1] 1 1 1 1
> 
> 
> # 連鎖平衡のときのハプロタイプ頻度
> H<-array(c(F[1,]),rep(maxNa,1))
> for(i in 2:Nm){
+  H<-array(c(H)%*%t(F[i,]),rep(maxNa,i))
+ }
> #H<-array(F[1,]%*%t(F[2,]),rep(maxNa,2))
> #H<-array(c(H)%*%t(F[3,]),rep(maxNa,3))
> 
> # 確認
> for(i in 1:Nm){
+  print(apply(H,i,sum))
+ }
[1] 0.40045337 0.20421325 0.32550353 0.06982986 0.00000000
[1] 0.3457009 0.1557042 0.4985949 0.0000000 0.0000000
[1] 0.34726459 0.16189860 0.04297426 0.18060354 0.26725901
[1] 0.28985514 0.24791290 0.14146643 0.27727000 0.04349553
> F
          [,1]      [,2]       [,3]       [,4]       [,5]
[1,] 0.4004534 0.2042133 0.32550353 0.06982986 0.00000000
[2,] 0.3457009 0.1557042 0.49859489 0.00000000 0.00000000
[3,] 0.3472646 0.1618986 0.04297426 0.18060354 0.26725901
[4,] 0.2898551 0.2479129 0.14146643 0.27727000 0.04349553
> 
> apply(H,1,sum)
[1] 0.40045337 0.20421325 0.32550353 0.06982986 0.00000000
> apply(H,2,sum)
[1] 0.3457009 0.1557042 0.4985949 0.0000000 0.0000000
> f1
[1] 0.3662530 0.4191157 0.2146312
> f2
[1] 0.007674897 0.331545544 0.062147874 0.598631686
> 
> # HWEのときのディプロタイプ頻度
> 
> 
> B2<-array(c(H)%*%t(c(H)),rep(maxNa,Nm*2))
> #確認
> sum(B2)-1
[1] 0
> for(i in 1:(Nm*2)){
+  print(apply(B2,i,sum))
+ }
[1] 0.40045337 0.20421325 0.32550353 0.06982986 0.00000000
[1] 0.3457009 0.1557042 0.4985949 0.0000000 0.0000000
[1] 0.34726459 0.16189860 0.04297426 0.18060354 0.26725901
[1] 0.28985514 0.24791290 0.14146643 0.27727000 0.04349553
[1] 0.40045337 0.20421325 0.32550353 0.06982986 0.00000000
[1] 0.3457009 0.1557042 0.4985949 0.0000000 0.0000000
[1] 0.34726459 0.16189860 0.04297426 0.18060354 0.26725901
[1] 0.28985514 0.24791290 0.14146643 0.27727000 0.04349553
> F
          [,1]      [,2]       [,3]       [,4]       [,5]
[1,] 0.4004534 0.2042133 0.32550353 0.06982986 0.00000000
[2,] 0.3457009 0.1557042 0.49859489 0.00000000 0.00000000
[3,] 0.3472646 0.1618986 0.04297426 0.18060354 0.26725901
[4,] 0.2898551 0.2479129 0.14146643 0.27727000 0.04349553
>

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------