Allelic associationのこと



Allelic associationとは、あるサンプル集合について、2箇所の多型のアレルの間に相関があることを言う。一般集団においては、連鎖不平衡による相関もあり、これはDNA上の物理的距離との関係が強い。また、集団に構造がある場合には、構造の存在によりAllelic associationが認められる。このような場合のAllelic associationは物理的な距離とは関係ないので、異なる染色体上に存在する2多型の間にも認められる。

連鎖不平衡によるAllelic associationが存在する2多型における検定統計量の間には相関が生じる。したがって、連鎖不平衡にある多型を含むような多マーカー検定においては、そのすべてのマーカーが相互に独立であるという前提でのマルチプルテスティング補正は保守的にすぎる。したがって、独立とみなせるマーカーの数がわかれば、それに基づいて、マルチプルテスティング補正を行うのが妥当である。連鎖不平衡はゲノムワイドに存在するが、ケース群・コントロール群間で差がないとみなせるとき、これら多マーカーでの検定統計量には、理論的な分布をとる。

他方、集団構造がある場合には、たまたま、ケースサンプルとコントロールサンプルとの構造が一致すれば、検定統計量は理論的な分布をとるが、サンプリングバイアスによりケース群とコントロール群の構造に差が生じてくると、統計量は理論値をよりも偏った値をとる(P値に換算して小さくなる)。そのもっとも極端な場合は、集団が2亜集団から成り立ち、ケースのサンプルがすべて1つの亜集団由来でコントロールのサンプルがすべてもう片方の亜集団由来となってしまったような場合である。このような場合、個々のマーカーにおいて、ケース・コントロールが均一であるという帰無仮説に基づく検定統計量の偏りは最も大きくなる。

今、構造化のある集団からのサンプルについて、複数のマーカーを用いて分割表検定を行う場合を考える。2通りの極端な仮説について検討する。

(1)構造化のない集団からのランダムなサンプルであるとする

(2)構造化が存在し、かつ、構造化は2亜集団によって形成され、ケース・コントロールのサンプルがまさに、異なる亜集団からそれぞれすべてサンプルされた(しまった)とする

総人数N、ケース人数N1、コントロール人数N2、因子の保有数X、非保有数Yのような周辺度数を持つ分割表a,b,c,d(a+b=N1,c+d=N2,a+c=X,b+d=Y)が得られたとする。

(1)の場合に照らすと、¥frac{N1! N2! X! Y!}{N!a!b!c!d!}なる確率であることがわかる。この式を(2)の場合と比べるために、少し変更する。今、N1,N2の条件でa,b,c,dを観測する確率は、帰無仮説における因子Xの保有率の最尤推定p_0=¥frac{X}{N}を用いて、¥frac{N1!N2!}{a!b!c!d!}p_0^{X}(1-p_0){Y}と表せる。今、N,N1,N2,X,Yなる周辺度数を持つすべての分割表、a',b',c',d'について¥frac{N1!N2!}{a’!b’!c’!d’!}p_0^{X}(1-p_0){Y}を計算してやる。これらの総和が1となるように補正をした値が、それぞれの分割表を得る確率である。このうち、a,b,c,dを観測する確率以下であるものを足し合わせた値がフィッシャーの正確確率検定P値である。

(2)の場合を考えよう。今、対立仮説として、ケースとコントロールとは、全くことなる2亜群の正確な代表であるから、それぞれに異なる因子保有頻度を想定するべきである。その最尤推定値はp_1=¥frac{a}{N1}p_2=¥frac{c}{N2}である。このような場合に、分割表a,b,c,dを観測する確率は¥frac{N1!N2!}{a!b!c!d!}p_1^{a}(1-p_1){b}p_2^{c}(1-p_2)^{d}に比例する。同様に、それ以外の分割表の観測確率は¥frac{N1!N2!}{a’!b’!c’!d’!}p_1^{a}(1-p_1){b}p_2^{c}(1-p_2)^{d}に比例する。すべての分割表が網羅されているので、この総和が1となるように補正してやると、観測データから求めた対立仮説における最尤推定値をもとにした、それぞれの分割表の観測確率が求められた。

今、対立仮説が成り立つかもしれないと思いつつも、帰無仮説の棄却検定(普通のケースコントロール分割表検定)を行ったとする。フィッシャーの正確確率検定を行ったとすれば、(1)で記したような値のP値が算出されるが、その確率は、対立仮説が成り立っているものとすると、その値は、フィッシャーのP値そのものではなくなってくる(帰無仮説が成り立っているときには、あるフィッシャーのP値を観測する確率はそのPの値そのものである)。冒頭で、構造化があるときには、検定統計量がP値として小さいほうにずれると書いた。今、対立仮説に基づいたフィッシャーのP値の観測確率から、フィッシャーのP値の期待値を求めると、確かに、0.5よりも小さくなることがわかる。

そんな具合を確かめるためのjavaソースとエクセルはこちら

3,7,7,3なる2x2表が観測されたとき、これと同じ周辺度数をもつ分割表は11通りあり、この11枚の分割表にから得られるフィッシャーのPの値は 1.08251E-05 0.001093334 0.023014138 0.178895408 0.656281799 1 の6通りあることが示されている。そして、その6通りのExactP値につき、対立仮説の最尤推定保有率に基づく観測確率は0.004129263 0.075843702 0.282102601 0.368865727 0.214242907 0.0548158であり、帰無仮説に基づく観測確率が 1.08251E-05 0.001082509 0.021920804 0.15588127 0.477386391 0.343718201 であるのに比べ、小さなP値での観測確率が多き無なっていることがわかる。それを反映して、ExactP値の期待値は、対立仮説のときの 0.267983221 に対し、帰無仮説の 0.685410316 という差も生じている。




/*
* 作成日: 2007/01/19
*
* TODO この生成されたファイルのテンプレートを変更するには次へジャンプ:
* ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート
*/
package StatUtils;

/**
* @author ryamada
*
* TODO この生成された型コメントのテンプレートを変更するには次へジャンプ:
* ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート
*/
public class AlternativeFisher {
public int[][] data;
public int[][] marginal;
public int max;
//public double[] logFactList;
public int[][][] tables;
public double[] prob;
public double[][] ordif;

public ProbabilityDensity[] pd;

public static void main(String[] args) {
int[][] d = {{93,7},{97,3}};

int[][] varRange = Fisher.RangeFor2xN(d);
for(int i=0;i<varRange.length;i++){
for(int j=0;j<varRange[i].length;j++){
System.out.println(""+i+" "+j+" "+varRange[i][j]);
}
}

//int[][] marginal=Fisher.marginal(d);
AlternativeFisher cf = new AlternativeFisher(d);
ProbabilityDensity.PrintProbabilityDensity(cf.pd[0],"\t","\n");
ProbabilityDensity.PrintProbabilityDensity(cf.pd[1],"\t","\n");
//PrintCubicFisher(cf,"\t","\n");
}

public AlternativeFisher(int[][] d){
data = MiscUtil.DeepCopyInt2(d);
marginal = Fisher.marginal(d);
for(int i=0;i<marginal.length;i++){
//for(int j=0;j<marginal[i].length;j++){
//System.out.println("marginal "+marginal[i][j]);
//}

}
max = 0;
for(int i=0;i<marginal[0].length;i++){
max+=marginal[0][i];
}
double[] logFactList = Fisher.serialLogfact(max);
int[][][] varRange2=Fisher.RangeForMxNFromMarginal(marginal);
tables = Fisher.makeTables(marginal,varRange2);

double constant = 0;
constant += logFactList[max];
for(int i=0;i<marginal.length;i++){
for(int j=0;j<marginal[i].length;j++){
constant -= logFactList[marginal[i][j]];
}

}
double[] pertable = new double[tables.length];
prob = new double[tables.length];
for(int i=0;i<pertable.length;i++){
double tmp = 0;
for(int j=0;j<tables[i].length;j++){
for(int k=0;k<tables[i][j].length;k++){
//System.out.println("table cell "+ tables[i][j][k] + " "+logFactList[tables[i][j][k]]);
tmp += logFactList[tables[i][j][k]];
}
}
prob[i] = -(constant + tmp);
prob[i]=Math.exp(prob[i]);
//System.out.println("probtable "+prob[i]);
}
double[][] fraction = new double[d.length][d[0].length];
for(int i=0;i<fraction.length;i++){
for(int j=0;j<fraction[i].length;j++){
fraction[i][j]=(double)d[i][j]/(double)marginal[0][i];
//System.out.println("fraction " + fraction[i][j]);
}
}
double[] probAlternative = new double[tables.length];
for(int i=0;i<probAlternative.length;i++){
probAlternative[i] = 0;
for(int j=0;j<d.length;j++){
for(int k=0;k<d[j].length;k++){
if(fraction[j][k]!=0){
probAlternative[i] += tables[i][j][k]*Math.log(fraction[j][k]);
}

}
}
probAlternative[i] = Math.exp(probAlternative[i]);
//System.out.println("ppp "+ probAlternative[i]);
//probAlternative[i] += pertable[i];
//probAlternative[i] = Math.exp(probAlternative[i]);
//System.out.println("ppp "+ probAlternative[i]);
}
MiscUtil.standard(probAlternative);
for(int i=0;i<probAlternative.length;i++){
probAlternative[i]*=prob[i];
//System.out.println(" "+probAlternative[i]);
}
MiscUtil.standard(probAlternative);
//for(int i=0;i<probAlternative.length;i++){
//probAlternative[i]*=prob[i];
//System.out.println(""+probAlternative[i]);
//}
int[][] Porder = MiscUtil.bublSortRetInt2(prob);
//for(int i=0;i<Porder.length;i++){
//System.out.println("Porder "+Porder[i][1] + " "+prob[Porder[i][1]]);
//}
double[] cumulProb=new double[Porder.length];
double[] cumulAltP=new double[Porder.length];
for(int i=0;i<prob.length;i++){
for(int j=Porder[i][1];j<prob.length;j++){
cumulProb[j]+=prob[i];
cumulAltP[j]+=probAlternative[i];
}
}
/*
for(int i=0;i<cumulProb.length;i++){
System.out.println("cumulProb "+ cumulProb[i]);
}
for(int i=0;i<cumulAltP.length;i++){
System.out.println("cumulAltP "+ cumulAltP[i]);
}
*/
double[] NRprobExactP = {cumulProb[0]};
double[] NRprobAltP = {cumulAltP[0]};
for(int i=1;i<cumulProb.length;i++){
if(cumulProb[i]!=cumulProb[i-1]){
NRprobExactP = MiscUtil.AddElemDouble1(NRprobExactP,cumulProb[i]);
NRprobAltP = MiscUtil.AddElemDouble1(NRprobAltP,cumulAltP[i]);
}
}
for(int i=0;i<NRprobExactP.length;i++){
NRprobExactP[i]/=NRprobExactP[NRprobExactP.length-1];
}
pd=new ProbabilityDensity[2];
int type=1;
pd[0]=new ProbabilityDensity(NRprobExactP,NRprobExactP,type);
pd[1]=new ProbabilityDensity(NRprobExactP,NRprobAltP,type);
}