libStatGen Software 1
PedigreeTwin.cpp
1/*
2 * Copyright (C) 2010 Regents of the University of Michigan
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18#include "Pedigree.h"
19#include "Error.h"
20
21#include <stdio.h>
22
23bool Pedigree::TwinCheck()
24{
25 bool fail = false;
26 IntArray mzTwins;
27
28 for (int f = 0; f < familyCount; f++)
29 {
30 mzTwins.Clear();
31
32 for (int i = families[f]->first, j; i <= families[f]->last; i++)
33 // Is this person an identical twin?
34 if (persons[i]->isMzTwin(*persons[i]))
35 {
36 // Have we got another identical sib yet?
37 for (j = 0; j < mzTwins.Length(); j++)
38 if (persons[i]->isMzTwin(*persons[mzTwins[j]]))
39 break;
40
41 // If not, add to list of twins
42 if (j == mzTwins.Length())
43 {
44 mzTwins.Push(i);
45 continue;
46 }
47
48 // Check that their genotypes are compatible and
49 // merge new twin's genotypes into original twin...
50 Person * original = persons[mzTwins[j]];
51 Person * twin = persons[i];
52
53 for (int m = 0; m < Person::markerCount; m++)
54 {
55 if (!original->markers[m].isKnown())
56 original->markers[m] = twin->markers[m];
57 else if (twin->markers[m].isKnown() &&
58 twin->markers[m] != original->markers[m])
59 printf("MZ Twins %s and %s in family %s have "
60 "different %s genotypes\n",
61 (const char *) original->pid,
62 (const char *) twin->pid,
63 (const char *) original->famid,
64 (const char *) Person::markerNames[m]),
65 fail = true;
66
67 if (twin->sex != original->sex)
68 printf("MZ Twins %s and %s in family %s have "
69 "different sexes\n",
70 (const char *) original->pid,
71 (const char *) twin->pid,
72 (const char *) original->famid),
73 fail = true;
74 }
75 }
76
77 if (mzTwins.Length() == 0) continue;
78
79 // In the second pass we copy merged twin genotypes
80 // from original twin to other twins
81 for (int i = families[f]->first, j; i <= families[f]->last; i++)
82 if (persons[i]->isMzTwin(*persons[i]))
83 {
84 for (j = 0; j < mzTwins.Length(); j++)
85 if (persons[i]->isMzTwin(*persons[mzTwins[j]]))
86 break;
87
88 if (mzTwins[j] == i) continue;
89
90 Person * original = persons[mzTwins[j]];
91 Person * twin = persons[i];
92
93 for (int m = 0; m < Person::markerCount; m++)
94 twin->markers[m] = original->markers[m];
95 }
96 }
97 return fail;
98}
99
100void Pedigree::MergeTwins()
101{
102 if (!haveTwins) return;
103
104 IntArray mzTwins, surplus;
105
106 printf("Merging MZ twins into a single individual...\n");
107
108 for (int f = 0; f < familyCount; f++)
109 {
110 mzTwins.Clear();
111
112 for (int i = families[f]->first, j; i <= families[f]->last; i++)
113 if (persons[i]->isMzTwin(*persons[i]))
114 {
115 // Have we got another identical sib yet?
116 for (j = 0; j < mzTwins.Length(); j++)
117 if (persons[i]->isMzTwin(*persons[mzTwins[j]]))
118 break;
119
120 // If not, add to list of twins
121 if (j == mzTwins.Length())
122 {
123 mzTwins.Push(i);
124 continue;
125 }
126
127 // Append name to first twins name
128 persons[mzTwins[j]]->pid += ((char) '$') + persons[i]->pid;
129
130 // Set the first twin to affected if at least one of the cotwins is affected
131 for (int j = 0; j < affectionCount; j++)
132 if (persons[i]->affections[j] == 2)
133 persons[mzTwins[j]]->affections[j] = 2;
134
135 surplus.Push(i);
136 }
137
138 // More than one representative of each twin-pair...
139 if (surplus.Length())
140 {
141 // Reassign parent names for each offspring
142 for (int i = families[f]->first, j; i < families[f]->last; i++)
143 if (!persons[i]->isFounder())
144 {
145 if (persons[i]->father->isMzTwin(*persons[i]->father))
146 {
147 for (j = 0; j < mzTwins.Length(); j++)
148 if (persons[i]->father->isMzTwin(*persons[mzTwins[j]]))
149 break;
150 persons[i]->fatid = persons[mzTwins[j]]->pid;
151 }
152 if (persons[i]->mother->isMzTwin(*persons[i]->mother))
153 {
154 for (j = 0; j < mzTwins.Length(); j++)
155 if (persons[i]->mother->isMzTwin(*persons[mzTwins[j]]))
156 break;
157 persons[i]->motid = persons[mzTwins[j]]->pid;
158 }
159 }
160
161 // Delete surplus individuals
162 while (surplus.Length())
163 {
164 int serial = surplus.Pop();
165
166 delete persons[serial];
167
168 for (count--; serial < count; serial++)
169 persons[serial] = persons[serial + 1];
170 }
171
172 // Resort pedigree
173 Sort();
174 }
175 }
176}
177
178
179
180