すべてのハプロタイプは世界に1本だけ配列のとき

Nh<-10 # 染色体本数
Nm<-20 # SNP数
# 適当にハプロタイプを作る
H<-matrix(sample(c(0,1),Nh*Nm,replace=TRUE),Nh,Nm)
  • 染色体が2本合わさってディプロタイプを作る
  • 染色体のプールからランダムにペアを作るとそれはHWE仮定
ShuffleAndDiploid<-function(h){
	Nh<-length(h[,1])
	H<-h[sample(1:Nh),]

	H1<-H[1:(Nh/2),]
	H2<-H[(Nh/2+1):Nh,]
	H1+H2
}
  • 多型性のないマーカーをシミュレーションで作ったら、除去してやりたい
RemoveNonPolymorphic<-function(h){
	m<-apply(h,2,min)
	M<-apply(h,2,max)
	h[,which(m!=M)]
}
  • 系統的に「癖」をつける
# ちょっと工夫して遠近関係に癖をつけよう
# 確率情報は非負
prob<-cos((1:Nm)*0.05)
prob<-prob-min(prob)
t<-sort(sample(1:Nm,Nh,prob=prob,replace=TRUE))

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/plott.jpeg

  • 系統的に作った染色体の並び方をでたらめにする
H<-H[sample(1:Nh),]
  • LDの様子を図にしよう
image(cor(H))

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/LD.jpeg

  • ハプロタイプ同士の遠近関係(距離)はマンハッタン距離を用いる、そしてNJ法で木表示しよう

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/htree.jpeg

# ハプロタイプが1本の線になっている場合
Nh<-200
Nm<-100
H<-matrix(0,Nh,Nm)
# ちょっと工夫して遠近感系に癖をつけよう
# 確率情報は非負
prob<-cos((1:Nm)*0.05)
prob<-prob-min(prob)
t<-sort(sample(1:Nm,Nh,prob=prob,replace=TRUE))
# 癖の具合を図で示す
plot(t)
# 癖に応じてハプロタイプを定める
for(i in 1:Nh){
	H[i,1:t[i]]<-1
}
# 非多型性マーカーを除去
H<-RemoveNonPolymorphic(H)
# マーカーが減っていることの確認
dim(H)
# LDの具合を図にしよう
image(cor(H))
# 木にしよう
library(ape)
plot(nj(dist(H,method="manhattan")),type="u")
# Diplotypeを作ろう
G<-ShuffleAndDiploid(H)
plot(nj(dist(G,method="manhattan")),type="u")
  • Diplotypeの木

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/gtree.jpeg

  • おまけ
genes<-matrix(0,1,2)
genes[1,1]<-1
genes[1,2]<-length(G[1,])
p<-c(rep(0,Nh/4),rep(1,Nh/4))
Ks<-3
CheckDists<-0:6
sss3<-SpherePowerGene(p=p,g=G,ptype=c(0,1),genes=genes,Ks=Ks,CheckDists=CheckDists,nperm=100)

matplot(t(sss3$powers),type="l")