2011-09-29

■

サンプル数がある程度少ないとき、世界全体でのアレル頻度分布と、少ないサンプル数でのアレル頻度分布は変わるかもしれない。

実際のプログラミングについて考えると、最初100個RVを設定していても、多分100種類も観察されていなくて、重要なのはどちらかというと何個観察されているか、だと思う。

アレル頻度モデルは、下の図に近似させるが、サンプル200人なら実際どの程度それぞれの頻度のRVを持っているかを調べるのは重要かもしれない。
というのは、実際200人持ってきたとき、０．０１％のRVなんてほぼ見つからないからだ。
観察されるRVを100個、と決めるようにしないと、よくない、はず。

このアレル頻度モデルでやると、MSTの辺の和は常に長くなる感じがする。

この図のka=300くらいなので、黒そう。白いほど検出力がいい。

2011-09-28

■

http://www.hsph.harvard.edu/~xlin/SKAT-paper-and-supplement.pdf
のアレル頻度モデルのもととなった実際のデータ。RV93個。

アレル頻度かなり低い。これだけ低いと、RVをもってて1,2個だから、上手く結果がでないかもしれない。

大体ka=300くらい。kaとkrを変えてやってみたとき、kaが大きいとMSTの長さが長くなったような。

2011-09-26

■

レアバリアントが高リスクのとき、フェノタイプと関係性があるとでて、レアバリアントが低リスクのとき、フェノタイプと関係性がないという結果が出したい。
この高リスク、低リスクは、CVに比べて、というものなので、CVについてもう少し考えて組み込んで、リサンプリングの人数のところからやり直す。

多分CVに与えるリスクを上手いこと設定すると上手くいく感じがする。

今までのを無駄にしないなら、リスクを0でなく、RVに与えるもののうち最小のもの、という風にすると、モデルがあまり変わらず、いい結果が出そう。

2011-09-26

メモ

CVはリスク0で、アレル頻度20~50%。４個。
RVは高リスクで、アレル頻度0~10%。40個。

CVはリスク0だから、フェノタイプにまったく寄与せず、遺伝要因に占める割合は0。RVは、持っているとフェノタイプに大きく寄与する。

というモデルにおいて、MST検定を行うと、検出力が高い。→RVとCVのセットがフェノタイプに関係している。→しかし今CVは関係していない　誤検出？

でも、CVの部分は、実際のRVセットについて調べるときに勝手にリスク０で作れば、RVが関係しているかわかる。

ネガコンでCVを加えてみて結果がどうなうか見ればいい！

CV:RV=1:10のとき

ネガコン（リスクがアレル頻度にかかわらず±１）のときも検出力が高いが、RVのリスクが±１というのは、量的形質に関係しているから、検出できる方がいいはず。

リスク0のCVをいれると、RVのリスクが高い時、結果が出やすいモデルの形に近づく感じがする。

ネガコンを考え直そう。

RVのリスクが0に近い、フェノタイプに関係がない時がネガコン。

2011-09-26

■

CV4個で固定してRVの数を変えてみるとこうなった。

こうなる理由がつかめない。

モデルとしては、4個、と固定するより、RVの数によって変わる方が説明しやすそう。たとえば、

グラフの変化が小さい順にka=20,40,60,80
上のモデルのような確率分布に従う時、検出されやすい。というような。

今CV:RV=1:10でやっているけれど、CV1個でRV10個のときなどが上手くいかなさそう。。。やはり意味合いを考えねば。。。

2011-09-26

■

#バリアントの数が60のときの、CNを変える

Date<-"09261"

source("MST.R")

Np<-100						#MST何回か
HN<-50000

df<-200
cnnu<-c(20,30,40,50,60,70)
Rnu<-length(cnnu)
syu<-matrix(0,Np,Rnu)	#boxplot用
colnames(syu)<-cnnu

for (nnp in 1:Rnu){
CN<-cnnu[nnp]
RN<-60
ka<-40
#アレル頻度
af<-afb<-0
for (i in 1:df){
	afd<-(i-1)/df
	if (afd==0) afd<-0.001
	afb<-c(afb,rep(afd,round(RN*(pexp(i/df,ka)-pexp((i-1)/df,ka)))))
}
af<-afb[-1]

RN<-length(af)
risk<-rep(0,RN)
kr<-40
risk<-exp(-kr*af)*sign(rnorm(RN, mean=0, sd=1))

caf<-crisk<-rep(0,CN)
caf<-runif(CN,0.2,0.5)

af<-c(af,caf)
risk<-c(risk,crisk)

N<-length(af)
ir<-al<-matrix(0,HN,N)			#各ローカスにおけるリスク、持っているアレル
wei<-wei2<-P<-sr<-rep(0,HN)	#重み、フェノタイプ、リスクの総和
iaf<-rep(0,N)		#真のアレル頻度


eeee<-as.character(CN)
ssss<-paste(Date,"CN",eeee,".png")
png(ssss)          #デバイスドライバの用意
hist(af)
dev.off()

#個人のアレルの有無振り分け iさんについて
for(i in 1:HN){
	a1<-a2<-tmp1<-tmp2<-rep(0,N)
	a1<-runif(N)					#ＲＶの有無を決定
	a2<-runif(N)
	tmp1[which(a1<af)]<-1
	tmp2[which(a2<af)]<-1
	al[i,]<-tmp1+tmp2				#アレルをいくつもつか
	sr[i]<-sum(risk*al[i,])		#リスク総和
}

#for (i in 1:N) iaf[i]<-mean(al[,i])/2	#アレル頻度確かめ
#plot(af,iaf)                      #グラフを描画


#eeee<-as.character(df)
ssss<-paste(Date,"risk CN",eeee,".png")
png(ssss)          #デバイスドライバの用意
plot(density(sr))
dev.off()

#sdfrac<-0.98			#遺伝率
#P<-sr*sdfrac+rnorm(N,sd=sqrt(var(sr))*(1-sdfrac))
P<-sr    
#plot(sr,P)

#----------------ＭＳＴ・量的形質----------------
#par(ask=FALSE)

Ncaco<-200

for (counte in 1:Np){
	Cmst<-sample(1:HN,Ncaco,replace=FALSE)
	X<-cbind(P[c(Cmst)],al[c(Cmst),])
	
	MST<-StatsMST(X,perm=TRUE,Nperm=100,tobeshuffled=c(1),dist.method="manhattan")
	plot(sort(MST$PermL),ylim=range(c(MST$L,MST$PermL)))
	abline(h=MST$L,col=2)
	title(CN*1000+counte)
	ll<-MST$L
	cat1<-seq(from=0,to=1,by=0.01)
	cat2<-quantile(MST$PermL,cat1)
	cat3<-cat2-ll
	syu[counte,nnp]<-length(cat3[cat3<0])

}
save.image("09261vari.RData")
}