1,555
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to search830 - MarkovModel::Impute() (view source)
Revision as of 14:50, 25 September 2013
, 14:50, 25 September 2013Created page with "<source lang="cpp"> void MarkovModel::Impute(char * major, char * observed, float * probs, char ** haplotypes, float ** freqs, int position) { d..."
<source lang="cpp">
void MarkovModel::Impute(char * major, char * observed, float * probs,
char ** haplotypes, float ** freqs, int position)
{
double P[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
for (int i = 0; i < states; i++)
P[haplotypes[i][position]] += probs[i];
double ptotal = P[1] + P[2] + P[3] + P[4];
double pmajor = P[major[position]];
int mle = 1;
for (int i = 2; i <= 4; i++)
if (P[i] >= P[mle])
mle = i;
char labels[] = {0, 'a', 'c', 'g', 't'};
imputedDose[position] += imputedHap[position] = (pmajor / ptotal);
imputedAlleles[position] = labels[mle];
// What would we have imputed if we didn't observe a genotype at the current position?
// This assumes that freqs[i][position] stores the frequency of allele i at position (for i > 0)
// and is equal to 1.0 when the current genotype is not observed (i = 0)
double fmatch = 1.0 / (1. - E[position] + E[position] * freqs[observed[position]][position]);
double fmismatch = 1.0 / (E[position] * freqs[observed[position]][position]);
for (int i = 1; i <= 4; i++)
if (observed[position] == i)
P[i] *= fmatch;
else
P[i] *= fmismatch;
ptotal = P[1] + P[2] + P[3] + P[4];
pmajor = P[major[position]];
leaveOneOut[position] = pmajor / ptotal;
}
</source>
void MarkovModel::Impute(char * major, char * observed, float * probs,
char ** haplotypes, float ** freqs, int position)
{
double P[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
for (int i = 0; i < states; i++)
P[haplotypes[i][position]] += probs[i];
double ptotal = P[1] + P[2] + P[3] + P[4];
double pmajor = P[major[position]];
int mle = 1;
for (int i = 2; i <= 4; i++)
if (P[i] >= P[mle])
mle = i;
char labels[] = {0, 'a', 'c', 'g', 't'};
imputedDose[position] += imputedHap[position] = (pmajor / ptotal);
imputedAlleles[position] = labels[mle];
// What would we have imputed if we didn't observe a genotype at the current position?
// This assumes that freqs[i][position] stores the frequency of allele i at position (for i > 0)
// and is equal to 1.0 when the current genotype is not observed (i = 0)
double fmatch = 1.0 / (1. - E[position] + E[position] * freqs[observed[position]][position]);
double fmismatch = 1.0 / (E[position] * freqs[observed[position]][position]);
for (int i = 1; i <= 4; i++)
if (observed[position] == i)
P[i] *= fmatch;
else
P[i] *= fmismatch;
ptotal = P[1] + P[2] + P[3] + P[4];
pmajor = P[major[position]];
leaveOneOut[position] = pmajor / ptotal;
}
</source>
