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

データの保存

R
  • Rで何かしら、シミュレーションをしているとする
  • パラメタをいくつも振って実行しているとする
  • ２つの形で結果を残すことを考える
    • パッと見て結果について了解できる図を納めた１枚のファイル
    • 必要なデータの保管、全般
  • ファイル名の管理
    • パラメタで振った値をファイル名に反映させる
# ５種類のパラメタがあり、それぞれ値を１個以上持つ
Ls<-c(100)
Ns<-c(1000)
prevs<-c(0.01,0.05,0.1,0.2,0.5)
heritabilities<-c(0.01,0.05,0.1,0.2,0.3,0.4,0.5)
frisks<-c(0.1,0.25,0.5,1)

# ５個のパラメタの値の総当たりループ
for(L in Ls){
	for(N in Ns){
		for(prev in prevs){
			for(heritability in heritabilities){
				for(frisk in frisks){
# 処理。。。
# 処理。。。
# 処理。。。
# 処理。。。
 
# データ書き出し
# ファイル名の中にパラメタを込める
# 注意："="を使うと、".Rdata"ファイルとしてきちんとデータを持てないなど、制約があることに注意
# 画像ファイル
outfile<-paste("out","L=",L,"N_",N,"prev_",prev,"herit_",heritability,"frisk_",frisk,".pdf",sep="")
# ".Rdata"ファイル
outfile2<-paste("out","L_",L,"N_",N,"prev_",prev,"herit_",heritability,"frisk_",frisk,".Rdata",sep="")
# 画像ファイルへの書き出しのために…
pdf(file=outfile)
# １枚のpdfファイルに収めるために領域を分割
par(mfcol=c(3,3))
# プロットを繰り返す(３ｘ３＝９回まで)
par(mfcol=c(1,1))
dev.off()
# 以上で画像ファイルは終了

# データの保管
# オブジェクトをいくつか指定して納める場合(LSG,LSGdich,LSGcasecontの３つの場合が以下)
save(LSG,LSGdich,LSGcasecont,file = outfile2)

# いちいち面倒くさいので、全部残すとなったら、
save.image(file = outfile2)


				}
			}
		}
	}
}
  • おまけ(実際に延々と条件を変えて、出力させた例)
rm(list = ls())

RareVariantAccumulationTest<-function(P,G,rank=TRUE,Niter=1000){
	if(rank){
		R<-rank(P)
	}else{
		R<-P
	}
	T<-sum(R*G)
	Tperm<-rep(0,Niter)
	for(i in 1:Niter){
		Tperm[i]<-sum(sample(R)*G)
	}
	LSG<-c(length(which(Tperm<T)),length(which(Tperm==T)),length(which(Tperm>T)))/Niter
	p<-LSG[2]+LSG[3]
	return(list(p=p,LSG=LSG,Niter=Niter,rank=rank))
}


dichtomous<-TRUE
casecont<-TRUE
Ntrial<-100

Ls<-c(100)
Ns<-c(1000)
prevs<-c(0.01,0.05,0.1,0.2,0.5)
heritabilities<-c(0.01,0.05,0.1,0.2,0.3,0.4,0.5)
frisks<-c(0.1,0.25,0.5,1)


for(L in Ls){
	for(N in Ns){
		for(prev in prevs){
			for(heritability in heritabilities){
				for(frisk in frisks){
					


afcoef<-0.1
#afs<-abs(rnorm(L))*afcoef
afs<-rexp(L,1)
afs<-afs/max(afs)
afs[which(afs>=1)]<-min(afs) # afの上限は１
#plot(sort(afs))
afs<-sort(afs)
# RiskLocus閾値
LRcoef<-1
# LocusRisk
LocusRisk<-rep(0,L)



riskLoci<-which(afs<=LRcoef)
LocusRisk[riskLoci]<-sample(c(0,1),length(riskLoci),replace=TRUE,prob=c(1-frisk,frisk))
#LocusRisk<-LocusRisk*sample(c(-1,1),L,replace=TRUE)
LocusRisk<-LocusRisk*rnorm(L)
#LocusRisk<-LocusRisk/afs
#LocusRisk<-LocusRisk*rnorm(L)/afs
#LocusRisk<-LocusRisk*rnorm(L)/(afs*(1-afs))

#LocusRisk<-LocusRisk/(afs*(1-afs)) * rep(c(-1,1),L/2)


LSG<-LSGdich<-LSGcasecont<-matrix(0,Ntrial,3)
for(x in 1:Ntrial){


G<-matrix(0,N,L) # ジェノタイプ


for(i in 1:L){
	af<-afs[i]
	tmp<-sample(c(0,1,2),N,prob=c((1-af)^2,2*af*(1-af),af^2),replace=TRUE)
	G[,i]<-tmp
}
#G



IndRisk<-apply(t(G)*LocusRisk,2,sum)


#plot(IndRisk)

#plot(afs,apply(G,2,mean)/2)

#P<-rnorm(N) # フェノタイプ


P<-IndRisk+rnorm(N,sd=sqrt(var(IndRisk)*(1-heritability)/heritability))



RVATout<-RareVariantAccumulationTest(P=P,G=G,rank=TRUE,Niter=1000)
LSG[x,]<-RVATout$LSG

if(dichtomous){
	P01<-rep(0,N)
	P01[order(P,decreasing=TRUE)[1:(N*prev)]]<-1
	RVATout<-RareVariantAccumulationTest(P=P01,G=G,rank=TRUE,Niter=1000)
	LSGdich[x,]<-RVATout$LSG
	if(casecont){
		Nsample<-min(length(which(P01==0)),length(which(P01==1)))
		cases<-sample(which(P01==1),Nsample)
		conts<-sample(which(P01==0),Nsample)
		LSGcasecont[x,]<-RareVariantAccumulationTest(P=P01[c(cases,conts)],G=G[c(cases,conts),],rank=TRUE,Niter=1000)$LSG
	}

}

}

outfile<-paste("out","L=",L,"N_",N,"prev_",prev,"herit_",heritability,"frisk_",frisk,".pdf",sep="")
outfile2<-paste("out","L_",L,"N_",N,"prev_",prev,"herit_",heritability,"frisk_",frisk,".Rdata",sep="")
pdf(file=outfile)
par(mfcol=c(3,3))
plot(sort(afs))

plot(afs,LocusRisk)

plot(sort(IndRisk))
plot(sort(P))
plot(IndRisk,P)

matplot(LSG,type="b")
matplot(LSGdich,type="b")
matplot(LSGcasecont,type="b")
plot(LSGdich[,1],LSGcasecont[,1])
par(mfcol=c(1,1))
dev.off()

save(LSG,LSGdich,LSGcasecont,file = outfile2)

					
					
					
				}
			}
		}
	}
}