830 - Haplotyper::RandomSetup()
From Genome Analysis Wiki
Jump to navigationJump to searchvoid Haplotyper::RandomSetup(Random * rand)
{
// accesses two arrays...
// genotypes[i][j] includes genotypes for each individual
// i ranges from 0 .. individuals - 1
// j ranges from 0 .. markers - 1
// coded as 0 (missing), 1, 2, 3
//
// haplotypes[i][j] includes haplotypes for each individual
// i ranges from 0 .. 2 * individuals - 1
// coded as 0 or 1, no missing data
for (int j = 0; j < markers; j++)
{
int alleles = 0, mac = 0;
for (int i = 0; i < individuals; i++)
if (genotypes[i][j] != GENOTYPE_MISSING)
alleles += 2, mac += genotypes[i][j] - 1;
if (alleles == 0)
{
for (int i = 0; i < individuals; i++)
haplotypes[i * 2][j] = haplotypes[i * 2 + 1][j] = 0;
continue;
}
float freq = mac / (float) alleles;
for (int i = 0; i < individuals; i++)
switch (genotypes[i][j])
{
case 0: /* MISSING GENOTYPE */
haplotypes[i * 2][j] = rand->Next() < freq;
haplotypes[i * 2 + 1][j] = rand->Next() < freq;
break;
case 1: /* HOMOZYGOUS FOR ALLELE 1 */
haplotypes[i * 2][j] = 0;
haplotypes[i * 2 + 1][j] = 0;
break;
case 2: /* HETEROZYGOUS */
{
int bit = rand->Binary();
haplotypes[i * 2][j] = bit;
haplotypes[i * 2 + 1][j] = bit ^ 1;
}
break;
case 3: /* HOMOZYGOUS FOR ALLELE 2 */
haplotypes[i * 2][j] = 1;
haplotypes[i * 2 + 1][j] = 1;
break;
}
}
}