18#include "PedigreeAlleleFreq.h"
19#include "QuickIndex.h"
24int CountAlleles(
Pedigree & ,
int marker)
28 MarkerInfo * info = Pedigree::GetMarkerInfo(marker);
30 return info->CountAlleles();
33void LumpAlleles(
Pedigree & ped,
int marker,
double threshold,
bool reorder)
36 int alleles = ped.CountAlleles(marker);
38 if (alleles < 2)
return;
40 MarkerInfo * info = PedigreeGlobals::GetMarkerInfo(marker);
42 if (alleles < info->freq.Length())
43 alleles = info->freq.Length() - 1;
49 for (
int i = 0; i < ped.count; i++)
51 counts[int(ped[i].markers[marker][0])]++;
52 counts[int(ped[i].markers[marker][1])]++;
57 for (
int i = 1; i <= alleles; i++)
59 int thresh = int(total * threshold);
67 for (
int i = 1; i < info->freq.Length(); i++)
68 if (counts[i] == 0 && info->freq[i] > 0.0)
69 counts[i] = 1, total++;
78 for (
int i = 0; i < index.Length(); i++)
79 counts[index[i]] = i + 1;
81 total = counts.Sum(1, counts.Length() - 1);
86 counts[0] = total + 1;
95 for (
int i = 0; i <= alleles; i++)
96 if (counts[index[i]] > thresh)
103 if (counts[index[i]] > 0)
105 counts[index[i]] = N;
109 for (
int i = 0; i < ped.count; i++)
111 Alleles & current = ped[i].markers[marker];
112 current[0] = counts[current[0]];
113 current[1] = counts[current[1]];
119 info->alleleLabels.Clear();
120 info->alleleNumbers.Clear();
122 for (
int i = 0; i < N; i++)
124 if (oldLabels.Length() <= index[i])
125 info->alleleLabels.Push(label = index[i]);
127 info->alleleLabels.Push(oldLabels[index[i]]);
129 if (i) info->alleleNumbers.SetInteger(info->alleleLabels.Last(), i);
133 if (info->freq.Length())
137 info->freq.Dimension(N);
140 for (
int i = 1; i < N; i++)
142 info->freq[i] = freq[index[i]];
146 if ((1.0 - info->freq.Sum()) > 1e-10)
151 info->freq.Dimension(N + 1);
152 info->freq[N] = 1.0 - info->freq.Sum();
158 info->alleleLabels.Push(
"OTHER");
159 info->alleleNumbers.SetInteger(
"OTHER", info->alleleLabels.Length());
163bool EstimateFrequencies(
Pedigree & ped,
int marker,
int estimator)
165 int alleleCount = CountAlleles(ped, marker);
173 for (
int i = 0; i < ped.count; i++)
177 all[ped[i].markers[marker][0]]++;
178 if (!ped.chromosomeX || ped[i].sex != SEX_MALE)
179 all[ped[i].markers[marker][1]]++;
180 if (!ped[i].isFounder())
continue;
181 founder[ped[i].markers[marker][0]]++;
182 if (!ped.chromosomeX || ped[i].sex != SEX_MALE)
183 founder[ped[i].markers[marker][1]]++;
186 MarkerInfo * info = ped.GetMarkerInfo(marker);
188 if (info->freq.dim > 0)
191 if (alleleCount >= info->freq.dim)
192 error(
"For marker %s, input files define %d alleles, but at least\n"
193 "one other allele (named '%s') occurs in the pedigree\n",
194 (
const char *) info->name, info->freq.dim - 1,
195 (
const char *) info->GetAlleleLabel(alleleCount));
197 for (
int i = 1; i <= alleleCount; i++)
198 if (all[i] > 0 && info->freq[i] <= 0.0)
199 error(
"Although allele %s for marker %s has frequency zero,\n"
200 "it occurs %d times in the pedigree",
201 (
const char *) info->GetAlleleLabel(i), (
const char *) info->name, all[i]);
211 info->freq.Dimension(3);
213 info->freq[1] = 0.99999;
214 info->freq[2] = 0.00001;
218 info->freq.Dimension(alleleCount + 1);
221 if (estimator == FREQ_FOUNDERS && founder.Sum() > founder[0])
225 for (
int i = 1; i <= alleleCount; i++)
226 if (founder[i] == 0 && all[i] > 0)
230 double factor = 1.0 / (founder.Sum() - founder[0]);
232 for (
int i = 1; i <= alleleCount; i++)
233 info->freq[i] = founder[i] * factor;
235 else if (estimator == FREQ_ALL || estimator == FREQ_FOUNDERS)
238 double factor = 1.0 / (all.Sum() - all[0]);
240 for (
int i = 1; i <= alleleCount; i++)
241 info->freq[i] = all[i] * factor;
243 else if (estimator == FREQ_EQUAL)
248 int alleles = all.CountIfGreater(0);
249 double freq = 1.0 / alleles;
252 for (
int i = 0; i <= alleleCount; i++)
253 info->freq[i] = all[i] ? freq : 0.0;