← Back to branch summary

~ubuntu-branches/debian/sid/libvcflib/sid

« back to all changes in this revision

Viewing changes to src/var.cpp

  • Committer: Package Import Robot
  • Author(s): Andreas Tille
  • Date: 2016-09-16 15:52:29 UTC
  • Revision ID: package-import@ubuntu.com-20160916155229-24mxrntfylvsshsg
Tags: upstream-1.0.0~rc1+dfsg
Import upstream version 1.0.0~rc1+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/var.cpp
 
1
 
 
2
#include "var.hpp"
 
3
 
 
4
genotype::~genotype(){}
 
5
 
 
6
zvar::~zvar(){}
 
7
 
 
8
void zvar::setPopName(string  popName){
 
9
  name = popName;
 
10
}
 
11
 
 
12
pooled::~pooled(){}
 
13
 
 
14
double pooled::bound(double v){
 
15
  if(v <= 0.00001){
 
16
    return  0.00001;
 
17
  }
 
18
  if(v >= 0.99999){
 
19
    return 0.99999;
 
20
  }
 
21
  return v;
 
22
}
 
23
 
 
24
gl::~gl(){}
 
25
 
 
26
gl::gl(void){
 
27
  nalt  = 0;
 
28
  nref  = 0;
 
29
  af    = 0;
 
30
  nhomr = 0;
 
31
  nhoma = 0;
 
32
  nhet  = 0;
 
33
  ngeno = 0;
 
34
  fis   = 0;
 
35
  hfrq  = 0;
 
36
 
 
37
  alpha = 0.01;
 
38
  beta  = 0.01;
 
39
}
 
40
 
 
41
pl::~pl(){}
 
42
 
 
43
pl::pl(void){
 
44
  nalt  = 0;
 
45
  nref  = 0;
 
46
  af    = 0;
 
47
  nhomr = 0;
 
48
  nhoma = 0;
 
49
  nhet  = 0;
 
50
  ngeno = 0;
 
51
  fis   = 0;
 
52
  hfrq  = 0;
 
53
 
 
54
  alpha = 0.01;
 
55
  beta  = 0.01;
 
56
 
 
57
}
 
58
 
 
59
gp::~gp(){}
 
60
 
 
61
gp::gp(void){
 
62
  nalt  = 0;
 
63
  nref  = 0;
 
64
  af    = 0;
 
65
  nhomr = 0;
 
66
  nhoma = 0;
 
67
  nhet  = 0;
 
68
  ngeno = 0;
 
69
  fis   = 0;
 
70
  hfrq  = 0;
 
71
 
 
72
  alpha = 0.01;
 
73
  beta  = 0.01;
 
74
 
 
75
}
 
76
 
 
77
gt::~gt(){}
 
78
 
 
79
gt::gt(void){
 
80
  nalt  = 0;
 
81
  nref  = 0;
 
82
  af    = 0;
 
83
  nhomr = 0;
 
84
  nhoma = 0;
 
85
  nhet  = 0;
 
86
  ngeno = 0;
 
87
  fis   = 0;
 
88
  hfrq  = 0;
 
89
 
 
90
  alpha = 0.01;
 
91
  beta  = 0.01;
 
92
 
 
93
}
 
94
 
 
95
pooled::pooled(void){
 
96
 
 
97
  npop   = 0;
 
98
  afsum  = 0;
 
99
  nalt   = 0;
 
100
  nref   = 0;
 
101
  af     = 0;
 
102
  npop   = 0;
 
103
  ntot   = 0;
 
104
 
 
105
  alpha  = 0.01;
 
106
  beta   = 0.01;
 
107
 
 
108
}
 
109
 
 
110
// polymorphism for GL and PL
 
111
 
 
112
double gt::unphred(map< string, vector<string> > & geno, int index){
 
113
  return -1;
 
114
}
 
115
 
 
116
double gl::unphred(map< string, vector<string> > & geno, int index){
 
117
 
 
118
  double unphreded = atof(geno["GL"][index].c_str());
 
119
  return unphreded;
 
120
}
 
121
 
 
122
double gp::unphred(map< string, vector<string> > & geno, int index){
 
123
 
 
124
  double unphreded = atof(geno["GP"][index].c_str());
 
125
  return log(unphreded) ;
 
126
}
 
127
 
 
128
double pl::unphred(map< string, vector<string> > & geno, int index){
 
129
 
 
130
   double unphreded = atof(geno["PL"][index].c_str());
 
131
 
 
132
   unphreded = log(pow(10, (-unphreded/10)));
 
133
     
 
134
   return unphreded;
 
135
}
 
136
 
 
137
void pooled::loadPop(vector< map< string, vector<string> > > & group, string seqid, long int position){
 
138
  vector< map< string, vector<string> > >::iterator targ_it = group.begin();
 
139
 
 
140
 
 
141
  for(; targ_it != group.end(); targ_it++){
 
142
 
 
143
    string genotype = (*targ_it)["GT"].front();
 
144
 
 
145
    if(genotype == "./."){
 
146
      continue;
 
147
    }
 
148
 
 
149
    string allelecounts = (*targ_it)["AD"].front();
 
150
 
 
151
    vector<string> ac   = (*targ_it)["AD"];
 
152
    
 
153
    npop += 1;
 
154
    
 
155
 
 
156
    double af = atof(ac[1].c_str()) / ( atof(ac[0].c_str()) + atof(ac[1].c_str()) );
 
157
 
 
158
    if(atof(ac[1].c_str()) == 0){
 
159
      af = 0;
 
160
    }
 
161
    if(atof(ac[0].c_str()) == 0){
 
162
      af = 1;
 
163
    }
 
164
        
 
165
    afsum += af;
 
166
                                 
 
167
    afs.push_back(af);
 
168
 
 
169
    nrefs.push_back(atof(ac[0].c_str()));
 
170
    nalts.push_back(atof(ac[1].c_str()));
 
171
 
 
172
    nref += atof(ac[0].c_str());
 
173
    ntot += atof(ac[0].c_str());
 
174
    nalt += atof(ac[1].c_str());
 
175
    ntot += atof(ac[1].c_str());
 
176
  }
 
177
  if(npop < 1){
 
178
    af = -1;
 
179
  }
 
180
  else{
 
181
    af = afsum / npop;
 
182
  }
 
183
}
 
184
 
 
185
void genotype::estimatePosterior(void){
 
186
 
 
187
  int ng = genoIndex.size();
 
188
 
 
189
  for(int i = 0 ; i < ng; i++){
 
190
 
 
191
    if(genoIndex[i] == -1){
 
192
      continue;
 
193
    }
 
194
 
 
195
    double aa = genoLikelihoods[i][0] ;
 
196
    double ab = genoLikelihoods[i][1] ;
 
197
    double bb = genoLikelihoods[i][2] ;
 
198
 
 
199
    alpha += exp(ab);
 
200
    beta  += exp(ab);
 
201
 
 
202
    alpha += 2 * exp(aa);
 
203
    beta  += 2 * exp(bb);
 
204
  }
 
205
  
 
206
  
 
207
}
 
208
 
 
209
void pooled::estimatePosterior(void){
 
210
  
 
211
  if(npop < 2){
 
212
    cerr << "FATAL: not enough pooled populations in the target or background\n";
 
213
    exit(1);
 
214
  }
 
215
  
 
216
  double xbar = af;
 
217
  double ss = 0;
 
218
  
 
219
  for(int i = 0 ; i < npop; i++){
 
220
    ss += pow(( afs[i] - xbar),2);
 
221
  }
 
222
  
 
223
  double var = (1/(npop-1))*ss;
 
224
  
 
225
  xbar = bound(xbar);
 
226
  if(var < 0.01){
 
227
    var = 0.01;
 
228
  }
 
229
  if(var < xbar*(1-xbar)){
 
230
    alpha = xbar * (((xbar*(1-xbar))/var) -1);
 
231
    beta  = (1 - xbar) * (((xbar*(1-xbar))/var) -1);
 
232
  }
 
233
  else{
 
234
    alpha = -1;
 
235
    beta  = -1;
 
236
  }
 
237
}
 
238
 
 
239
void genotype::loadPop( vector< map< string, vector<string> > >& group, string seqid, long int position){
 
240
 
 
241
  seqid = seqid;
 
242
  pos   = position  ;
 
243
 
 
244
  vector< map< string, vector<string> > >::iterator targ_it = group.begin();
 
245
 
 
246
  for(; targ_it != group.end(); targ_it++){
 
247
        
 
248
    string genotype = (*targ_it)["GT"].front();
 
249
 
 
250
    gts.push_back(genotype);
 
251
 
 
252
    vector<double> phreds;
 
253
    vector<double> phredsCDF;
 
254
 
 
255
    double sum  = 0;
 
256
    if(genotype != "./."){
 
257
 
 
258
      double pa  =  unphred( (*targ_it), 0)  ;
 
259
      double pab =  unphred( (*targ_it), 1)  ;
 
260
      double pbb =  unphred( (*targ_it), 2)  ;
 
261
 
 
262
      double norm = log(exp(pa) + exp(pab) + exp(pbb))  ;
 
263
 
 
264
      phreds.push_back(pa  - norm);
 
265
      phreds.push_back(pab - norm);
 
266
      phreds.push_back(pbb - norm);
 
267
 
 
268
      
 
269
       // cerr << exp(pa  - norm) << endl;
 
270
       // cerr << exp(pab - norm) << endl;
 
271
       // cerr << exp(pbb - norm) << endl;
 
272
       // cerr << endl;
 
273
 
 
274
      sum += exp(pa - norm);
 
275
      //      cerr << sum << endl;
 
276
      phredsCDF.push_back(sum);
 
277
      sum += exp(pab - norm);
 
278
      //      cerr << sum << endl;
 
279
      phredsCDF.push_back(sum);
 
280
      sum += exp(pbb - norm);
 
281
      //      cerr << sum << endl;
 
282
      //      cerr << endl;
 
283
      phredsCDF.push_back(sum);
 
284
    }
 
285
    else{
 
286
      phreds.push_back(log(1/3));
 
287
      phreds.push_back(log(1/3));
 
288
      phreds.push_back(log(1/3));
 
289
      phredsCDF.push_back(1/3);
 
290
      phredsCDF.push_back(2/3);
 
291
      phredsCDF.push_back(1);
 
292
    }
 
293
 
 
294
    genoLikelihoods.push_back(phreds);
 
295
    genoLikelihoodsCDF.push_back(phredsCDF);
 
296
 
 
297
    if(genotype == "./0" || genotype == "./1"){
 
298
      genotype = "./.";
 
299
    }
 
300
 
 
301
    while(1){
 
302
      if(genotype == "./."){
 
303
        genoIndex.push_back(-1);
 
304
        break;
 
305
      }
 
306
      if(genotype == "0/0"){
 
307
        ngeno += 1;
 
308
        nhomr += 1;
 
309
        nref  += 2;
 
310
        genoIndex.push_back(0);
 
311
        break;
 
312
      }
 
313
      if(genotype == "0/1"){
 
314
        ngeno += 1;
 
315
        nhet  += 1;
 
316
        nref  += 1;
 
317
        nalt  += 1;
 
318
        genoIndex.push_back(1);
 
319
        break;
 
320
      }
 
321
      if(genotype == "1/0"){
 
322
        ngeno += 1;
 
323
        nhet  += 1;
 
324
        nref  += 1;
 
325
        nalt  += 1;
 
326
        genoIndex.push_back(1);
 
327
        break;
 
328
      }
 
329
      if(genotype == "1/1"){
 
330
        ngeno += 1;
 
331
        nhoma += 1;
 
332
        nalt  += 2;
 
333
        genoIndex.push_back(2);
 
334
        break;
 
335
      }
 
336
      if(genotype == "0|0"){
 
337
        ngeno += 1;
 
338
        nhomr += 1;
 
339
        nref  += 2;
 
340
        genoIndex.push_back(0);
 
341
        break;
 
342
      }
 
343
      if(genotype == "0|1"){
 
344
        ngeno += 1;
 
345
        nhet  += 1;
 
346
        nref  += 1;
 
347
        nalt  += 1;
 
348
        genoIndex.push_back(1);
 
349
        break;
 
350
      }
 
351
      if(genotype == "1|0"){
 
352
        ngeno += 1;
 
353
        nhet  += 1;
 
354
        nref  += 1;
 
355
        nalt  += 1;
 
356
        genoIndex.push_back(1);
 
357
        break;
 
358
      }
 
359
      if(genotype == "1|1"){
 
360
        ngeno += 1;
 
361
        nhoma += 1;
 
362
        nalt  += 2;
 
363
        genoIndex.push_back(2);
 
364
        break;
 
365
      }
 
366
      cerr << "FATAL: unknown genotype: " << genotype << endl;
 
367
      exit(1);
 
368
    }
 
369
  }
 
370
  if(nalt == 0 && nref == 0){
 
371
    af = -1;
 
372
  }
 
373
  else{
 
374
    af  = (nalt / (nref + nalt));
 
375
 
 
376
    if(nhet > 0){
 
377
      fis = ( 1 - ((nhet/ngeno) / (2*af*(1 - af))));
 
378
    }
 
379
    else{
 
380
      fis = 1;
 
381
    }
 
382
    if(fis < 0){
 
383
      fis = 0.00001;
 
384
    }
 
385
  }
 
386
  hfrq = nhet / ngeno;
 
387
  npop = ngeno;
 
388
}