関連SNPの統計量空間の狭さ:ラベルパーミュテーション用エクセル2



昨日、ラベルパーミュテーションをエクセルで簡易に行ってみた。

今日は、相互に独立でない2SNP(連鎖不平衡など)において同様にエクセルでラベルパーミュテーションし、その統計量を比較・プロットしてみた。

これは、2SNPのジェノタイプの近さと、2SNPからの関連統計量がとりうる統計量2次元空間の広さ、密度の具合についてのおおまかなアイディアを得るためのツールである。

エクセルはこちら(重めです)

Javaソース版はこちら

"org.apache.commons.math.distribution.ChiSquaredDistributionImpl"については、こちらを参照


public class SNPpairSim {

int[][][] two3x3;
int N;
int[] mergL;

int numperm;
/**
* @param args
*/
public static void main(String[] args) {
// TODO 自動生成されたメソッド・スタブ
int[][][] two3x3 = {{{50,4,1},{26,4,1},{12,2,0}},
{{40,4,1},{14,4,1},{5,2,0}}};

int numperm=1000;
SNPpairSim sss = new SNPpairSim(two3x3,numperm);

}
public SNPpairSim(int[][][] two3x3_,int numperm_){
two3x3=StatUtils.MiscUtil.DeepCopyInt3(two3x3_);
N=NtotalFrom2x3x3(two3x3);
numperm=numperm_;
int[] sh = new int[N];
int[][] tmpgen = new int[sh.length][2];
int counter=0;
for(int i=0;i<two3x3.length;i++){
for(int j=0;j<two3x3[i].length;j++){
for(int k=0;k<two3x3[i][j].length;k++){
for(int l=0;l<two3x3[i][j][k];l++){
tmpgen[counter][0]=j;
tmpgen[counter][1]=k;
counter++;
}

}
}
}
mergL=mergLFrom2x3x3(two3x3);
for(int i=0;i<N;i++){
if(i<mergL[0]){
sh[i]=0;
}else{
sh[i]=1;
}
}
double[] psOr = Ps(tmpgen,sh);
System.out.println("OR\t"+psOr[0]+"\t"+psOr[1]);
int x=100000000;
int seed = (int)(Math.random()*x);
Utils.MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);
for(int i=0;i<numperm;i++){
DoubleLdPlot.PermutationOP.shuffleMZ(sh,mz);
double[] psPerm = Ps(tmpgen,sh);
System.out.println("Perm\t"+psPerm[0]+"\t"+psPerm[1]);
}

}
public static double[] Ps(int[][] tmpgen,int[] sh){
double[] ret = new double[2];


double[][] twoBYthree1 = new double[2][3];
double[][] twoBYthree2 = new double[2][3];

for(int i=0;i<tmpgen.length;i++){
twoBYthree1[sh[i]][tmpgen[i][0]]++;
twoBYthree2[sh[i]][tmpgen[i][1]]++;
}

String tmp="";
for(int i=0;i<twoBYthree1.length;i++){
for(int j=0;j<twoBYthree1[i].length;j++){
tmp+=twoBYthree1[i][j]+"\t";
}
tmp+="\n";
}
//System.out.println(tmp);
tmp="";
for(int i=0;i<twoBYthree2.length;i++){
for(int j=0;j<twoBYthree2[i].length;j++){
tmp+=twoBYthree2[i][j]+"\t";
}
tmp+="\n";
}
//System.out.println(tmp);
ret[0] = BioBankPerm2.SingleSNPCalculator.trendChi(twoBYthree1);
ret[1] = BioBankPerm2.SingleSNPCalculator.trendChi(twoBYthree2);

int df=1;//自由度
org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 =
new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df);
try{
ret[0] = 1-chidf1.cumulativeProbability(ret[0]);
ret[1] = 1-chidf1.cumulativeProbability(ret[1]);

} catch (org.apache.commons.math.MathException ex) {

}


return ret;
}
public static int NtotalFrom2x3x3(int[][][] two3x3){
int N=0;
for(int i=0;i<two3x3.length;i++){
for(int j=0;j<two3x3[i].length;j++){
for(int k=0;k<two3x3[i][j].length;k++){
N+=two3x3[i][j][k];
}

}
}
return N;
}
public static int[] mergLFrom2x3x3(int[][][] two3x3){
int[] ret = new int[2];
for(int i=0;i<two3x3.length;i++){
for(int j=0;j<two3x3[i].length;j++){
for(int k=0;k<two3x3[i][j].length;k++){
ret[i]+=two3x3[i][j][k];
}
}
}
return ret;
}
public static int[] mergLFrom2x3(int[][] twoBYthree){
int[] mergL = {twoBYthree[0][0]+twoBYthree[0][1]+twoBYthree[0][2],
twoBYthree[1][0]+twoBYthree[1][1]+twoBYthree[1][2]};

return mergL;
}
public static int[] mergCFrom2x3(int[][] twoBYthree){
int[] mergC = {twoBYthree[0][0]+twoBYthree[1][0],
twoBYthree[0][1]+twoBYthree[1][1],
twoBYthree[0][2]+twoBYthree[1][2]};
return mergC;
}

}