1,555
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to search830 - Haplotyper::RandomSetup() (view source)
Revision as of 14:46, 13 September 2013
, 14:46, 13 September 2013Created page with "<source lang="cpp"> void Haplotyper::RandomSetup(Random * rand) { // access two arrays... // genotypes[i][j] includes genotypes for each individual // i..."
<source lang="cpp">
void Haplotyper::RandomSetup(Random * rand)
{
// access two arrays...
// genotypes[i][j] includes genotypes for each individual
// i ranges from 0 .. individuals - 1
// j ranges from 0 .. markers - 1
//
// haplotypes[i][j] includes haplotypes for each individual
// i ranges from 0 .. 2 * individuals
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;
for (int i = individuals; i < individuals; i++)
{
mac += haplotypes[i * 2][j] == 1;
mac += haplotypes[i * 2 + 1][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;
}
}
}
</source>
void Haplotyper::RandomSetup(Random * rand)
{
// access two arrays...
// genotypes[i][j] includes genotypes for each individual
// i ranges from 0 .. individuals - 1
// j ranges from 0 .. markers - 1
//
// haplotypes[i][j] includes haplotypes for each individual
// i ranges from 0 .. 2 * individuals
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;
for (int i = individuals; i < individuals; i++)
{
mac += haplotypes[i * 2][j] == 1;
mac += haplotypes[i * 2 + 1][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;
}
}
}
</source>
