-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathDistAlleleCounter.h
149 lines (123 loc) · 4.13 KB
/
DistAlleleCounter.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
/*
* DistAlleleCounter
* Date: Aug-15-2012
* Author : Gabriel Renaud gabriel.reno@gmail.com
*
*/
#ifndef DistAlleleCounter_h
#define DistAlleleCounter_h
#include <iostream>
#include <ostream>
#include <math.h>
#include <limits>
#include <vector>
#include <sstream>
#include "libgab.h"
//#include <iostream>
//#include <fstream>
using namespace std;
class DistAlleleCounter{
private:
public:
unsigned int count[16];
//vector<string> dimers ;
static const string dimers[16];
void reinitializedCounters();
void addAllelePair(int indexDimer);
string headerForCount() const;
unsigned int getIdent() const;
unsigned int getMutations() const;
DistAlleleCounter(); //constructor
/* double lowerConf (const unsigned int shortBranch,const unsigned int commonBranch) const; */
/* double highConf (const unsigned int shortBranch,const unsigned int commonBranch) const; */
/* friend ostream & operator<<(ostream & os, const DistAlleleCounter & ct); */
string getHeader(string prefixToAdd="");
/* double distRefSam () const ; */
/* double distSamRef () const ; */
double dist(const string & dnaDistMode) const;
DistAlleleCounter & operator+= (const DistAlleleCounter & other){
for(int i=0;i<16;i++)
count[i] += other.count[i] ;
/* counterSame += other.counterSame; */
/* counterCommon += other.counterCommon; */
/* counterReference += other.counterReference; */
/* counterSample += other.counterSample; */
return *this;
}
DistAlleleCounter & operator-= (const DistAlleleCounter & other){
for(int i=0;i<16;i++)
count[i] -= other.count[i] ;
/* counterSame -= other.counterSame; */
/* counterCommon -= other.counterCommon; */
/* counterReference -= other.counterReference; */
/* counterSample -= other.counterSample; */
return *this;
}
string printWithLabel() const;
friend ostream& operator<<(ostream& os, const DistAlleleCounter & ct){
/* for(int i=0;i<15;i++) */
/* os<<ct.counter[i]<<"\t"; */
/* os<<ct.counter[15]<<"\t"; */
unsigned int mutations=0;
unsigned int ident =0;
for(int i=0;i<16;i++){
if( i == 0 || i == 5 || i == 10 || i == 15 ){
ident+= ct.count[i] ;
}else{
mutations+= ct.count[i];
}
}
for(int i=0;i<15;i++){
os<<ct.count[i]<<"\t";
}
os<<ct.count[15]<<"\t"<<ident<<"\t"<<mutations<<"\t"<<double(mutations) /double(mutations+ident);
/* os<<ct.counterSame<<"\t" */
/* <<ct.counterCommon<<"\t" */
/* <<ct.counterReference<<"\t" */
/* <<ct.counterSample<<"\t"; */
/* if( (ct.counterReference+ct.counterCommon) != 0) */
/* os<<float(ct.counterReference)/float(ct.counterReference+ct.counterCommon)<<"\t"<<ct.lowerConf(ct.counterReference,ct.counterCommon)<<"\t"<<ct.highConf(ct.counterReference,ct.counterCommon); */
/* else */
/* os<<"nan"<<"\t"<<"nan"<<"\t"<<"nan"; */
/* os<<"\t"; */
/* if( (ct.counterSample+ct.counterCommon) != 0) */
/* os<<float(ct.counterSample)/float(ct.counterSample+ct.counterCommon)<<"\t"<<ct.lowerConf(ct.counterSample,ct.counterCommon)<<"\t"<<ct.highConf(ct.counterSample,ct.counterCommon); */
/* else */
/* os<<"nan"<<"\t"<<"nan"<<"\t"<<"nan"; */
return os;
}
//os<<ct.count[15]<<"\t"<<ident<<"\t"<<mutations<<"\t"<<double(mutations) /double(mutations+ident);
friend istream &operator>>(istream &is , DistAlleleCounter &ct){
//cerr<<"DstatCounter in"<<endl;
//double dst;
for(int i=0;i<=15;i++){
is>>ct.count[i];
}
double ident; double mutations;
is>>ident>>mutations;
double dist;//<<"\t"<<double(mutations) /double(mutations+ident);
string distStr;
char c;
is.get(c);
if(c != '\t'){
cerr<<"A single dist should have 19 fields"<<endl;
}
while (is.get(c)){ // loop getting single characters
if(c == '\t') break;
distStr+=c;
//cerr<<"#"<<c<<"#"<<endl;
}
if(distStr == "-nan"){
dist = -1.0*numeric_limits<double>::quiet_NaN();
}else{
dist = destringify<double>(distStr);
}
/* //cerr<<"#"; */
/* for(int i=0;i<=15;i++){ */
/* cerr<<ct.count[i]<<" "; */
/* } */
/* cerr<<"i:"<<ident<<" m:"<<mutations<<" d:"<<dist<<"# "; */
return is;
}
};
#endif