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

いくつかの特徴を反映させたモデル

R 関節リウマチ
  • (一昨日)まで、疾患の病像のモデル化を考えている
  • その本当の発端は、「遺伝疫学(病気の遺伝因子の学問)」がテーマだからであるが、
  • もう少し、小さな枠組みでの発端は、以下の諸点になる
    • GEEとそれが使う一般化線形モデルから、「少々の複雑なパラメタ構成は解析の対象になるだろう」、という目算と、
    • 「疾患は、たくさんの因子が集まって決める多くの状態を分類したものとみなしたい」、という気持ちと、
    • 「その分類にあたっては、多くの場合、１次元情報に射影している」、という(おそらく正しい)類推と、
    • 「病気を特徴づける症状・症候は、数え上げることができる病変の数として現れる(炎症関節数、荒廃糸球体数、脳内変性性プラーク数・・・)」、という事実
  • 反映させたい条件
    • もっともらしい有病率(リウマチなら１％弱)
    • もっともらしい発病率の年齢変化(リウマチなら、高齢で好発)
    • もっともらしい病変数分布(リウマチなら、炎症関節数０(リウマチではない）がもっとも多いが、リウマチ患者で炎症関節数の分布をとると、数個(複数)の関節が炎症している例が多く、１個は少ない(分類基準では否定的)、少数例では、非常に多くの関節が炎症を起こしている
    • 病変箇所は変化しうる(ある病変が治っても、別の病変が新規に現れることで、病変箇所数は０にならないかもしれないが、病変は移動しうる)
    • 一度病変となったら、治らない(不可逆変化)もありえるものとする
    • 集団としては、有病率が年齢の関数になっているとしても、個人に着目すると、さまざまな経過をとりうる(一度きりのエピソード、有症状状態が遷延するエピソード、症状が消えたり出たりする、長引いたからと言って、必ずしも治らないというわけでもない)
  • どうしてこんなことを「このメモ」で扱いたいか、と言えば、「適切なモデルがあれば、そのモデルの中のどんなパラメタに多型がアレル特異的な機能差をもたらしているのかを仮説にする」ことができ、さらには、「その仮説に照らして、ジェノタイプ-フェノタイプ データの当てはまりの良さを計量する」ことができるから。
  • 最後に。モデルはできるだけ単純なのが、もちろんよい。
  • この話しのうち、病的状態になるのが、T回の同確率事象の積み重ねであると考えるとすると、それはポアソン過程で加わる障害がT回積み重なるまでの時間の分布であって、それはガンマ分布であらわされる(こちらを参照)
    • ガンマ分布を描くと
plot(function(x)dgamma(x,shape=10,rate=0.1),from=0,to=300)
  • ポアッソン過程を"rpois"でシミュレーションしてT回積み重ね時間の分布を描いて(ガンマ分布様であることを確かめるには
RP<-rpois(1000000,1)
RP2<-cumsum(RP)

N<-10
res<-c(0)
tmpi<-0
for(i in 1:length(RP2)){
	if(RP2[i]>N){
		tmp<-i-tmpi
		res<-c(res,tmp)
		tmpi<-i
		RP2<-RP2-N
	}
}
res<-res[2:length(res)]
length(res)
hist(res)
# 疾病状態は１本道状の木とする
# その状態数はNsとする
# 離散的経時変化を考えNt時点を考える
# 行列Sは状態推移を表す
# 状態は１本道の隣接状態にのみ移れるものとする
# S[i,j]は状態jから状態iに移る確率
# 単純化したモデルとして次のようなモデルとする
# 第一段階と最終段階以外のすべての段階においてS[,j]は同じ
# r S[j,j](そのままの確率)を指定
# s S[j-1,j]とPs[j+1,j]との比率を指定
# s=1:完全修復 s=0:非可逆 s=0.5:酔歩50:50
#################
# 病変はNdを最大数とする離散的な箇所数で表し、それぞれの病変は独立とする
# 病変は、疾病状態が定める関数によって新規に発生する
# 病変はつねに、修復される(されない場合は不可逆病変)
# 新規病変発生率関数
# たとえば
# fs<-rep(0,Ns)
# fs[(T+1):Ns]<-p
# では、j>Tでは一定値pj<=Tでは0
# 別の関数
# f(s)=min(1,(s-T)p)
# では、sはステージの段階番号、Tは閾値。
# Tより大のステージでは、線形に発生率が上昇するモデル
# f(s)をあらかじめ指定しているやり方(StageMultiLegion())と
# f(s)を与えるやりかた(StageMultiLegion2())とを作成してある
# 修復は確率g(s)でなされるものとする
# sによらずに一定とするモデルとして、g(s)=qとすることにする
# q=0の場合は不可逆病変であることを意味する
# 病変数がkであるような状態をB(k)とする　
# 病変数の変化は２段階で考えることとする
# 第１段階は、新規の病変部が発生するのみで、修復されない段階(b1 step)
# 第２段階は、新規の病変発生はなく、修復されるのみの段階(b2 step)
# b1 stepでは、病変数iの状態から病変数jの状態へは、
# 以下の確率で推移することとする
# when j>i f(i->j)=choose(Nd-i,Nd-(j-i)) f(s)^(Nd-(j-i)) (1-f(s))^j
# when j=i 1-sum(pi->j(j>i))
# when j<i 0
# b2 stepでは、
# when j>i 0
# when j=i 1-sum(pi->j(j<i))
# when j<i p(i->j)=choose(i,(i-j)) g(s)^(i-j) (1-g(s))^j
############################

Ns<-10 # No. stages
Nt<-200 # No. discrete time points
Nd<-20
T<-6
p<-0.2 
q<-0.1 # 病変の修復力
r<-0.5
s<-0.0  # 病気の状態の修復力 s=1 完全修復 s=0 非可逆 s=0.5 酔歩50:50


StageMultiLegion<-function(Ns=Ns,Nt=Nt,Nd=Nd,T=T,p=p,q=q,r=r,s=s){



fs<-((1:Ns)-T)*p
fs[which(fs>1)]<-1
fs[which(fs<0)]<-0
gs<-rep(q,Ns)

S<-matrix(0,Ns,Ns)
diag(S)<-r



for(i in 1:(Ns-1)){
	S[i,i+1]<-(1-r)*s
	S[i+1,i]<-(1-r)*(1-s)
}


S[1,1]<-1-S[2,1]
S[Ns,Ns]<-1-S[Ns-1,Ns]
apply(S,2,sum)

B1<-B2<-array(0,c(Nd+1,Nd+1,Ns))

for(i1 in 1:Ns){
	for(i2 in 1:(Nd)){
		for(i3 in (i2+1):(Nd+1)){
			# B1にとって：i3-1:後の病変数、i2-1:前の病変数
			# B2にとって：i3-1:前の病変数、i2-1:後の病変数
			B1[i3,i2,i1]<-choose(Nd-(i2-1),i3-i2)*fs[i1]^(i3-i2)*(1-fs[i1])^(Nd-(i3-1))
			B2[i2,i3,i1]<-choose((i3-1),i3-i2)*gs[i1]^(i3-i2)*(1-gs[i1])^(i2-1)
		}
	}
	tmpsum1<-apply(B1[,,i1],2,sum)
	tmpsum2<-apply(B2[,,i1],2,sum)
	diag(B1[,,i1])<-1-tmpsum1
	diag(B2[,,i1])<-1-tmpsum2
}

serialMultMat<-function(B,X){
	D<-dim(X)
	R<-matrix(0,D[1],D[2])
	for(i in 1:D[2]){
		R[,i]<-B[,,i]%*%X[,i]
	}
	R
}

#ResS<-matrix(0,Ns,Nt+1)
ResB<-array(0,c(Ns,Nd+1,Nt+1))


#ResS[1,1]<-1
ResB[1,1,1]<-1

for(i in 2:(Nt+1)){
	ResB[,,i]<-S%*%t(serialMultMat(B2,serialMultMat(B1,t(ResB[,,i-1]))))
	
}

# ステージ別比率の経時変化
# 病変数別比率の経時変化

ResStage<-matrix(0,Ns,Nt+1)
ResNd<-matrix(0,Nd+1,Nt+1)

for(i in 1:(Nt+1)){
	ResStage[,i]<-apply(ResB[,,i],1,sum)
	ResNd[,i]<-apply(ResB[,,i],2,sum)
}
par(mfcol=c(1,2))
matplot(t(ResStage),type="l")
matplot(t(ResNd),type="l")
par(mfcol=c(1,1))

list(shortsummary=list(Ns=Ns,Nt=Nt,Nd=Nd,T=T,p=p,q=q,r=r,s=s,S=S,B1=B1,B2=B2,fs=fs,gs=gs),ResB=ResB,ResStage=ResStage,ResNd=ResNd)

}
printResStageResNd<-function(x){
par(mfcol=c(1,2))
matplot(t(x$ResStage),type="l")
matplot(t(x$ResNd),type="l")
par(mfcol=c(1,1))
}
printShortSummary<-function(x){
	print(x$shortsummary)
}


StageMultiLegion2<-function(Ns=Ns,Nt=Nt,Nd=Nd,fs=fs,gs=gs,r=r,s=s){


S<-matrix(0,Ns,Ns)
diag(S)<-r



for(i in 1:(Ns-1)){
	S[i,i+1]<-(1-r)*s
	S[i+1,i]<-(1-r)*(1-s)
}


S[1,1]<-1-S[2,1]
S[Ns,Ns]<-1-S[Ns-1,Ns]
apply(S,2,sum)

B1<-B2<-array(0,c(Nd+1,Nd+1,Ns))

for(i1 in 1:Ns){
	for(i2 in 1:(Nd)){
		for(i3 in (i2+1):(Nd+1)){
			# B1にとって：i3-1:後の病変数、i2-1:前の病変数
			# B2にとって：i3-1:前の病変数、i2-1:後の病変数
			B1[i3,i2,i1]<-choose(Nd-(i2-1),i3-i2)*fs[i1]^(i3-i2)*(1-fs[i1])^(Nd-(i3-1))
			B2[i2,i3,i1]<-choose((i3-1),i3-i2)*gs[i1]^(i3-i2)*(1-gs[i1])^(i2-1)
		}
	}
	tmpsum1<-apply(B1[,,i1],2,sum)
	tmpsum2<-apply(B2[,,i1],2,sum)
	diag(B1[,,i1])<-1-tmpsum1
	diag(B2[,,i1])<-1-tmpsum2
}

serialMultMat<-function(B,X){
	D<-dim(X)
	R<-matrix(0,D[1],D[2])
	for(i in 1:D[2]){
		R[,i]<-B[,,i]%*%X[,i]
	}
	R
}

#ResS<-matrix(0,Ns,Nt+1)
ResB<-array(0,c(Ns,Nd+1,Nt+1))


#ResS[1,1]<-1
ResB[1,1,1]<-1

for(i in 2:(Nt+1)){
	ResB[,,i]<-S%*%t(serialMultMat(B2,serialMultMat(B1,t(ResB[,,i-1]))))
	
}

# ステージ別比率の経時変化
# 病変数別比率の経時変化

ResStage<-matrix(0,Ns,Nt+1)
ResNd<-matrix(0,Nd+1,Nt+1)

for(i in 1:(Nt+1)){
	ResStage[,i]<-apply(ResB[,,i],1,sum)
	ResNd[,i]<-apply(ResB[,,i],2,sum)
}
par(mfcol=c(1,2))
matplot(t(ResStage),type="l")
matplot(t(ResNd),type="l")
par(mfcol=c(1,1))

list(shortsummary=list(Ns=Ns,Nt=Nt,Nd=Nd,T=T,p=p,q=q,r=r,s=s,S=S,B1=B1,B2=B2,fs=fs,gs=gs),ResB=ResB,ResStage=ResStage,ResNd=ResNd)

}
printResStageResNd<-function(x){
par(mfcol=c(1,2))
matplot(t(x$ResStage),type="l")
matplot(t(x$ResNd),type="l")
par(mfcol=c(1,1))
}
printShortSummary<-function(x){
	print(x$shortsummary)
}


#################
# 実行例
#################
Ns<-10 # No. stages
Nt<-100 # No. discrete time points
Nd<-20
T<-8
p<-0.03 
q<-0.02
r<-0.9
s<-0.5  # 修復力 s=1 完全修復 s=0 非可逆 s=0.5 酔歩50:50

outSML<-StageMultiLegion(Ns=Ns,Nt=Nt,Nd=Nd,T=T,p=p,q=q,r=r,s=s)
printShortSummary(outSML)
printResStageResNd(outSML)
matplot(t(outSML$ResNd[2:(Nd+1),]),type="l")


plot(outSML$ResNd[2:(Nd+1),Nt+1])

#########################

Ns<-10 # No. stages
Nt<-100 # No. discrete time points
Nd<-20

T<-8
p<-0.4 
q<-0.7
r<-0.85
s<-0.5  # 修復力 s=1 完全修復 s=0 非可逆 s=0.5 酔歩50:50

fs<-rep(0,Ns)
fs[(T+1):Ns]<-p

gs<-rep(q,Ns)

outSML2<-StageMultiLegion2(Ns=Ns,Nt=Nt,Nd=Nd,fs=fs,gs=gs,r=r,s=s)
printShortSummary(outSML2)
printResStageResNd(outSML2)
matplot(t(outSML2$ResNd[2:(Nd+1),]),type="l")


plot(outSML2$ResNd[2:(Nd+1),Nt+1])

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


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