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 
23 bool 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 
100 void 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 
IntArray
Definition: IntArray.h:24
Person
Definition: PedigreePerson.h:32