Loading...
Searching...
No Matches
GenotypeLists.cpp
/*
* Copyright (C) 2010 Regents of the University of Michigan
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "GenotypeLists.h"
// When the next line is uncommented, the genotype elimination routines
// produce a lot of output useful for debugging
// #define DEBUG_ELIMINATOR
GenotypeList::GenotypeList()
{
ignore = false;
}
bool GenotypeList::EliminateGenotypes(Pedigree & ped, Family * family, int marker)
{
// First, allocate a genotype list for the family
GenotypeList * list = new GenotypeList [family->count];
// Next, update the possible allele lists for each individual
InitializeList(list, ped, family, marker);
// Then, do multiple rounds of elimination until a problem is found
// or no changes are made
#ifdef DEBUG_ELIMINATOR
Print(list, ped, family, marker);
#endif
while (PairwiseCheck(list, ped, family) || FamilyCheck(list, ped, family))
#ifdef DEBUG_ELIMINATOR
Print(list, ped, family, marker)
#endif
;
for (int i = 0; i < family->count; i++)
if (!list[i].ignore && list[i].allele1.Length() == 0)
{
printf("%s - Family %s has a subtle genotype inconsistency\n",
(const char *) ped.markerNames[marker], (const char *) family->famid);
delete [] list;
return false;
}
delete [] list;
return true;
}
void GenotypeList::InitializeList(GenotypeList * list, Pedigree & ped, Family * family, int marker)
{
for (int i = family->count - 1; i >= 0; i--)
{
Person & person = ped[family->path[i]];
int id = person.traverse;
bool maleX = person.sex == SEX_MALE && ped.chromosomeX;
#ifdef DEBUG_ELIMINATOR
printf("Initializing genotype list for %s ...\n", (const char *) person.pid);
#endif
// If an individual is genotyped ...
if (person.markers[marker].isKnown())
{
// Their genotype list starts with just one entry!
list[id].Dimension(1);
list[id].SetGenotype(0, person.markers[marker][0], person.markers[marker][1]);
list[id].alleles.Clear();
list[id].alleles.Push(person.markers[marker][0]);
list[id].alleles.PushIfNew(person.markers[marker][1]);
list[id].ignore = false;
// "Heterozygous" males have no possible genotypes
if (maleX && person.markers[marker].isHeterozygous())
list[id].Dimension(0);
}
else if (list[id].alleles.Length())
if (person.sex == SEX_MALE && ped.chromosomeX)
{
// Males only carry one X chromosome
list[id].Dimension(list[id].alleles.Length() + 1);
for (int i = 0, out = 0; i < list[id].alleles.Length(); i++)
list[id].SetGenotype(out++, list[id].alleles[i], list[id].alleles[i]);
list[id].SetGenotype(list[id].alleles.Length(), -1, -1);
list[id].ignore = false;
}
else
{
// Build the genotype list based on the available allele lists
int count = list[id].alleles.Length() * (list[id].alleles.Length() + 3) / 2 + 1;
list[id].Dimension(count);
for (int i = 0, out = 0; i < list[id].alleles.Length(); i++)
{
// Allow for all pairs of "transmitted" alleles
for (int j = 0; j <= i; j++)
list[id].SetGenotype(out++, list[id].alleles[i], list[id].alleles[j]);
// Allow for an unstransmitted allele
list[id].SetGenotype(out++, list[id].alleles[i], -1);
}
// Allow for a pair of untransmitted alleles
list[id].SetGenotype(count - 1, -1, -1);
list[id].ignore = false;
}
else
list[id].ignore = true;
// If the individual is a founder this is all there is to it
if (i < family->founders) continue;
// If the individual is not a founder, update the parental genotype lists...
int fatid = person.father->traverse;
int motid = person.mother->traverse;
for (int i = 0; i < list[id].alleles.Length(); i++)
{
list[motid].alleles.PushIfNew(list[id].alleles[i]);
if (!maleX) list[fatid].alleles.PushIfNew(list[id].alleles[i]);
}
}
}
bool GenotypeList::PairwiseCheck(GenotypeList * list, Pedigree & ped, Family * family)
{
#ifdef DEBUG_ELIMINATOR
printf("Checking Relative Pairs ...\n");
#endif
bool changed = false;
for (int i = family->count - 1; i >= family->founders; i--)
{
Person & person = ped[family->path[i]];
int id = person.traverse;
int fatid = person.father->traverse;
int motid = person.mother->traverse;
bool maleX = person.sex == SEX_MALE && ped.chromosomeX;
if (list[id].ignore) continue;
// Check if genotypes are consistent with paternal genotypes
for (int i = 0; i < list[id].allele1.Length(); i++)
{
int al1 = list[id].allele1[i];
int al2 = list[id].allele2[i];
// Remove offspring genotypes incompatible with parental genotypes
if ((maleX && !list[motid].Matches(al1) && al1 != -1) ||
(!maleX && !(al1 == -1 && al2 == -1) &&
!(list[fatid].Matches(al1) && (al2 == -1 || list[motid].Matches(al2))) &&
!((al2 == -1 || list[fatid].Matches(al2)) && list[motid].Matches(al1))))
{
list[id].Delete(i--);
changed = true;
}
}
// The offspring genotype list allows for a wild-card untransmitted allele
// so any single parental genotype is possible
if (list[id].Matches(-1))
continue;
// Check if genotypes are consistent with offspring genotypes
for (int i = 0; i < list[motid].allele1.Length(); i++)
{
int al1 = list[motid].allele1[i];
int al2 = list[motid].allele2[i];
// Remove genotypes incompatible with offspring genotypes
if (!list[id].Matches(al1) &&
!list[id].Matches(al2))
{
list[motid].Delete(i--);
changed = true;
}
}
// Males don't affect genotype lists for their fathers
if (maleX) continue;
// Check if genotypes are consistent with offspring genotypes
for (int i = 0; i < list[fatid].allele1.Length(); i++)
{
int al1 = list[fatid].allele1[i];
int al2 = list[fatid].allele2[i];
// Remove genotypes incompatible with offspring genotypes
if (!list[id].Matches(al1) &&
!list[id].Matches(al2))
{
list[fatid].Delete(i--);
changed = true;
}
}
#ifdef DEBUG_ELIMINATOR
printf("Done checking individual %s\n", (const char *) person.pid);
Print(list, ped, family, 0);
#endif
}
return changed;
}
bool GenotypeList::FamilyCheck(GenotypeList * list, Pedigree & ped, Family * family)
{
#ifdef DEBUG_ELIMINATOR
printf("Checking Nuclear Families ...\n");
#endif
bool changed = false;
for (int i = family->count - 1; i >= family->founders; i--)
{
Person & person = ped[family->path[i]];
int fatid = person.father->traverse;
int motid = person.mother->traverse;
// Only go through the loop once per sibship
if (person.sibs[0] != &person || list[fatid].ignore || list[motid].ignore)
continue;
#ifdef DEBUG_ELIMINATOR
printf("Checking Sibship with %s ...\n", (const char *) person.pid);
#endif
// Reset checked genotypes for the mother, father and child
list[fatid].checked = 0;
list[motid].checked = 0;
for (int i = 0; i < person.sibCount; i++)
list[person.sibs[i]->traverse].checked = 0;
// Go through each of the paternal genotypes
changed |= TrimParent(list, person, fatid, motid);
// Go through each of maternal genotypes
changed |= TrimParent(list, person, motid, fatid);
// Sort out the unchecked offspring genotypes ...
for (int i = 0; i < person.sibCount; i++)
{
int sibid = person.sibs[i]->traverse;
bool maleX = person.sibs[i]->sex == SEX_MALE && ped.chromosomeX;
// For dealing with male X chromosomes, the pairwise check is all we need
if (maleX) continue;
for (int j = list[sibid].checked; j < list[sibid].allele1.Length(); j++)
changed |= Cleanup(list, person, motid, fatid, sibid, j);
}
#ifdef DEBUG_ELIMINATOR
// Print(list, ped, family, 0);
#endif
}
return changed;
}
bool GenotypeList::Matches(int genotype, int allele)
{
return allele1[genotype] == allele || allele2[genotype] == allele;
}
bool GenotypeList::Matches(int allele)
{
return allele1.Find(allele) != -1 || allele2.Find(allele) != -1;
}
int GenotypeList::SaveGenotype(int genotype)
{
if (checked > genotype)
return genotype;
if (checked != genotype)
{
allele1.Swap(genotype, checked);
allele2.Swap(genotype, checked);
}
return checked++;
}
bool GenotypeList::CheckTrio(GenotypeList * list, int fatid, int motid, int child,
int i, int j, int k)
{
// TODO: add tests for this code...
return (list[fatid].Matches(i, list[child].allele1[k]) &&
(list[motid].Matches(j, list[child].allele2[k]) || list[child].allele2[k] == -1)) ||
((list[fatid].Matches(i, list[child].allele2[k]) || list[child].allele2[k] == -1) &&
list[motid].Matches(j, list[child].allele1[k])) ||
(list[child].allele1[k] == -1 && list[child].allele2[k] == -1);
}
void GenotypeList::Dimension(int genotypes)
{
allele1.Dimension(genotypes);
allele2.Dimension(genotypes);
}
void GenotypeList::SetGenotype(int genotype, int al1, int al2)
{
allele1[genotype] = al1;
allele2[genotype] = al2;
}
void GenotypeList::Delete(int genotype)
{
allele1.Delete(genotype);
allele2.Delete(genotype);
}
bool GenotypeList::TrimParent(GenotypeList * list, Person & person, int motid, int fatid)
{
bool trimmed = false;
while (list[motid].checked < list[motid].allele1.Length())
{
int current = list[motid].allele1.Length() - 1;
bool saved = false;
// Pair it with each possible paternal genotype
for (int i = list[fatid].allele1.Length() - 1; i >= 0; i--)
{
int matches = 0;
// Find out if the pairing is compatible with at least one genotype for each child
for (int j = 0; j < person.sibCount; j++)
{
int sibid = person.sibs[j]->traverse;
int maleX = person.sibs[j]->sex == SEX_MALE && person.chromosomeX;
// Since we have done the pairwise check, there is nothing more
// to do for males ...
if (list[sibid].ignore || maleX)
{
matches++;
continue;
}
for (int k = list[sibid].allele1.Length() - 1; k >= 0; k--)
if (CheckTrio(list, motid, fatid, sibid, current, i, k))
{
matches++;
break;
}
if (matches != j + 1)
break;
}
// Save maternal and paternal genotypes, mark all compatible sibling genotypes
if (matches == person.sibCount)
{
for (int j = 0; j < person.sibCount; j++)
{
int sibid = person.sibs[j]->traverse;
for (int k = list[sibid].checked; k < list[sibid].allele1.Length(); k++)
if (CheckTrio(list, motid, fatid, sibid, current, i, k))
list[sibid].SaveGenotype(k);
}
list[motid].SaveGenotype(current);
list[fatid].SaveGenotype(i);
saved = true;
break;
}
}
if (!saved)
{
list[motid].Delete(current);
trimmed = true;
}
}
return trimmed;
}
bool GenotypeList::Cleanup(GenotypeList * list, Person & person, int motid, int fatid, int child, int geno)
{
for (int current = 0; current < list[motid].allele1.Length(); current++)
for (int i = list[fatid].allele1.Length() - 1; i >= 0; i--)
if (CheckTrio(list, motid, fatid, child, current, i, geno))
{
int matches = 0;
// Find out if the pairing is compatible with at least one genotype for each child
for (int j = 0; j < person.sibCount; j++)
{
int sibid = person.sibs[j]->traverse;
int maleX = person.sibs[j]->sex == SEX_MALE && person.chromosomeX;
// After completing the pairwise check, all males are guaranteed
// to be compatible with their mothers
if (list[sibid].ignore || maleX)
{
matches++;
continue;
}
for (int k = list[sibid].allele1.Length() - 1; k >= 0; k--)
if (CheckTrio(list, motid, fatid, sibid, current, i, k))
{
matches++;
break;
}
if (matches != j + 1)
break;
}
// Update list of compatible sibling genotypes
if (matches == person.sibCount)
for (int j = 0; j < person.sibCount; j++)
{
int sibid = person.sibs[j]->traverse;
for (int k = list[sibid].checked; k < list[sibid].allele1.Length(); k++)
if (CheckTrio(list, motid, fatid, sibid, current, i, k))
list[sibid].SaveGenotype(k);
return false;
}
}
list[child].Delete(geno);
return true;
}
void GenotypeList::Print(GenotypeList * list, Pedigree & ped, Family * family, int marker)
{
MarkerInfo * info = ped.GetMarkerInfo(marker);
for (int i = 0; i < family->count; i++)
{
printf("%s - ", (const char *) ped[family->path[i]].pid);
for (int j = 0; j < list[i].allele1.Length(); j++)
{
if (list[i].allele1[j] == -1)
printf("*/");
else
printf("%s/", (const char *) info->GetAlleleLabel(list[i].allele1[j]));
if (list[i].allele2[j] == -1)
printf("* ");
else
printf("%s ", (const char *) info->GetAlleleLabel(list[i].allele2[j]));
}
printf("\n");
}
printf("\n");
}
Family
Definition PedigreeFamily.h:28
GenotypeList
Definition GenotypeLists.h:24
MarkerInfo
Definition PedigreeGlobals.h:30
Pedigree
Definition Pedigree.h:33
Person
Definition PedigreePerson.h:32