連鎖不平衡の強さに応じた2SNPハプロタイプ頻度の特定



データシミュレーションなどをするときには、LDの強さに応じてハプロタイプ頻度を与える必要がある。

一番、頻用するのは、2SNP間のLDの強さに応じて4ハプロタイプ頻度を与えるような場合である。

LDインデックスの1つであるr^2は、2SNPのアレル頻度が与えられているときに、その2アレル間の相関の強さを0−1で表すものである。

今、第1SNPのアレル頻度をp、第2SNPのそれをqとすると

¥delta = ¥sqrt{r^2¥times pq(1-p)(1-q)}とすれば、h1=pq-¥delta,h2=p(1-q)+¥delta,h3=(1-p)q+¥delta,h4=(1-p)(1-q)-¥deltaとなる。r^2,p,qの与え方によっては、適当なハプロタイプ頻度が存在しないこともあること、逆に2パターンのハプロタイプ頻度セットがありえることに注意する。

上限・下限は¥frac{(r^2)(1-p)}{(r^2)(1-p)+p}¥le q ¥le ¥frac{1-p}{(r^2)p+(1-p)},または¥frac{(r^2)p}{(r^2)p+(1-p)}¥le q ¥le ¥frac{p}{(r^2) (1-p) + p}

r^2に基づく4ハプロタイプ頻度の算出にあたって、qのとりうる範囲の上限・下限においては、r^2が1にならないが、D’

1となっている。また、qのとりうる範囲が2つの離れた範囲となるか、隣接(連結)した範囲となるかはrsq=¥frac{p}{1-p},rsq=¥frac{1-p}{p}を境としている。

これらを確かめるエクセルはこちら


public static int LDtransGenotype(double p,double[] afqrange,int x,double d,MersenneTwisterFast mz){
/*
* pは与えられたSNP1のアレル頻度
* xはgenotypeの番地
* dはr^2
*
*/
int ret = x;
double[][] rangeQ = RangeQForHapFreq(p,d);

rangeQ = NarrowRangeDual(rangeQ,afqrange);
double[] length = new double[2];
length[0] = rangeQ[0][1]-rangeQ[0][0];
length[1] = rangeQ[1][1]-rangeQ[1][0];
double q=1-p;
int type = 0;
if(length[0]==0){//lengty[0]=length[1]=0
q=1-p;
}else{
for(int i=0;i<length.length;i++){
if(length[i]<0){
length[i]=0;
}
}
double sumlen = length[0]+length[1];
double rand = mz.nextDouble();

if(rand<length[0]/sumlen){
type = 0;
}else{
type =1;
}
q = rangeQ[type][0] + mz.nextDouble()*(rangeQ[type][1]-rangeQ[type][0]);
}
//System.out.println("p"+"\t"+p);
//System.out.println("r"+"\t"+d);
//System.out.println("q\t" + q);
//double q = afqrange[0] + mz.nextDouble()*(afqrange[1]-afqrange[0]);//SNP2のアレル頻度
double k = CoefForCalcHapFreq(p,q,d);
//System.out.println("k\t" + k);
double[] hfreq = new double[4];
if(type==0){
hfreq[0]=p*q-k;
hfreq[1]=p*(1-q)+k;
hfreq[2]=(1-p)*q+k;
hfreq[3]=(1-p)*(1-q)-k;
}else if(type==1){
hfreq[0]=p*q+k;
hfreq[1]=p*(1-q)-k;
hfreq[2]=(1-p)*q-k;
hfreq[3]=(1-p)*(1-q)+k;
}


double[] cumulhfreq = new double[4];
cumulhfreq[0]=hfreq[0];
for(int i=1;i<hfreq.length;i++){
cumulhfreq[i]=cumulhfreq[i-1]+hfreq[i];
}
double[] r = {mz.nextDouble(),mz.nextDouble()};
/*
for(int i=0;i<hfreq.length;i++){
System.out.println("hfreq\t"+i+"\t"+hfreq[i]);
}
for(int i=0;i<cumulhfreq.length;i++){
System.out.println("cumulhfreq\t"+i+"\t"+cumulhfreq[i]);
}
*/

//a1 a2はxから出す2アレルを表す
int[] a ={0,0};
int[] h ={0,0};//2chromosomesのハプロタイプID
if(x==0){
a[0]=0;a[1]=0;
}else if(x==1){
a[0]=0;a[1]=1;
}else if(x==2){
a[0]=1;a[1]=1;
}
for(int i=0;i<a.length;i++){
double thres=0;
if(a[i]==0){
//h1かh2
thres = hfreq[0]/p;
}else if(a[i]==1){//h3 or h4
thres = hfreq[2]/(1-p);
}
if(a[i]==0){
if(r[i]<thres){
h[i]=0;
}else{
h[i]=1;
}
}else if(a[i]==1){
if(r[i]<thres){
h[i]=2;
}else{
h[i]=3;
}
}
}
if((h[0]==0 || h[0]==2)&&(h[1]==0 || h[1]==2)){
ret = 0;
}else if((h[0]==1 || h[0]==3)&&(h[1]==1 || h[1]==3)){
ret = 2;
}else{
ret = 1;
}
return ret;
}

public static double CoefForCalcHapFreq(double p,double q,double r){
double ret=Math.pow(r*p*q*(1-p)*(1-q), 0.5);
return ret;
}
public static double[][] RangeQForHapFreq(double p,double r){
double[][] ret = new double[2][2];
ret[0][0] = r*(1-p)/((r*(1-p)+p));
ret[0][1] = (1-p)/(r*p+(1-p));
ret[1][0] = r*p/(r*p+(1-p));
ret[1][1] = p/(r*(1-p)+p);
return ret;
}
public static double[] NarrowRange(double[] a,double[] b){
/*
* aのセグメントのうち
* bを満足する範囲を返す
*/
double[] ret = {a[0],a[1]};
if(a[0]<b[0]){
ret[0]=b[0];
}
if(a[1]>b[1]){
ret[1]=b[1];
}
//System.out.println("narrowed "+ret[0] + " " + ret[1]);
return ret;

}
public static double[][] NarrowRangeDual(double[][] a,double[] b){
double[][] ret = new double[2][2];
for(int i=0;i<ret.length;i++){
ret[i] = NarrowRange(a[i],b);
}
return ret;
}