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