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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
|
//
// opticluster.cpp
// Mothur
//
// Created by Sarah Westcott on 4/20/16.
// Copyright (c) 2016 Schloss Lab. All rights reserved.
//
#include "opticluster.h"
OptiCluster::OptiCluster(OptiData* mt, ClusterMetric* met, long long ns) : Cluster() {
matrix = mt; metric = met; truePositives = 0; trueNegatives = 0; falseNegatives = 0; falsePositives = 0; numSingletons = ns;
}
/***********************************************************************/
//randomly assign sequences to OTUs
int OptiCluster::initialize(double& value, bool randomize, string initialize) {
try {
numSeqs = matrix->getNumSeqs();
truePositives = 0;
falsePositives = 0;
falseNegatives = 0;
trueNegatives = 0;
bins.resize(numSeqs); //place seqs in own bin
vector<long long> temp;
bins.push_back(temp);
seqBin[numSeqs] = -1;
insertLocation = numSeqs;
Utils util;
if (initialize == "singleton") {
//put everyone in own bin
for (int i = 0; i < numSeqs; i++) { bins[i].push_back(i); }
//maps randomized sequences to bins
for (int i = 0; i < numSeqs; i++) {
seqBin[i] = bins[i][0];
randomizeSeqs.push_back(i);
}
if (randomize) { util.mothurRandomShuffle(randomizeSeqs); }
//for each sequence (singletons removed on read)
for (map<long long, long long>::iterator it = seqBin.begin(); it != seqBin.end(); it++) {
if (it->second == -1) { }
else {
long long numCloseSeqs = (matrix->getNumClose(it->first)); //does not include self
falseNegatives += numCloseSeqs;
}
}
falseNegatives /= 2; //square matrix
trueNegatives = numSeqs * (numSeqs-1)/2 - (falsePositives + falseNegatives + truePositives); //since everyone is a singleton no one clusters together. True negative = num far apart
}else {
//put everyone in first bin
for (int i = 0; i < numSeqs; i++) {
bins[0].push_back(i);
seqBin[i] = 0;
randomizeSeqs.push_back(i);
}
if (randomize) { util.mothurRandomShuffle(randomizeSeqs); }
//for each sequence (singletons removed on read)
for (map<long long, long long>::iterator it = seqBin.begin(); it != seqBin.end(); it++) {
if (it->second == -1) { }
else {
long long numCloseSeqs = (matrix->getNumClose(it->first)); //does not include self
truePositives += numCloseSeqs;
}
}
truePositives /= 2; //square matrix
falsePositives = numSeqs * (numSeqs-1)/2 - (trueNegatives + falseNegatives + truePositives);
}
value = metric->getValue(truePositives, trueNegatives, falsePositives, falseNegatives);
return value;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "initialize");
exit(1);
}
}
/***********************************************************************/
/* for each sequence with mutual information (close)
* remove from current OTU and calculate MCC when sequence forms its own OTU or joins one of the other OTUs where there is a sequence within the `threshold` (no need to calculate MCC if the paired sequence is already in same OTU and no need to try every OTU - just those where there's a close sequence)
* keep or move the sequence to the OTU where the `metric` is the largest - flip a coin on ties */
bool OptiCluster::update(double& listMetric) {
try {
//for each sequence (singletons removed on read)
for (int i = 0; i < randomizeSeqs.size(); i++) {
if (m->getControl_pressed()) { break; }
map<long long, long long>::iterator it = seqBin.find(randomizeSeqs[i]);
long long seqNumber = it->first;
long long binNumber = it->second;
if (binNumber == -1) { }
else {
double tn, tp, fp, fn;
double bestMetric = -1;
double bestBin, bestTp, bestTn, bestFn, bestFp;
tn = trueNegatives; tp = truePositives; fp = falsePositives; fn = falseNegatives;
//close / far count in current bin
vector<double> results = getCloseFarCounts(seqNumber, binNumber);
double cCount = results[0]; double fCount = results[1];
//metric in current bin
bestMetric = metric->getValue(tp, tn, fp, fn); bestBin = binNumber; bestTp = tp; bestTn = tn; bestFp = fp; bestFn = fn;
//if not already singleton, then calc value if singleton was created
if (!((bins[binNumber].size()) == 1)) {
//make a singleton
//move out of old bin
fn+=cCount; tn+=fCount; fp-=fCount; tp-=cCount;
double singleMetric = metric->getValue(tp, tn, fp, fn);
if (singleMetric > bestMetric) {
bestBin = -1; bestTp = tp; bestTn = tn; bestFp = fp; bestFn = fn;
bestMetric = singleMetric;
}
}
set<long long> binsToTry;
set<long long> closeSeqs = matrix->getCloseSeqs(seqNumber);
for (set<long long>::iterator itClose = closeSeqs.begin(); itClose != closeSeqs.end(); itClose++) { binsToTry.insert(seqBin[*itClose]); }
//merge into each "close" otu
for (set<long long>::iterator it = binsToTry.begin(); it != binsToTry.end(); it++) {
tn = trueNegatives; tp = truePositives; fp = falsePositives; fn = falseNegatives;
fn+=cCount; tn+=fCount; fp-=fCount; tp-=cCount; //move out of old bin
results = getCloseFarCounts(seqNumber, *it);
fn-=results[0]; tn-=results[1]; tp+=results[0]; fp+=results[1]; //move into new bin
double newMetric = metric->getValue(tp, tn, fp, fn); //score when sequence is moved
//new best
if (newMetric > bestMetric) { bestMetric = newMetric; bestBin = (*it); bestTp = tp; bestTn = tn; bestFp = fp; bestFn = fn; }
}
bool usedInsert = false;
if (bestBin == -1) { bestBin = insertLocation; usedInsert = true; }
if (bestBin != binNumber) {
truePositives = bestTp; trueNegatives = bestTn; falsePositives = bestFp; falseNegatives = bestFn;
//move seq from i to j
bins[bestBin].push_back(seqNumber); //add seq to bestbin
bins[binNumber].erase(remove(bins[binNumber].begin(), bins[binNumber].end(), seqNumber), bins[binNumber].end()); //remove from old bin i
}
if (usedInsert) { insertLocation = findInsert(); }
//update seqBins
seqBin[seqNumber] = bestBin; //set new OTU location
}
}
listMetric = metric->getValue(truePositives, trueNegatives, falsePositives, falseNegatives);
if (m->getDebug()) { ListVector* list = getList(); list->print(cout); delete list; }
return 0;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "update");
exit(1);
}
}
/***********************************************************************/
vector<double> OptiCluster::getCloseFarCounts(long long seq, long long newBin) {
try {
vector<double> results; results.push_back(0); results.push_back(0);
if (newBin == -1) { } //making a singleton bin. Close but we are forcing apart.
else { //merging a bin
for (int i = 0; i < bins[newBin].size(); i++) {
if (seq == bins[newBin][i]) {} //ignore self
else if (!matrix->isClose(seq, bins[newBin][i])) { results[1]++; } //this sequence is "far away" from sequence i - above the cutoff
else { results[0]++; } //this sequence is "close" to sequence i - distance between them is less than cutoff
}
}
return results;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "getCloseFarCounts");
exit(1);
}
}
/***********************************************************************/
vector<double> OptiCluster::getStats( double& tp, double& tn, double& fp, double& fn) {
try {
double singletn = matrix->getNumSingletons() + numSingletons;
double tempnumSeqs = numSeqs + singletn;
tp = truePositives;
fp = falsePositives;
fn = falseNegatives;
tn = tempnumSeqs * (tempnumSeqs-1)/2 - (falsePositives + falseNegatives + truePositives); //adds singletons to tn
vector<double> results;
Sensitivity sens; double sensitivity = sens.getValue(tp, tn, fp, fn); results.push_back(sensitivity);
Specificity spec; double specificity = spec.getValue(tp, tn, fp, fn); results.push_back(specificity);
PPV ppv; double positivePredictiveValue = ppv.getValue(tp, tn, fp, fn); results.push_back(positivePredictiveValue);
NPV npv; double negativePredictiveValue = npv.getValue(tp, tn, fp, fn); results.push_back(negativePredictiveValue);
FDR fdr; double falseDiscoveryRate = fdr.getValue(tp, tn, fp, fn); results.push_back(falseDiscoveryRate);
Accuracy acc; double accuracy = acc.getValue(tp, tn, fp, fn); results.push_back(accuracy);
MCC mcc; double matthewsCorrCoef = mcc.getValue(tp, tn, fp, fn); results.push_back(matthewsCorrCoef);
F1Score f1; double f1Score = f1.getValue(tp, tn, fp, fn); results.push_back(f1Score);
return results;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "getStats");
exit(1);
}
}
/***********************************************************************/
ListVector* OptiCluster::getList() {
try {
ListVector* list = new ListVector();
ListVector* singleton = matrix->getListSingle();
if (singleton != nullptr) { //add in any sequences above cutoff in read. Removing these saves clustering time.
for (int i = 0; i < singleton->getNumBins(); i++) {
if (singleton->get(i) != "") {
list->push_back(singleton->get(i));
}
}
delete singleton;
}
for (int i = 0; i < bins.size(); i++) {
if (bins[i].size() != 0) {
string otu = matrix->getName(bins[i][0]);
for (int j = 1; j < bins[i].size(); j++) {
otu += "," + matrix->getName(bins[i][j]);
}
list->push_back(otu);
}
}
return list;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "getList");
exit(1);
}
}
/***********************************************************************/
long long OptiCluster::getNumBins() {
try {
long long singletn = matrix->getNumSingletons();
for (int i = 0; i < bins.size(); i++) {
if (bins[i].size() != 0) {
singletn++;
}
}
return singletn;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "getNumBins");
exit(1);
}
}
/***********************************************************************/
long long OptiCluster::findInsert() {
try {
//initially there are bins for each sequence (excluding singletons removed on read)
for (long long i = 0; i < bins.size(); i++) {
if (m->getControl_pressed()) { break; }
if (bins[i].size() == 0) { return i; } //this bin is empty
}
return -1;
}
catch(exception& e) {
m->errorOut(e, "OptiCluster", "findInsert");
exit(1);
}
}
/***********************************************************************/
|
