■

• 昨日の続き
• ケース・コントロールの２群が構造化のある集団からサンプリングされ、その構成にずれが生じているようなときに、そのサンプル２群のずれの方向にアレル頻度の勾配があるマーカーでのテストと、そのずれの方向とは直行する方向にアレル頻度の勾配があるマーカーとでは、統計量のインフレーションの具合が異なる
• 集団の方向とマーカーの頻度勾配の方向との一致がよければよいほど、統計量のずれは大きくなる
• それを図示するための、データシミュレーションRコマンドとその結果のプロットを３次元で見やすくするためのgnuplotをWindowsで使う話
• gnuplotのWindows版
  • 入手とインストールと起動
    • 入手先はこちらから、CVS version Windows binaries で。今回はそこの "MinGW"をダウンロードしました
    • 解凍して、"(解凍先)\gp43-winbin\gp43-winbin\gnuplot\binaries"にできるいくつかのexeファイルのうち、
      • gnuplot.exeはコマンドラインが立ち上がるもの
      • wgnuplot.exeがWindows用のGUIが立ち上がるもの
• 初期設定
  • 初めはフォントがつぶれています。GUIの一番上のgnuplotのロゴを左クリックするとOptionsが出るので、そこから"Choose Font"として、好きなフォントの種類とサイズにします。そのうえで、ロゴを左クリック、Optionsから、Updateすれば、次回から、同じフォントで起動します。 こちらにも記載がある通り、GUIの画面を右クリックしても同じことができます。
• ファイルにある３列データの3Dplot
  • GUIのツールバーの"ChDir"("change directory")からデータファイルのあるディレクトリへ移動します
  • そのファイルが、次のようなタブ区切りの３列データで、第１、第２、第３カラムのそれぞれが、x,y,z軸の値であるとします。ファイル名は、"data.txt"とします
  • "gnuplot>" なるプロンプトに　『splot "data.txt"』なるコマンドを打つと、掲載図（レコード数が下記の提示より多いです）が現れます。
  • ツールバーから、"Plot" → "3D plot" としてやると、コマンドである "splot "が表示されます。その後、" "にファイル名を囲んで指定してやればよいです
  • gnuplotがよいのは、ぐるぐる回して、適切なアングルが選びやすいこと。掲載図も２アングルからの絵を並べています
• 保存
  • データの保存は、"File"→"Save"でOK。
  • "File"→"Exit"でいったん終了したあと、再度wgnuplot.exeを立ち上げ、保存したファイルを"File"→"Open"で再開できます
  • 図の保存がいまいち、不親切（やり方に気付いていないだけかも知れないですけれど→やはり方法ありました）。グラフの図を表示させた状態で、そのウィンドウの左上のロゴを左クリックして"Option"→"Copy to Clipboard"で、適当なWindowsアプリケーションに貼り付けて使います。ちなみに、掲載図はMSWordにはりつけたら、プロットがおかしくなったので、Windowsアクセサリのペイントに貼り付けてJPG保存したもの
  • 図の保存はterminal　を画像用かが他のに変えて何か別の図でもなんでも描く

>gnuplot set terminal png
#epsに出力するなら　"set terminal postscript". terminal の種類は　"set terminal"と打つと例示される
>gnuplot set output "test.png"
>splot "test.data"
>set terminal win
>replot
#Windowsの場合はこうして出力端末を標準に戻すことで指定pngファイルを閉じて図ファイルの作成を完成させる。Windowsでない場合は"set terminal x11"

• 実行環境その他
  • ファイルの保存、リロードなどに際しては、日本語フォントを使用したディレクトリ名・ファイル名などはトラブルを起こすようなので、日本語フォントなしの環境での実行・保存・再利用をすることが適当と思われます

0.189604286484527 0.198897580425052 0.836431930510805
0.135755359159011 0.559375416575752 0.413132845638937
0.00755818691411057 -0.328968111745793 1.43470780629034
0.274352627154823 -0.213331161945922 0.582662914478237
-0.632028814438994 -0.249819188878758 14.0025669018830
-0.307528382571158 0.254545640532727 2.76805895225835
0.947152457611384 0.0468104113695391 16.0841937389343
0.356146989566693 -0.880018367045698 6.95354566166122
-0.854894270865624 0.377966552087912 19.8345082245713
0.337330414055211 -0.500199432051971 4.64227844223897
0.0830013799240398 0.0624558343260513 0.487522308239964

• 上述のデータの作り方
  • 集団・マーカーの存在空間次元を与え、そこに２集団のサンプリング確率密度を多次元正規分布で取る。マーカーはアレル頻度がある方向について0から1まで単調増加するように設定する。変化が急であればそれだけサンプリング領域においてもアレル頻度が変化し、変化を緩やかにすれば、サンプリング領域内ではほとんど一定になる。
  • マーカーの頻度変化方向はディリクレ分布からの乱数を使って適当にばらつかせた、単位ベクトルで与え、アレル変化の変化の強さを別の係数で与えている
  • サンプルのディプロタイプは、サンプルの位置における各マーカーのアレル頻度においてHWEを仮定している
  • 最終的には、マーカーの頻度変化ベクトルの軸ごとの成分の強さを算出し、それと、２群間における２ｘ３ディプロタイプテーブルでのトレンド検定ｐ値との関係を見ることで、マーカー頻度変化方向と集団サンプリング領域のずれベクトルとの関係を図示して説明することを主な目的としている。

#空間次元
K<-4
#マーカー数
Nx<-100
#集団数
Np<-2
#集団別のサンプル数
Ns<-rep(0,Np)
Ns[1]<-100
Ns[2]<-100
#２集団間の違い(第一軸にのみ２群の違いを入れている例）
popdif<-1
coef<-rep(0,K)
#for(i in 1:K){
#coef[i]<-popdif*i
#}
coef[1]<-popdif;
#マーカーの頻度勾配
markerV<-10

#マーカーのアレル頻度分布の平均値と分散
innerproduct<-function(a){return(a%*%a)}

LocationDirichlet<-function(Nx,K){
library(MCMCpack)
ret<-rdirichlet(Nx,rep(0.1,K))
#ret<-ret-1/K
dist<-sqrt(apply(ret,1,innerproduct))
for(i in 1:Nx){
ret[i,]<-ret[i,]/dist[i]
for(j in 1:K){
randomsign<-runif(1)
if(randomsign<0.5){
ret[i,j]<-(-1)*ret[i,j]
}

}

}
return(ret)
}

Mx<-LocationDirichlet(Nx,K)
Mx<-matrix(rnorm(Nx*K,0,1),nrow=Nx)
#Mx<-matrix(runif(Nx*K),nrow=Nx)
Vx<-matrix(10*(runif(Nx*K))^2,nrow=Nx)
#V2x<-runif(Nx)*1
V2x<-rep(markerV,Nx)
#マーカーのリスク係数
Rx<-rnorm(Nx)

#集団の分布の平均値と分散
#Mp<-matrix(runif(Np*K),nrow=Np)
Mp<-matrix(rep(0,Np*K),nrow=Np)

#Mp[1,]<-rep(-1*coef,K)
#Mp[2,]<-rep(1*coef,K)
for(i in 1:K){
Mp[1,i]<-(-1)*coef[i]
Mp[2,i]<-1*coef[i]
}

#Mp[1,2]<-(-1)*coef
#Mp[2,2]<-1*coef
#Mp<-matrix(rnorm(Np*K,0,0.01),nrow=Np)

Vp<-matrix((runif(Np*K))^2,nrow=Np)
Vp[1,]<-rep(1,K)
Vp[2,]<-rep(1,K)
#Vp<-matrix(runif(Np*K),nrow=Np)



#サンプルの位置を決める
Location<-function(ns,mp,vp){
ret<-matrix(rep(0,ns*length(mp)),nrow=ns)
for(i in 1:length(mp)){
ret[,i]<-rnorm(ns,mp[i],sqrt(vp[i]^2))
}
return(ret)
}

Ls1<-Location(Ns[1],Mp[1,],Vp[1,])
Ls2<-Location(Ns[2],Mp[2,],Vp[2,])

#サンプルの所在地でのマーカーのアレル頻度を決める
AlleleFreq<-function(L,M,V){
I<-length(L[,1])
J<-length(M[,1])
K<-length(L[1,])
ret<-matrix(rep(1,I*J),nrow=I)
for(i in 1:I){
for(j in 1:J){
for(k in 1:K){
 ret[i,j]<-ret[i,j]*pnorm(L[i,k],M[j,k],V[j,k])

}
}
}
return(ret)
}
#アレル頻度単調増減
AlleleFreq2<-function(L,M,V2){
I<-length(L[,1])
J<-length(M[,1])
K<-length(L[1,])
ret<-matrix(rep(0,I*J),nrow=I)
for(i in 1:I){
for(j in 1:J){
DifVec<-L[i,]-M[j,]
Dist<-DifVec%*%M[j,]
ret[i,j]<-atan(Dist/V2[j])/pi+0.5
}
}
return(ret)
}



#AF1<-AlleleFreq(Ls1,Mx,Vx)
#AF2<-AlleleFreq(Ls2,Mx,Vx)

AF1<-AlleleFreq2(Ls1,Mx,V2x)
AF2<-AlleleFreq2(Ls2,Mx,V2x)


#アレル頻度から、ディプロタイプをサンプリングする
SampleDiplotypeFromAlleleFreq<-function(af){
df1<-af*af
df2<-2*af*(1-af)
df3<-(1-af)*(1-af)
type<-c(0,1,2)
I<-length(af[,1])
J<-length(af[1,])
ret<-matrix(rep(0,I*J),nrow=I)
for(i in 1:I){
for(j in 1:J){
rrr<-runif(1)
if(rrr<df1[i,j]){
ret[i,j]<-0
}else if(rrr<df1[i,j]+df2[i,j]){
ret[i,j]<-1
}else{
ret[i,j]<-2
}
#ret[i,j]<-sample(type,1,replace=T,c(df1[i,j],df2[i,j],df3[1,j]))
}
}
return(ret)
}


GT1<-SampleDiplotypeFromAlleleFreq(AF1)
GT2<-SampleDiplotypeFromAlleleFreq(AF2)

#サンプルごとのリスクを計算する
R1<-GT1%*%Rx
R2<-GT2%*%Rx

#リスクのヒストグラムを描く
Hall<-hist(c(R1,R2))
hist(R1,breaks=Hall$breaks,ylim=c(0,max(Hall$counts)),col="blue",density=40,main="",xlab="",ylab="")
par(new=T)
hist(R2,breaks=Hall$breaks,ylim=c(0,max(Hall$counts)),col="red",density=40,main="",xlab="",ylab="")

#集団別のリスク平均
mean(R1)
mean(R2)

#リスクの累積分布を描く
maxrisk<-max(R1,R2)
minrisk<-min(R1,R2)
plot(sort(R1),ylim=c(minrisk,maxrisk),col="blue")
par(new=T)
plot(sort(R2),ylim=c(minrisk,maxrisk),col="red")

#サンプルの２次元位置分布を第一軸・第二軸で描く
max1<-max(Ls1[,1],Ls2[,1])
min1<-min(Ls1[,1],Ls2[,1])

max2<-max(Ls1[,2],Ls2[,2])
min2<-min(Ls1[,2],Ls2[,2])

plot(Ls1[,1],Ls1[,2],col="blue",xlim=c(min1,max1),ylim=c(min2,max2))
par(new =T)
plot(Ls2[,1],Ls2[,2],col="red",xlim=c(min1,max1),ylim=c(min2,max2))

#カウントの仕方(これは結局使わない)
count1<-apply(GT1,2,table)
count2<-apply(GT2,2,table)

#ケースコントロール集計

Rbind<-rbind(GT1,GT2)
phenotype<-c(rep(0,Ns[1]),rep(1,Ns[2]))

#それをトレンド検定
Trend<-function(Rbind,phenotype){
ret<-rep(0,length(Rbind[1,]))
for(i in 1:length(Rbind[1,])){
table<-table(Rbind[,i],phenotype)
w<-seq(from=0,to=length(table[,1])-1,by=1)
trendresult<-prop.trend.test(table[,1],table[,1]+table[,2],w)
ret[i]<-trendresult$p.value
}
return(ret)
}

trendP<-Trend(Rbind,phenotype)

#トレンド検定の累積分布を描く
plot(sort(trendP))
thomas<-seq(from=1/(length(trendP)+1),to=length(trendP)/(length(trendP)+1),by=1/(length(trendP)+1))

plot(thomas,sort(trendP),xlim=c(0,1),ylim=c(0,1))
plot(log(thomas),log(sort(trendP)),xlim=c(min(log(trendP)),0),ylim=c(min(log(trendP)),0))

#plot(AF1,AF2)

min<-min(Ls1,Ls2)
max<-max(Ls1,Ls2)
lim<-c(min,max)
plot(Mx[,1],Mx[,2],xlim=lim,ylim=lim)
par(new=T)
plot(Ls1[,1],Ls1[,2],col="blue",xlim=lim,ylim=lim)
par(new =T)
plot(Ls2[,1],Ls2[,2],col="red",xlim=lim,ylim=lim)

#plot(atan(Mx[,2]/Mx[,1]),trendP)

#plot(atan(Mx[,2]/Mx[,1]),log(trendP))


MxLen<-apply(Mx,1,innerproduct)
MxLen<-sqrt(MxLen)
divideByVector<-function(mat,vec){
ret<-mat
for(i in 1:length(vec)){
ret[i,]<-mat[i,]/vec[i]
}
return(ret)
}

tmpAcos<-divideByVector(Mx,MxLen)

Acos<-acos(tmpAcos)/pi-0.5

plot(Acos[,1],log(trendP),xlim=c(-0.5,0.5))
plot(Acos[,2],log(trendP),xlim=c(-0.5,0.5))


plot(Mx[,1],trendP)

plot(sort(trendP))

#トレンド検定の累積分布を描く
plot(sort(trendP))
thomas<-seq(from=1/(length(trendP)+1),to=length(trendP)/(length(trendP)+1),by=1/(length(trendP)+1))

plot(thomas,sort(trendP),xlim=c(0,1),ylim=c(0,1))
plot(log(thomas),log(sort(trendP)),xlim=c(min(log(trendP)),0),ylim=c(min(log(trendP)),0))


library(scatterplot3d)
scatterplot3d(Mx[,1],Mx[,2],-log(trendP),highlight.3d=T,type="h",angle=75)

binder<-cbind(Mx[,1],Mx[,2,-log(trendP))
write.table(binder,file="data.txt",row.names=FALSE,col.names=FALSE)