■
関連する記事はこちら。
まず、集団の２SNP、４ハプロタイプの初期頻度を与え、集団サイズを固定したうえで、次世代のハプロタイプ型は当該世代のハプロタイプ構成頻度からランダムに選ぶことで浮動をシミュレートする。
シミュレート世代数を与え、各世代のハプロタイプ頻度と、それを正四面体表現するための座標に置換する。
これを後掲のjavaで行い、座標データを出力したファイルをRの関数(こちらも後掲)で処理して、３次元表示する。
public class DrawLDTetrahedron {
/*
* Tetrahedron
*/
public static double[] v1 = {0,0,Math.sqrt(6)/4};
public static double[] v2 = {Math.sqrt(3)/3, 0, -Math.sqrt(6)/12};
public static double[] v3 = {-Math.sqrt(3)/6, 0.5, -Math.sqrt(6)/12};
public static double[] v4 = {-Math.sqrt(3)/6, -0.5, -Math.sqrt(6)/12};
public static double[] v12 = {(v1[0]+v2[0])/2,(v1[1]+v2[1])/2,(v1[2]+v2[2])/2};
public static double[] v13 = {(v1[0]+v3[0])/2,(v1[1]+v3[1])/2,(v1[2]+v3[2])/2};
public static double[] v14 = {(v1[0]+v4[0])/2,(v1[1]+v4[1])/2,(v1[2]+v4[2])/2};
public static double[] v23 = {(v2[0]+v3[0])/2,(v2[1]+v3[1])/2,(v2[2]+v3[2])/2};
public static double[] v24 = {(v2[0]+v4[0])/2,(v2[1]+v4[1])/2,(v2[2]+v4[2])/2};
public static double[] v34 = {(v3[0]+v4[0])/2,(v3[1]+v4[1])/2,(v3[2]+v4[2])/2};
public static double[] vec13_24 ={(v13[0]-v24[0])/2,(v13[1]-v24[1])/2,(v13[2]-v24[2])/2};
public static double[] vec12_34 ={(v12[0]-v34[0])/2,(v12[1]-v34[1])/2,(v12[2]-v34[2])/2};
public static double[] vec14_23 ={(v14[0]-v23[0])/2,(v14[1]-v23[1])/2,(v14[2]-v23[2])/2};
public static double[] LocInTetrahedron(double[] four){
double[] ret = new double[3];
double[] psi = PsiFromHapFour(four);
ret[0]=psi[0]*vec13_24[0]+psi[1]*vec12_34[0]+psi[2]*vec14_23[0];
ret[1]=psi[0]*vec13_24[1]+psi[1]*vec12_34[1]+psi[2]*vec14_23[1];
ret[2]=psi[0]*vec13_24[2]+psi[1]*vec12_34[2]+psi[2]*vec14_23[2];
return ret;
}
public static double[] PsiFromHapFour(double[] four){
double[] ret = {(four[0]+four[1]-four[2]-four[3]),(four[0]+four[2]-four[1]-four[3]),(four[0]+four[3]-four[1]-four[2])};
return ret;
}
}
public static void main(String[] args) throws IOException {
String dir = "C:\\Documents and Settings\\Yamada\\デスクトップ\\June8";
String outfile = "DriftSimHapFreq2.txt";
String outfile2 = "DriftSimPhi2.txt";
File file = new File(dir,outfile);
File file2 = new File(dir,outfile2);
BufferedWriter[] bw = new BufferedWriter[2];
bw[0]=null;
bw[1]=null;
bw[0] = new BufferedWriter(new FileWriter(file));
bw[1] = new BufferedWriter(new FileWriter(file2));
double[] init = {0.25,0.25,0.25,0.25};
double[] loc = Test.LocInTetrahedron(init);
String st = ""+init[0]+"\t"+init[1]+"\t"+init[2]+"\t"+init[3]+"\n";
String st2 = "" + loc[0]+"\t"+loc[1]+"\t"+loc[2]+"\n";
InOut4.out3File(bw[0],st);
InOut4.out3File(bw[1],st2);
System.out.println(""+"\t"+init[0]+"\t"+init[1]+"\t"+init[2]+"\t"+init[3]);
init=MiscUtil.accumstandard(init);
//System.out.println(""+"\t"+init[0]+"\t"+init[1]+"\t"+init[2]+"\t"+init[3]);
/*
MiscUtil.standard(init);
double[] accum = new double[init.length];
accum[accum.length-1]=1;
for(int i=accum.length-2;i>-1;i--){
accum[i]=accum[i+1]-init[i+1];
}
*/
int numchrom=1000;
int numgeneration=10000;
int seedx = 1000000;
int seed = (int)(Math.random()*seedx);
Utils.MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);
for(int i=0;i<numgeneration;i++){
double[] tmp = new double[init.length];
for(int j=0;j<numchrom;j++){
double r = mz.nextDouble();
for(int k=0;k<init.length;k++){
if(r<=init[k]){
tmp[k]++;
break;
}
}
}
init=MiscUtil.accumstandard(tmp);
double[] tmp2 = {init[0],init[1]-init[0],init[2]-init[1],init[3]-init[2]};
System.out.println(""+i+"\t"+tmp2[0]+"\t"+tmp2[1]+"\t"+tmp2[2]+"\t"+tmp2[3]);
st = tmp2[0]+"\t"+tmp2[1]+"\t"+tmp2[2]+"\t"+tmp2[3]+"\n";
loc = Test.LocInTetrahedron(tmp2);
st2 = "" + loc[0]+"\t"+loc[1]+"\t"+loc[2]+"\n";
InOut4.out3File(bw[0],st);
InOut4.out3File(bw[1],st2);
}
bw[0].close();
bw[1].close();
}
Rの関数
plotTetrahedron<-function(angle=40,scaley=1,File){
library(scatterplot3d)
v1<-c(0,0,0.612372436)/2
v2<-c(0.577350269, 0, -0.204124145)/2
v3<-c(-0.288675135, 0.5, -0.204124145)/2
v4<-c(-0.288675135, -0.5, -0.204124145)/2
kizami<-seq(0,1,0.01)
L12x<-v1[1]*kizami+v2[1]*(1-kizami)
L12y<-v1[2]*kizami+v2[2]*(1-kizami)
L12z<-v1[3]*kizami+v2[3]*(1-kizami)
L13x<-v1[1]*kizami+v3[1]*(1-kizami)
L13y<-v1[2]*kizami+v3[2]*(1-kizami)
L13z<-v1[3]*kizami+v3[3]*(1-kizami)
L14x<-v1[1]*kizami+v4[1]*(1-kizami)
L14y<-v1[2]*kizami+v4[2]*(1-kizami)
L14z<-v1[3]*kizami+v4[3]*(1-kizami)
L23x<-v2[1]*kizami+v3[1]*(1-kizami)
L23y<-v2[2]*kizami+v3[2]*(1-kizami)
L23z<-v2[3]*kizami+v3[3]*(1-kizami)
L24x<-v2[1]*kizami+v4[1]*(1-kizami)
L24y<-v2[2]*kizami+v4[2]*(1-kizami)
L24z<-v2[3]*kizami+v4[3]*(1-kizami)
L34x<-v3[1]*kizami+v4[1]*(1-kizami)
L34y<-v3[2]*kizami+v4[2]*(1-kizami)
L34z<-v3[3]*kizami+v4[3]*(1-kizami)
LimX<-c(-0.3,0.7)
LimY<-c(-0.6,0.6)
LimZ<-c(-0.3,0.7)
ScaleY=scaley
AAngle <- angle
scatterplot3d(L12x,L12y,L12z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="x",ylab="y",zlab="z")
par(new=T)
scatterplot3d(L13x,L13y,L13z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(L14x,L14y,L14z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(L23x,L23y,L23z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(L24x,L24y,L24z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(L34x,L34y,L34z,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
data<-read.table(file=File)
scatterplot3d(data$V1,data$V2,data$V3,xlim = LimX, ylim = LimY, zlim = LimZ,type="l",angle=AAngle,scale.y=ScaleY,highlight.3d=TRUE,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(data$V1[1],data$V2[1],data$V3[1],xlim = LimX, ylim = LimY, zlim = LimZ,type="p",angle=AAngle,scale.y=ScaleY,xlab="",ylab="",zlab="")
par(new=T)
scatterplot3d(data$V1[length(data)/3],data$V2[length(data)/3],data$V3[length(data)/3],xlim = LimX, ylim = LimY, zlim = LimZ,type="p",angle=AAngle,scale.y=ScaleY,xlab="",ylab="",zlab="")
}
4ハプロタイプ頻度から正四面体座標への変換関数(R)
LocInTetrahedron<-function(four){
psi<-c(four[1]+four[2]-four[3]-four[4],four[1]+four[3]-four[2]-four[4],four[1]+four[4]-four[2]-four[3])
v1<-c(0,0,0.612372436)
v2<-c(0.577350269, 0, -0.204124145)
v3<-c(-0.288675135, 0.5, -0.204124145)
v4<-c(-0.288675135, -0.5, -0.204124145)v12 <-c((v1[1]+v2[1])/2,(v1[2]+v2[2])/2,(v1[3]+v2[3])/2)
v13 <-c((v1[1]+v3[1])/2,(v1[2]+v3[2])/2,(v1[3]+v3[3])/2)
v14 <-c((v1[1]+v4[1])/2,(v1[2]+v4[2])/2,(v1[3]+v4[3])/2)
v23 <-c((v2[1]+v3[1])/2,(v2[2]+v3[2])/2,(v2[3]+v3[3])/2)
v24 <-c((v2[1]+v4[1])/2,(v2[2]+v4[2])/2,(v2[3]+v4[3])/2)
v34 <-c((v3[1]+v4[1])/2,(v3[2]+v4[2])/2,(v3[3]+v4[3])/2)
vec13_24 <-c((v13[1]-v24[1])/2,(v13[2]-v24[2])/2,(v13[3]-v24[3])/2)
vec12_34 <-c((v12[1]-v34[1])/2,(v12[2]-v34[2])/2,(v12[3]-v34[3])/2)
vec14_23 <-c((v14[1]-v23[1])/2,(v14[2]-v23[2])/2,(v14[3]-v23[3])/2)
ret<-c(psi[1]*vec13_24[1]+psi[2]*vec12_34[1]+psi[3]*vec14_23[1],psi[1]*vec13_24[2]+psi[2]*vec12_34[2]+psi[3]*vec14_23[2],psi[1]*vec13_24[3]+psi[2]*vec12_34[3]+psi[3]*vec14_23[3])return(ret)
}
