5 * Created by Pat Schloss on 8/31/11.
6 * Copyright 2011 Patrick D. Schloss. All rights reserved.
11 #include "sequence.hpp"
13 #define MIN_DELTA 1.0e-6
17 #define MIN_TAU 1.0e-4
18 #define MIN_WEIGHT 0.1
21 /**************************************************************************************************/
22 int seqNoise::getSequenceData(string sequenceFileName, vector<string>& sequences){
25 ifstream sequenceFile;
26 m->openInputFile(sequenceFileName, sequenceFile);
28 while(!sequenceFile.eof()){
30 if (m->control_pressed) { break; }
32 Sequence temp(sequenceFile); m->gobble(sequenceFile);
34 if (temp.getName() != "") {
35 sequences.push_back(temp.getAligned());
43 m->errorOut(e, "seqNoise", "getSequenceData");
47 /**************************************************************************************************/
48 int seqNoise::addSeq(string seq, vector<string>& sequences){
50 sequences.push_back(seq);
54 m->errorOut(e, "seqNoise", "addSeq");
58 /**************************************************************************************************/
59 //no checks for file mismatches
60 int seqNoise::getRedundantNames(string namesFileName, vector<string>& uniqueNames, vector<string>& redundantNames, vector<int>& seqFreq){
62 string unique, redundant;
64 m->openInputFile(namesFileName, namesFile);
66 for(int i=0;i<redundantNames.size();i++){
68 if (m->control_pressed) { break; }
70 namesFile >> uniqueNames[i]; m->gobble(namesFile);
71 namesFile >> redundantNames[i]; m->gobble(namesFile);
73 seqFreq[i] = m->getNumNames(redundantNames[i]);
80 m->errorOut(e, "seqNoise", "getRedundantNames");
84 /**************************************************************************************************/
85 int seqNoise::addRedundantName(string uniqueName, string redundantName, vector<string>& uniqueNames, vector<string>& redundantNames, vector<int>& seqFreq){
88 uniqueNames.push_back(uniqueName);
89 redundantNames.push_back(redundantName);
90 seqFreq.push_back(m->getNumNames(redundantName));
95 m->errorOut(e, "seqNoise", "addRedundantName");
99 /**************************************************************************************************/
100 int seqNoise::getDistanceData(string distFileName, vector<double>& distances){
104 m->openInputFile(distFileName, distFile);
111 for(int i=0;i<numSeqs;i++){
113 if (m->control_pressed) { break; }
115 distances[i * numSeqs + i] = 0.0000;
119 for(int j=0;j<i;j++){
120 distFile >> distances[i * numSeqs + j];
121 distances[j * numSeqs + i] = distances[i * numSeqs + j];
129 catch(exception& e) {
130 m->errorOut(e, "seqNoise", "getDistanceData");
135 /**************************************************************************************************/
136 int seqNoise::getListData(string listFileName, double cutOff, vector<int>& otuData, vector<int>& otuFreq, vector<vector<int> >& otuBySeqLookUp){
140 m->openInputFile(listFileName, listFile);
144 if(listFile.peek() == 'u'){ m->getline(listFile); }
147 listFile >> threshold;
149 if(threshold < cutOff){
150 m->getline(listFile);
154 otuFreq.resize(numOTUs, 0);
156 for(int i=0;i<numOTUs;i++){
158 if (m->control_pressed) { return 0; }
167 for(int j=0;j<otu.size();j++){
172 int index = atoi(number.c_str());
179 int index = atoi(number.c_str());
186 otuBySeqLookUp.resize(numOTUs);
188 int numSeqs = otuData.size();
190 for(int i=0;i<numSeqs;i++){
191 if (m->control_pressed) { return 0; }
192 otuBySeqLookUp[otuData[i]].push_back(i);
194 for(int i=0;i<numOTUs;i++){
195 if (m->control_pressed) { return 0; }
196 for(int j=otuBySeqLookUp[i].size();j<numSeqs;j++){
197 otuBySeqLookUp[i].push_back(0);
206 catch(exception& e) {
207 m->errorOut(e, "seqNoise", "getListData");
212 /**************************************************************************************************/
213 int seqNoise::updateOTUCountData(vector<int> otuFreq,
214 vector<vector<int> > otuBySeqLookUp,
215 vector<vector<int> > aanI,
218 vector<int>& cumCount
221 int numOTUs = otuFreq.size();
225 for(int i=0;i<numOTUs;i++){
228 if (m->control_pressed) { return 0; }
230 for(int j=0;j<otuFreq[i];j++){
231 anP[count] = otuBySeqLookUp[i][j];
232 anI[count] = aanI[i][j];
240 catch(exception& e) {
241 m->errorOut(e, "seqNoise", "updateOTUCountData");
245 /**************************************************************************************************/
246 double seqNoise::calcNewWeights(
247 vector<double>& weights, //
248 vector<int> seqFreq, //
250 vector<int> cumCount, //
252 vector<int> otuFreq, //
253 vector<double> tau //
257 int numOTUs = weights.size();
258 double maxChange = -1;
262 for(int i=0;i<numOTUs;i++){
264 if (m->control_pressed) { return 0; }
266 double change = weights[i];
270 for(int j=0;j<otuFreq[i];j++){
272 int index1 = cumCount[i] + j;
273 int index2 = anI[index1];
275 double currentTau = tau[anP[index1]];
276 double freq = double(seqFreq[index2]);
278 weights[i] += currentTau * freq;
280 change = fabs(weights[i] - change);
282 if(change > maxChange){ maxChange = change; }
288 catch(exception& e) {
289 m->errorOut(e, "seqNoise", "calcNewWeights");
294 /**************************************************************************************************/
296 int seqNoise::calcCentroids(
300 vector<int>& centroids,
301 vector<int> cumCount,
302 vector<double> distances,///
308 int numOTUs = change.size();
309 int numSeqs = seqFreq.size();
311 for(int i=0;i<numOTUs;i++){
313 if (m->control_pressed) { return 0; }
316 double minFValue = 1e10;
319 double count = 0.00000;
321 int freqOfOTU = otuFreq[i];
323 for(int j=0;j<freqOfOTU;j++){
324 int index = cumCount[i] + j;
325 count += seqFreq[anI[index]]*tau[anP[index]];
328 if(freqOfOTU > 0 && count > MIN_COUNT){
330 vector<double> adF(freqOfOTU);
331 vector<int> anL(freqOfOTU);
333 for(int j=0;j<freqOfOTU;j++){
334 anL[j] = anI[cumCount[i] + j];
338 for(int j=0;j<freqOfOTU;j++){
339 int index = cumCount[i] + j;
340 double curTau = tau[anP[index]];
342 for(int k=0;k<freqOfOTU;k++){
343 double dist = distances[anL[j]*numSeqs + anL[k]];
345 adF[k] += dist * curTau * seqFreq[anL[j]];
349 for(int j=0;j<freqOfOTU;j++){
350 if(adF[j] < minFValue){
356 if(centroids[i] != anL[minFIndex]){
358 centroids[i] = anL[minFIndex];
361 else if(centroids[i] != -1){
369 catch(exception& e) {
370 m->errorOut(e, "seqNoise", "calcCentroids");
375 /**************************************************************************************************/
377 int seqNoise::checkCentroids(vector<double>& weights, vector<int> centroids){
379 int numOTUs = centroids.size();
380 vector<int> unique(numOTUs, 1);
382 double minWeight = MIN_WEIGHT;
383 for(int i=0;i<numOTUs;i++){
384 if (m->control_pressed) { return 0; }
385 if(weights[i] < minWeight){ unique[i] = -1; }
388 for(int i=0;i<numOTUs;i++){
389 if (m->control_pressed) { return 0; }
391 for(int j=i+1; j<numOTUs;j++){
393 if(centroids[i] == centroids[j]){
395 weights[i] += weights[j];
405 catch(exception& e) {
406 m->errorOut(e, "seqNoise", "checkCentroids");
411 /**************************************************************************************************/
413 int seqNoise::setUpOTUData(vector<int>& otuData, vector<double>& percentage, vector<int> cumCount, vector<double> tau, vector<int> otuFreq, vector<int> anP, vector<int> anI){
416 int numOTUs = cumCount.size();
417 int numSeqs = otuData.size();
419 vector<double> bestTau(numSeqs, 0);
420 vector<double> bestIndex(numSeqs, -1);
422 for(int i=0;i<numOTUs;i++){
423 if (m->control_pressed) { return 0; }
424 for(int j=0;j<otuFreq[i];j++){
426 int index1 = cumCount[i] + j;
427 double thisTau = tau[anP[index1]];
428 int index2 = anI[index1];
430 if(thisTau > bestTau[index2]){
431 bestTau[index2] = thisTau;
432 bestIndex[index2] = i;
437 for(int i=0;i<numSeqs;i++){
438 if (m->control_pressed) { return 0; }
439 otuData[i] = bestIndex[i];
440 percentage[i] = 1 - bestTau[i];
444 catch(exception& e) {
445 m->errorOut(e, "seqNoise", "setUpOTUData");
450 /**************************************************************************************************/
452 int seqNoise::finishOTUData(vector<int> otuData, vector<int>& otuFreq, vector<int>& anP, vector<int>& anI, vector<int>& cumCount, vector<vector<int> >& otuBySeqLookUp, vector<vector<int> >& aanI, vector<double>& tau){
454 int numSeqs = otuData.size();
455 int numOTUs = otuFreq.size();
458 otuFreq.assign(numOTUs, 0);
459 tau.assign(numSeqs, 1);
460 anP.assign(numSeqs, 0);
461 anI.assign(numSeqs, 0);
463 for(int i=0;i<numSeqs;i++){
464 if (m->control_pressed) { return 0; }
465 int otu = otuData[i];
468 otuBySeqLookUp[otu][otuFreq[otu]] = i;
469 aanI[otu][otuFreq[otu]] = i;
472 updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
475 catch(exception& e) {
476 m->errorOut(e, "seqNoise", "finishOTUData");
481 /**************************************************************************************************/
483 int seqNoise::getLastMatch(char direction, vector<vector<char> >& alignMoves, int i, int j, vector<int>& seqA, vector<int>& seqB){
485 char nullReturn = -1;
488 if (m->control_pressed) { return nullReturn; }
489 if(direction == 'd'){
490 if(seqA[i-1] == seqB[j-1]) { return seqA[i-1]; }
491 else { return nullReturn; }
494 else if(direction == 'l') { j--; }
497 direction = alignMoves[i][j];
502 catch(exception& e) {
503 m->errorOut(e, "seqNoise", "getLastMatch");
508 /**************************************************************************************************/
510 int seqNoise::countDiffs(vector<int> query, vector<int> ref){
512 //double MATCH = 5.0;
513 //double MISMATCH = -2.0;
516 vector<vector<double> > correctMatrix(4);
517 for(int i=0;i<4;i++){ correctMatrix[i].resize(4); }
519 correctMatrix[0][0] = 0.000000; //AA
520 correctMatrix[1][0] = 11.619259; //CA
521 correctMatrix[2][0] = 11.694004; //TA
522 correctMatrix[3][0] = 7.748623; //GA
524 correctMatrix[1][1] = 0.000000; //CC
525 correctMatrix[2][1] = 7.619657; //TC
526 correctMatrix[3][1] = 12.852562; //GC
528 correctMatrix[2][2] = 0.000000; //TT
529 correctMatrix[3][2] = 10.964048; //TG
531 correctMatrix[3][3] = 0.000000; //GG
533 for(int i=0;i<4;i++){
534 for(int j=0;j<i;j++){
535 correctMatrix[j][i] = correctMatrix[i][j];
539 int queryLength = query.size();
540 int refLength = ref.size();
542 vector<vector<double> > alignMatrix(queryLength + 1);
543 vector<vector<char> > alignMoves(queryLength + 1);
545 for(int i=0;i<=queryLength;i++){
546 if (m->control_pressed) { return 0; }
547 alignMatrix[i].resize(refLength + 1, 0);
548 alignMoves[i].resize(refLength + 1, 'x');
551 for(int i=0;i<=queryLength;i++){
552 if (m->control_pressed) { return 0; }
553 alignMatrix[i][0] = 15.0 * i;
554 alignMoves[i][0] = 'u';
557 for(int i=0;i<=refLength;i++){
558 if (m->control_pressed) { return 0; }
559 alignMatrix[0][i] = 15.0 * i;
560 alignMoves[0][i] = 'l';
563 for(int i=1;i<=queryLength;i++){
564 if (m->control_pressed) { return 0; }
565 for(int j=1;j<=refLength;j++){
568 nogap = alignMatrix[i-1][j-1] + correctMatrix[query[i-1]][ref[j-1]];
573 if(i == queryLength){ //terminal gap
574 left = alignMatrix[i][j-1];
577 if(ref[j-1] == getLastMatch('l', alignMoves, i, j, query, ref)){
584 left = alignMatrix[i][j-1] + gap;
589 if(j == refLength){ //terminal gap
590 up = alignMatrix[i-1][j];
594 if(query[i-1] == getLastMatch('u', alignMoves, i, j, query, ref)){
601 up = alignMatrix[i-1][j] + gap;
608 alignMoves[i][j] = 'd';
609 alignMatrix[i][j] = nogap;
612 alignMoves[i][j] = 'u';
613 alignMatrix[i][j] = up;
618 alignMoves[i][j] = 'l';
619 alignMatrix[i][j] = left;
622 alignMoves[i][j] = 'u';
623 alignMatrix[i][j] = up;
633 // string alignA = "";
634 // string alignB = "";
635 // string bases = "ACTG";
637 while(i > 0 && j > 0){
638 if (m->control_pressed) { return 0; }
639 if(alignMoves[i][j] == 'd'){
640 // alignA = bases[query[i-1]] + alignA;
641 // alignB = bases[ref[j-1]] + alignB;
643 if(query[i-1] != ref[j-1]) { diffs++; }
648 else if(alignMoves[i][j] == 'u'){
650 // alignA = bases[query[i-1]] + alignA;
651 // alignB = '-' + alignB;
658 else if(alignMoves[i][j] == 'l'){
659 if(i != queryLength){
660 // alignA = '-' + alignA;
661 // alignB = bases[ref[j-1]] + alignB;
670 // cout << diffs << endl;
671 // cout << alignA << endl;
672 // cout << alignB << endl;
677 catch(exception& e) {
678 m->errorOut(e, "seqNoise", "countDiffs");
684 /**************************************************************************************************/
686 vector<int> seqNoise::convertSeq(string bases){
688 vector<int> numbers(bases.length(), -1);
690 for(int i=0;i<bases.length();i++){
691 if (m->control_pressed) { return numbers; }
695 if(b == 'A') { numbers[i] = 0; }
696 else if(b=='C') { numbers[i] = 1; }
697 else if(b=='T') { numbers[i] = 2; }
698 else if(b=='G') { numbers[i] = 3; }
699 else { numbers[i] = 0; }
704 catch(exception& e) {
705 m->errorOut(e, "seqNoise", "convertSeq");
710 /**************************************************************************************************/
712 string seqNoise::degapSeq(string aligned){
714 string unaligned = "";
716 for(int i=0;i<aligned.length();i++){
718 if (m->control_pressed) { return ""; }
720 if(aligned[i] != '-' && aligned[i] != '.'){
721 unaligned += aligned[i];
727 catch(exception& e) {
728 m->errorOut(e, "seqNoise", "degapSeq");
733 /**************************************************************************************************/
735 int seqNoise::writeOutput(string fastaFileName, string namesFileName, string uMapFileName, vector<int> finalTau, vector<int> centroids, vector<int> otuData, vector<string> sequences, vector<string> uniqueNames, vector<string> redundantNames, vector<int> seqFreq, vector<double>& distances){
737 int numOTUs = finalTau.size();
738 int numSeqs = uniqueNames.size();
740 ofstream fastaFile(fastaFileName.c_str());
741 ofstream namesFile(namesFileName.c_str());
742 ofstream uMapFile(uMapFileName.c_str());
744 vector<int> maxSequenceAbund(numOTUs, 0);
745 vector<int> maxSequenceIndex(numOTUs, 0);
747 for(int i=0;i<numSeqs;i++){
748 if (m->control_pressed) { return 0; }
749 if(maxSequenceAbund[otuData[i]] < seqFreq[i]){
750 maxSequenceAbund[otuData[i]] = seqFreq[i];
751 maxSequenceIndex[otuData[i]] = i;
757 for(int i=0;i<numOTUs;i++){
758 if (m->control_pressed) { return 0; }
762 if(maxSequenceIndex[i] != centroids[i] && distances[maxSequenceIndex[i]*numSeqs + centroids[i]] == 0){
763 // cout << uniqueNames[centroids[i]] << '\t' << uniqueNames[maxSequenceIndex[i]] << '\t' << count << endl;
764 centroids[i] = maxSequenceIndex[i];
767 int index = centroids[i];
769 fastaFile << '>' << uniqueNames[index] << endl << sequences[index] << endl;
770 namesFile << uniqueNames[index] << '\t';
772 string refSeq = sequences[index];
773 string redundantSeqs = redundantNames[index];;
776 vector<freqData> frequencyData;
778 for(int j=0;j<numSeqs;j++){
779 if(otuData[j] == i && j != index){
780 frequencyData.push_back(freqData(j, seqFreq[j]));
783 sort(frequencyData.rbegin(), frequencyData.rend());
785 string refDegap = degapSeq(refSeq);
786 vector<int> rUnalign = convertSeq(refDegap);
788 uMapFile << "ideal_seq_" << count << '\t' << finalTau[i] << endl;
789 uMapFile << uniqueNames[index] << '\t' << seqFreq[index] << "\t0\t" << refDegap << endl;
792 for(int j=0;j<frequencyData.size();j++){
793 if (m->control_pressed) { return 0; }
794 redundantSeqs += ',' + redundantNames[frequencyData[j].index];
796 uMapFile << uniqueNames[frequencyData[j].index] << '\t' << seqFreq[frequencyData[j].index] << '\t';
798 string querySeq = sequences[frequencyData[j].index];
800 string queryDegap = degapSeq(querySeq);
801 vector<int> qUnalign = convertSeq(queryDegap);
803 int udiffs = countDiffs(qUnalign, rUnalign);
804 uMapFile << udiffs << '\t' << queryDegap << endl;
809 namesFile << redundantSeqs << endl;
819 catch(exception& e) {
820 m->errorOut(e, "seqNoise", "writeOutput");
825 /**************************************************************************************************
827 int main(int argc, char *argv[]){
830 sigma = atof(argv[5]);
832 double cutOff = 0.08;
835 double minDelta = 1e-6;
837 string sequenceFileName = argv[1];
838 string fileNameStub = sequenceFileName.substr(0,sequenceFileName.find_last_of('.')) + ".shhh";
840 vector<string> sequences;
841 getSequenceData(sequenceFileName, sequences);
843 int numSeqs = sequences.size();
845 vector<string> uniqueNames(numSeqs);
846 vector<string> redundantNames(numSeqs);
847 vector<int> seqFreq(numSeqs);
849 string namesFileName = argv[4];
850 getRedundantNames(namesFileName, uniqueNames, redundantNames, seqFreq);
852 string distFileName = argv[2];
853 vector<double> distances(numSeqs * numSeqs);
854 getDistanceData(distFileName, distances);
856 string listFileName = argv[3];
857 vector<int> otuData(numSeqs);
859 vector<vector<int> > otuBySeqLookUp;
861 getListData(listFileName, cutOff, otuData, otuFreq, otuBySeqLookUp);
863 int numOTUs = otuFreq.size();
865 vector<double> weights(numOTUs, 0);
866 vector<int> change(numOTUs, 1);
867 vector<int> centroids(numOTUs, -1);
868 vector<int> cumCount(numOTUs, 0);
870 vector<double> tau(numSeqs, 1);
871 vector<int> anP(numSeqs, 0);
872 vector<int> anI(numSeqs, 0);
873 vector<int> anN(numSeqs, 0);
874 vector<vector<int> > aanI = otuBySeqLookUp;
877 double maxDelta = 1e6;
879 while(numIters < minIter || ((maxDelta > minDelta) && (numIters < maxIter))){
881 updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
882 maxDelta = calcNewWeights(weights, seqFreq, anI, cumCount, anP, otuFreq, tau);
884 calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau);
885 checkCentroids(weights, centroids);
887 otuFreq.assign(numOTUs, 0);
891 for(int i=0;i<numSeqs;i++){
893 double norm = 0.0000;
894 double minWeight = MIN_WEIGHT;
895 vector<double> currentTau(numOTUs);
897 for(int j=0;j<numOTUs;j++){
898 if(weights[j] > minWeight && distances[i * numSeqs+centroids[j]] < offset){
899 offset = distances[i * numSeqs+centroids[j]];
903 for(int j=0;j<numOTUs;j++){
904 if(weights[j] > minWeight){
905 currentTau[j] = exp(sigma * (-distances[(i * numSeqs + centroids[j])] + offset)) * weights[j];
906 norm += currentTau[j];
909 currentTau[j] = 0.0000;
913 for(int j=0;j<numOTUs;j++){
914 currentTau[j] /= norm;
917 for(int j=0;j<numOTUs;j++){
919 if(currentTau[j] > MIN_TAU){
920 int oldTotal = total;
923 tau.resize(oldTotal+1);
924 tau[oldTotal] = currentTau[j];
925 otuBySeqLookUp[j][otuFreq[j]] = oldTotal;
926 aanI[j][otuFreq[j]] = i;
939 updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
941 vector<double> percentage(numSeqs);
942 setUpOTUData(otuData, percentage, cumCount, tau, otuFreq, anP, anI);
943 finishOTUData(otuData, otuFreq, anP, anI, cumCount, otuBySeqLookUp, aanI, tau);
945 change.assign(numOTUs, 1);
946 calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau);
949 vector<int> finalTau(numOTUs, 0);
950 for(int i=0;i<numSeqs;i++){
951 finalTau[otuData[i]] += int(seqFreq[i]);
954 writeOutput(fileNameStub, finalTau, centroids, otuData, sequences, uniqueNames, redundantNames, seqFreq, distances);
959 /**************************************************************************************************/