構造化集団のシミュレーション

Genomic controlのレビューを紹介した(記事はこちら)。その中で、構造化集団の構造化の程度をfで表し、それを用いて、集団のHWEから外れたジェノタイプ頻度の計算と、構成亜集団のアレル頻度のばらつき具合についての記載があった。

平均アレル頻度pに対して、fが与えられたときに、ジェノタイプ頻度は。

また、亜集団のアレル頻度は期待値p、分散の正規分布でシミュレートする話し。

MersenneTwisterFastのnextGaussian()関数を用いて次のようにシミュレートしてみる。

public class FstAlleleGenotypeFreq {
	public static void main(String[] args) {

		double fst = 0.01;

		double meanp = 0.3;

		for(int i=0;i<args.length;i++){

			if(args[i].equals("-meanp")){

				meanp = Double.parseDouble(args[i+1]);

			}

			if(args[i].equals("-fst")){

				fst = Double.parseDouble(args[i+1]);

				

				

			}

		}

		System.out.println("genotype frequency");

		

		double[] genotypefreq = GenFreqFst(meanp,fst);

		System.out.println(genotypefreq[0] + "\t" + genotypefreq[1] + "\t" + genotypefreq[2]);

		int numSim = 100;

		System.out.println("allele freq simulation");

		double[] randomallelefreq = new double[numSim];

		for(int i=0;i<numSim;i++){

			randomallelefreq[i] = AlleleFreqRandom(meanp,fst);

			System.out.println(randomallelefreq[i]);

		}

		double averageP=Calculator.mean(randomallelefreq);

		double var = Calculator.var(randomallelefreq);

		System.out.println("given meanp and mean of simulated allelefreq\t" + meanp + "\t"+averageP);

		double calculatedvarP=fst * meanp*(1-meanp);

		System.out.println("calculated varP and var of simulated allelefreq\t" + calculatedvarP + "\t"+var);

		

		

	}

	

	public static double[] GenFreqFst(double meanp,double fst){

		double[] ret = new double[3];

		ret[0] = fst * meanp + (1-fst)*meanp*meanp;

		ret[1] = (1-fst)*2*meanp*(1-meanp);

		ret[2] = fst * (1-meanp) + (1-fst)*(1-meanp)*(1-meanp);

		return ret;

	

	}

	

	public static double AlleleFreqRandom(double meanp,double fst){

		double ret = 0;

		int x =1000000;

		int seed = (int)((double)x*Math.random());

		MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);

		double std = Math.pow(meanp*(1-meanp)*fst,0.5);

		ret = meanp+std*mz.nextGaussian();

		return ret;

		

	}

}