1 // $Id: phenotypes.cpp,v 1.5 2011/03/07 06:08:49 bobgian Exp $
2
3 /*
4 Copyright 2006 Lucian Smith, Mary Kuhner, Jon Yamato and Joseph Felsenstein
5
6 This software is distributed free of charge for non-commercial use
7 and is copyrighted. Of course, we do not guarantee that the software
8 works, and are not responsible for any damage you may cause or have.
9 */
10
11 #include <cassert>
12
13 #include "errhandling.h"
14 #include "locuscell.h" //Have to include this to use Haplotypes
15 #include "mathx.h"
16 #include "phenotypes.h"
17 #include "registry.h"
18 #include "stringx.h"
19 #include "vectorx.h"
20 #include "xml_strings.h"
21
22 using namespace std;
23
24 //------------------------------------------------------------------------------------
25
26 typedef map<multiset<string>, pair<StringVec1d, DoubleVec1d> >::iterator PhenMapIter;
27 typedef map<multiset<string>, pair<StringVec1d, DoubleVec1d> >::const_iterator PhenMapConstIter;
28
29 //------------------------------------------------------------------------------------
30
Phenotypes(long regnum,string lname)31 Phenotypes::Phenotypes(long regnum, string lname)
32 : m_regionnum(regnum),
33 m_locusname(lname),
34 m_phenomap(),
35 m_hapsmap()
36 {
37 }
38
39 //Needed a blank copy in the locus--we replace it if we have a real one.
40 //
41 //Boy howdy does this all need to be refactored. Whew, it stinks.
Phenotypes(string lname)42 Phenotypes::Phenotypes(string lname)
43 : m_regionnum(FLAGLONG),
44 m_locusname(lname),
45 m_phenomap(),
46 m_hapsmap()
47 {
48 }
49
AddPhenotype(const StringVec1d & alleles,string name,double penetrance)50 void Phenotypes::AddPhenotype(const StringVec1d& alleles, string name, double penetrance)
51 {
52 assert (m_regionnum != FLAGLONG);
53 multiset<string> alleleSet = VecToSet(alleles);
54 PhenMapIter genotype = m_phenomap.find(alleleSet);
55 if (genotype == m_phenomap.end())
56 {
57 StringVec1d names;
58 DoubleVec1d penetrances;
59 m_phenomap.insert(make_pair(alleleSet, make_pair(names, penetrances)));
60 genotype = m_phenomap.find(alleleSet);
61 assert(genotype != m_phenomap.end());
62 }
63 (*genotype).second.first.push_back(name);
64 (*genotype).second.second.push_back(penetrance);
65
66 }
67
ChooseHaplotypes(const StringVec1d & alleles)68 Haplotypes Phenotypes::ChooseHaplotypes(const StringVec1d& alleles)
69 {
70 multiset<string> alleleSet = VecToSet(alleles);
71 PhenMapIter genotype = m_phenomap.find(alleleSet);
72 if (genotype == m_phenomap.end())
73 {
74 string msg = "Unable to find any phenotypes for the set of alleles \"";
75 for (size_t allele = 0; allele<alleles.size(); allele++)
76 {
77 msg += " " + alleles[allele];
78 }
79 msg += " \".";
80 throw data_error(msg);
81 }
82 DoubleVec1d penetrances = (*genotype).second.second; //for convenience
83 if (penetrances.size() == 1)
84 {
85 return GetHaplotypes((*genotype).second.first[0]);
86 }
87 double rand = registry.GetRandom().Float();
88 double sum = 0;
89 for (size_t which=0; which < penetrances.size(); which++)
90 {
91 sum += penetrances[which];
92 if (rand < sum)
93 {
94 return GetHaplotypes((*genotype).second.first[which]);
95 }
96 }
97 //Just in case
98 return GetHaplotypes(*(*genotype).second.first.rbegin());
99 }
100
GetHaplotypes(string name)101 Haplotypes Phenotypes::GetHaplotypes(string name)
102 {
103 if (m_hapsmap.size() == 0)
104 {
105 MakeHaplotypes();
106 }
107 map<string, Haplotypes>::iterator hap = m_hapsmap.find(name);
108 if (hap == m_hapsmap.end())
109 {
110 throw data_error("Unable to find a set of haplotypes for the phenotype named " + name + ".");
111 }
112 return (*hap).second;
113 }
114
GetPhenotypesXML(long nspaces) const115 StringVec1d Phenotypes::GetPhenotypesXML(long nspaces) const
116 {
117 StringVec1d xmllines;
118 if (m_phenomap.size() == 0)
119 {
120 return xmllines;
121 }
122 string line = MakeIndent(MakeTag(xmlstr::XML_TAG_PHENOTYPES), nspaces);
123 xmllines.push_back(line);
124 nspaces += INDENT_DEPTH;
125
126 for (PhenMapConstIter genotype=m_phenomap.begin(); genotype != m_phenomap.end(); genotype++)
127 {
128 line = MakeIndent(MakeTag(xmlstr::XML_TAG_GENOTYPE), nspaces);
129 xmllines.push_back(line);
130 nspaces += INDENT_DEPTH;
131
132 line = MakeTag(xmlstr::XML_TAG_ALLELES) + " "
133 + ToString(MultisetElemToString((*genotype).first))
134 + MakeCloseTag(xmlstr::XML_TAG_ALLELES);
135 xmllines.push_back(MakeIndent(line, nspaces));
136
137 for (size_t phenotype=0; phenotype<(*genotype).second.second.size(); phenotype++)
138 {
139 line = MakeIndent(MakeTag(xmlstr::XML_TAG_PHENOTYPE), nspaces);
140 xmllines.push_back(line);
141 nspaces += INDENT_DEPTH;
142
143 line = MakeTag(xmlstr::XML_TAG_PHENOTYPE_NAME) + " "
144 + (*genotype).second.first[phenotype] + " "
145 + MakeCloseTag(xmlstr::XML_TAG_PHENOTYPE_NAME);
146 xmllines.push_back(MakeIndent(line, nspaces));
147
148 line = MakeTag(xmlstr::XML_TAG_PENETRANCE) + " "
149 + ToString((*genotype).second.second[phenotype]) + " "
150 + MakeCloseTag(xmlstr::XML_TAG_PENETRANCE);
151 xmllines.push_back(MakeIndent(line, nspaces));
152
153 nspaces -= INDENT_DEPTH;
154 line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_PHENOTYPE), nspaces);
155 xmllines.push_back(line);
156 }
157
158 nspaces -= INDENT_DEPTH;
159 line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_GENOTYPE), nspaces);
160 xmllines.push_back(line);
161 }
162 nspaces -= INDENT_DEPTH;
163 line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_PHENOTYPES), nspaces);
164 xmllines.push_back(line);
165 return xmllines;
166 }
167
MakeHaplotypes()168 void Phenotypes::MakeHaplotypes()
169 {
170 for(PhenMapIter genotype=m_phenomap.begin(); genotype != m_phenomap.end(); genotype++)
171 {
172 ScaleToSumToOne((*genotype).second.second);
173 assert((*genotype).second.first.size() == (*genotype).second.second.size());
174 for (size_t phenNum=0; phenNum<(*genotype).second.second.size(); phenNum++)
175 {
176 multiset<string> alleles = (*genotype).first;
177 string phenName = (*genotype).second.first[phenNum];
178 double penetrance = (*genotype).second.second[phenNum];
179
180 map<string, Haplotypes>::iterator hap = m_hapsmap.find(phenName);
181 if (hap == m_hapsmap.end())
182 {
183 Haplotypes newhaps(m_regionnum, m_locusname);
184 m_hapsmap.insert(make_pair(phenName, newhaps));
185 hap = m_hapsmap.find(phenName);
186 assert(hap != m_hapsmap.end());
187 }
188 (*hap).second.AddHaplotype(alleles, penetrance);
189 }
190 }
191 }
192
VecToSet(StringVec1d vec)193 multiset<string> Phenotypes::VecToSet(StringVec1d vec)
194 {
195 multiset<string> retset;
196 for (size_t str=0; str<vec.size(); str++)
197 {
198 retset.insert(vec[str]);
199 }
200 return retset;
201 }
202
203 //____________________________________________________________________________________
204