PedigreeAlleleFreq.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 "PedigreeAlleleFreq.h"
#include "QuickIndex.h"
#include "Error.h"
#include <math.h>
int CountAlleles(Pedigree & /* ped */, int marker)
{
// With automatic recoding in the pedigree loader there
// is no need to iterate through the pedigree ...
MarkerInfo * info = Pedigree::GetMarkerInfo(marker);
return info->CountAlleles();
}
void LumpAlleles(Pedigree & ped, int marker, double threshold, bool reorder)
{
// find out how many alleles there are
int alleles = ped.CountAlleles(marker);
if (alleles < 2) return;
MarkerInfo * info = PedigreeGlobals::GetMarkerInfo(marker);
if (alleles < info->freq.Length())
alleles = info->freq.Length() - 1;
IntArray counts(alleles + 1);
counts.Zero();
// Count number of occurrences for each allele
for (int i = 0; i < ped.count; i++)
{
counts[int(ped[i].markers[marker][0])]++;
counts[int(ped[i].markers[marker][1])]++;
}
// Calculate treshold for lumping alleles
int total = 0;
for (int i = 1; i <= alleles; i++)
total += counts[i];
int thresh = int(total * threshold);
// If threshold is set at zero, we artificially increase
// counts for alleles that do not appear in the pedigree
// but whose frequencies are set > 0.0. This ensures that
// allele frequency data does not get discarded when simply
// recoding alleles (vs. lumping)
if (thresh == 0)
for (int i = 1; i < info->freq.Length(); i++)
if (counts[i] == 0 && info->freq[i] > 0.0)
counts[i] = 1, total++;
// If allele reordering is disabled, put in dummy allele
// counts so as to ensure that allele have desired ordering
if (!reorder)
{
QuickIndex index(info->alleleLabels);
index.Reverse();
for (int i = 0; i < index.Length(); i++)
counts[index[i]] = i + 1;
total = counts.Sum(1, counts.Length() - 1);
}
// Order all alleles according to their frequency
// Zero, which corresponds to missing values, stays put!
counts[0] = total + 1;
QuickIndex index(counts);
index.Reverse();
// recode alleles
// all alleles where frequency < thresh are labelled N
// use counts array to keep track of labels
int N = 0;
bool rare = false;
for (int i = 0; i <= alleles; i++)
if (counts[index[i]] > thresh)
{
counts[index[i]] = i;
N++;
}
else
{
if (counts[index[i]] > 0)
rare = true;
counts[index[i]] = N;
}
// This loop does the recoding
for (int i = 0; i < ped.count; i++)
{
Alleles & current = ped[i].markers[marker];
current[0] = counts[current[0]];
current[1] = counts[current[1]];
}
StringArray oldLabels(info->alleleLabels);
String label;
info->alleleLabels.Clear();
info->alleleNumbers.Clear();
for (int i = 0; i < N; i++)
{
if (oldLabels.Length() <= index[i])
info->alleleLabels.Push(label = index[i]);
else
info->alleleLabels.Push(oldLabels[index[i]]);
if (i) info->alleleNumbers.SetInteger(info->alleleLabels.Last(), i);
}
// Reorder allele frequencies if necessary
if (info->freq.Length())
{
Vector freq(info->freq);
info->freq.Dimension(N);
info->freq[0] = 0.0;
for (int i = 1; i < N; i++)
{
info->freq[i] = freq[index[i]];
freq[index[i]] = 0;
}
if ((1.0 - info->freq.Sum()) > 1e-10)
rare = true;
if (rare)
{
info->freq.Dimension(N + 1);
info->freq[N] = 1.0 - info->freq.Sum();
}
}
if (rare)
{
info->alleleLabels.Push("OTHER");
info->alleleNumbers.SetInteger("OTHER", info->alleleLabels.Length());
}
}
bool EstimateFrequencies(Pedigree & ped, int marker, int estimator)
{
int alleleCount = CountAlleles(ped, marker);
IntArray founder(alleleCount + 1);
IntArray all(alleleCount + 1);
founder.Zero();
all.Zero();
for (int i = 0; i < ped.count; i++)
{
// When counting alleles, note that males only carry one X chromosome
// and are arbitrarily scored as homozygous.
all[ped[i].markers[marker][0]]++;
if (!ped.chromosomeX || ped[i].sex != SEX_MALE)
all[ped[i].markers[marker][1]]++;
if (!ped[i].isFounder()) continue;
founder[ped[i].markers[marker][0]]++;
if (!ped.chromosomeX || ped[i].sex != SEX_MALE)
founder[ped[i].markers[marker][1]]++;
}
MarkerInfo * info = ped.GetMarkerInfo(marker);
if (info->freq.dim > 0)
{
// previous allele frequency information is available
if (alleleCount >= info->freq.dim)
error("For marker %s, input files define %d alleles, but at least\n"
"one other allele (named '%s') occurs in the pedigree\n",
(const char *) info->name, info->freq.dim - 1,
(const char *) info->GetAlleleLabel(alleleCount));
for (int i = 1; i <= alleleCount; i++)
if (all[i] > 0 && info->freq[i] <= 0.0)
error("Although allele %s for marker %s has frequency zero,\n"
"it occurs %d times in the pedigree",
(const char *) info->GetAlleleLabel(i), (const char *) info->name, all[i]);
return false;
}
else
{
if (alleleCount < 1)
{
// If no one is genotyped, default to two equifrequent allele
// since some programs do not like monomorphic markers
info->freq.Dimension(3);
info->freq[0] = 0.0;
info->freq[1] = 0.99999;
info->freq[2] = 0.00001;
return true;
}
info->freq.Dimension(alleleCount + 1);
info->freq.Zero();
if (estimator == FREQ_FOUNDERS && founder.Sum() > founder[0])
{
// Make sure the frequency of alleles occuring in the pedigree
// is never zero
for (int i = 1; i <= alleleCount; i++)
if (founder[i] == 0 && all[i] > 0)
founder[i] = 1;
// To get frequencies, just multiply counts by 1 / total_counts
double factor = 1.0 / (founder.Sum() - founder[0]);
for (int i = 1; i <= alleleCount; i++)
info->freq[i] = founder[i] * factor;
}
else if (estimator == FREQ_ALL || estimator == FREQ_FOUNDERS)
{
// To get frequencies, just multiply counts by 1 / total_counts
double factor = 1.0 / (all.Sum() - all[0]);
for (int i = 1; i <= alleleCount; i++)
info->freq[i] = all[i] * factor;
}
else if (estimator == FREQ_EQUAL)
// Assume all alleles have equal frequency
{
// Count the number of observed alleles
all[0] = 0;
int alleles = all.CountIfGreater(0);
double freq = 1.0 / alleles;
// Set equal frequencies for all occuring alleles
for (int i = 0; i <= alleleCount; i++)
info->freq[i] = all[i] ? freq : 0.0;
}
}
return true;
}
Alleles
Definition: PedigreeAlleles.h:24
IntArray
Definition: IntArray.h:24
MarkerInfo
Definition: PedigreeGlobals.h:30
Pedigree
Definition: Pedigree.h:33
QuickIndex
Definition: QuickIndex.h:28
StringArray
Definition: StringArray.h:24
String
Definition: StringBasics.h:39
Vector
Definition: MathVector.h:29