+++ /dev/null
-/*
- * mySeqNoise.cpp
- *
- *
- * Created by Pat Schloss on 8/31/11.
- * Copyright 2011 Patrick D. Schloss. All rights reserved.
- *
- */
-
-#include "seqnoise.h"
-#include "sequence.hpp"
-
-#define MIN_DELTA 1.0e-6
-#define MIN_ITER 20
-#define MAX_ITER 1000
-#define MIN_COUNT 0.1
-#define MIN_TAU 1.0e-4
-#define MIN_WEIGHT 0.1
-
-
-/**************************************************************************************************/
-int seqNoise::getSequenceData(string sequenceFileName, vector<string>& sequences){
- try {
-
- ifstream sequenceFile;
- m->openInputFile(sequenceFileName, sequenceFile);
-
- while(!sequenceFile.eof()){
-
- if (m->control_pressed) { break; }
-
- Sequence temp(sequenceFile); m->gobble(sequenceFile);
-
- if (temp.getName() != "") {
- sequences.push_back(temp.getAligned());
- }
- }
- sequenceFile.close();
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "getSequenceData");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int seqNoise::addSeq(string seq, vector<string>& sequences){
- try {
- sequences.push_back(seq);
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "addSeq");
- exit(1);
- }
-}
-/**************************************************************************************************/
-//no checks for file mismatches
-int seqNoise::getRedundantNames(string namesFileName, vector<string>& uniqueNames, vector<string>& redundantNames, vector<int>& seqFreq){
- try {
- string unique, redundant;
- ifstream namesFile;
- m->openInputFile(namesFileName, namesFile);
-
- for(int i=0;i<redundantNames.size();i++){
-
- if (m->control_pressed) { break; }
-
- namesFile >> uniqueNames[i]; m->gobble(namesFile);
- namesFile >> redundantNames[i]; m->gobble(namesFile);
-
- seqFreq[i] = m->getNumNames(redundantNames[i]);
- }
- namesFile.close();
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "getRedundantNames");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int seqNoise::addRedundantName(string uniqueName, string redundantName, vector<string>& uniqueNames, vector<string>& redundantNames, vector<int>& seqFreq){
- try {
-
- uniqueNames.push_back(uniqueName);
- redundantNames.push_back(redundantName);
- seqFreq.push_back(m->getNumNames(redundantName));
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "addRedundantName");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int seqNoise::getDistanceData(string distFileName, vector<double>& distances){
- try {
-
- ifstream distFile;
- m->openInputFile(distFileName, distFile);
-
- int numSeqs = 0;
- string name = "";
-
- distFile >> numSeqs;
-
- for(int i=0;i<numSeqs;i++){
-
- if (m->control_pressed) { break; }
-
- distances[i * numSeqs + i] = 0.0000;
-
- distFile >> name;
-
- for(int j=0;j<i;j++){
- distFile >> distances[i * numSeqs + j];
- distances[j * numSeqs + i] = distances[i * numSeqs + j];
- }
- }
-
- distFile.close();
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "getDistanceData");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-int seqNoise::getListData(string listFileName, double cutOff, vector<int>& otuData, vector<int>& otuFreq, vector<vector<int> >& otuBySeqLookUp){
- try {
-
- ifstream listFile;
- m->openInputFile(listFileName, listFile);
- double threshold;
- int numOTUs;
- string line = "";
- bool adjustCutoff = true;
-
- if(listFile.peek() == 'u'){ m->getline(listFile); }
-
- while(listFile){
- listFile >> threshold;
-
- if(threshold < cutOff){
- line = m->getline(listFile); m->gobble(listFile);
- }
- else{
- adjustCutoff = false;
- listFile >> numOTUs;
- otuFreq.resize(numOTUs, 0);
-
- for(int i=0;i<numOTUs;i++){
-
- if (m->control_pressed) { return 0; }
-
- string otu;
- listFile >> otu;
-
- int count = 0;
-
- string number = "";
-
- for(int j=0;j<otu.size();j++){
- if(otu[j] != ','){
- number += otu[j];
- }
- else{
- int index = atoi(number.c_str());
- otuData[index] = i;
- count++;
- number = "";
- }
- }
-
- int index = atoi(number.c_str());
- otuData[index] = i;
- count++;
-
- otuFreq[i] = count;
- }
-
- otuBySeqLookUp.resize(numOTUs);
-
- int numSeqs = otuData.size();
-
- for(int i=0;i<numSeqs;i++){
- if (m->control_pressed) { return 0; }
- otuBySeqLookUp[otuData[i]].push_back(i);
- }
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
- for(int j=otuBySeqLookUp[i].size();j<numSeqs;j++){
- otuBySeqLookUp[i].push_back(0);
- }
- }
-
- break;
- }
- }
-
- listFile.close();
-
- //the listfile does not contain a threshold greater than the cutoff so use highest value
- if (adjustCutoff) {
- istringstream iss (line,istringstream::in);
-
- iss >> numOTUs;
- otuFreq.resize(numOTUs, 0);
-
- for(int i=0;i<numOTUs;i++){
-
- if (m->control_pressed) { return 0; }
-
- string otu;
- iss >> otu;
-
- int count = 0;
-
- string number = "";
-
- for(int j=0;j<otu.size();j++){
- if(otu[j] != ','){
- number += otu[j];
- }
- else{
- int index = atoi(number.c_str());
- otuData[index] = i;
- count++;
- number = "";
- }
- }
-
- int index = atoi(number.c_str());
- otuData[index] = i;
- count++;
-
- otuFreq[i] = count;
- }
-
- otuBySeqLookUp.resize(numOTUs);
-
- int numSeqs = otuData.size();
-
- for(int i=0;i<numSeqs;i++){
- if (m->control_pressed) { return 0; }
- otuBySeqLookUp[otuData[i]].push_back(i);
- }
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
- for(int j=otuBySeqLookUp[i].size();j<numSeqs;j++){
- otuBySeqLookUp[i].push_back(0);
- }
- }
-
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "getListData");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-int seqNoise::updateOTUCountData(vector<int> otuFreq,
- vector<vector<int> > otuBySeqLookUp,
- vector<vector<int> > aanI,
- vector<int>& anP,
- vector<int>& anI,
- vector<int>& cumCount
- ){
- try {
- int numOTUs = otuFreq.size();
-
- int count = 0;
-
- for(int i=0;i<numOTUs;i++){
- cumCount[i] = count;
-
- if (m->control_pressed) { return 0; }
-
- for(int j=0;j<otuFreq[i];j++){
- anP[count] = otuBySeqLookUp[i][j];
- anI[count] = aanI[i][j];
-
- count++;
- }
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "updateOTUCountData");
- exit(1);
- }
-}
-/**************************************************************************************************/
-double seqNoise::calcNewWeights(
- vector<double>& weights, //
- vector<int> seqFreq, //
- vector<int> anI, //
- vector<int> cumCount, //
- vector<int> anP, //
- vector<int> otuFreq, //
- vector<double> tau //
- ){
- try {
-
- int numOTUs = weights.size();
- double maxChange = -1;
-
- cout.flush();
-
- for(int i=0;i<numOTUs;i++){
-
- if (m->control_pressed) { return 0; }
-
- double change = weights[i];
-
- weights[i] = 0.0000;
-
- for(int j=0;j<otuFreq[i];j++){
-
- int index1 = cumCount[i] + j;
- int index2 = anI[index1];
-
- double currentTau = tau[anP[index1]];
- double freq = double(seqFreq[index2]);
-
- weights[i] += currentTau * freq;
- }
- change = fabs(weights[i] - change);
-
- if(change > maxChange){ maxChange = change; }
- cout.flush();
- }
- return maxChange;
-
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "calcNewWeights");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int seqNoise::calcCentroids(
- vector<int> anI,
- vector<int> anP,
- vector<int>& change,
- vector<int>& centroids,
- vector<int> cumCount,
- vector<double> distances,///
- vector<int> seqFreq,
- vector<int> otuFreq,
- vector<double> tau
- ){
- try {
- int numOTUs = change.size();
- int numSeqs = seqFreq.size();
-
- for(int i=0;i<numOTUs;i++){
-
- if (m->control_pressed) { return 0; }
-
- int minFIndex = -1;
- double minFValue = 1e10;
-
- change[i] = 0;
- double count = 0.00000;
-
- int freqOfOTU = otuFreq[i];
-
- for(int j=0;j<freqOfOTU;j++){
- int index = cumCount[i] + j;
- count += seqFreq[anI[index]]*tau[anP[index]];
- }
-
- if(freqOfOTU > 0 && count > MIN_COUNT){
-
- vector<double> adF(freqOfOTU);
- vector<int> anL(freqOfOTU);
-
- for(int j=0;j<freqOfOTU;j++){
- anL[j] = anI[cumCount[i] + j];
- adF[j] = 0.0000;
- }
-
- for(int j=0;j<freqOfOTU;j++){
- int index = cumCount[i] + j;
- double curTau = tau[anP[index]];
-
- for(int k=0;k<freqOfOTU;k++){
- double dist = distances[anL[j]*numSeqs + anL[k]];
-
- adF[k] += dist * curTau * seqFreq[anL[j]];
- }
- }
-
- for(int j=0;j<freqOfOTU;j++){
- if(adF[j] < minFValue){
- minFIndex = j;
- minFValue = adF[j];
- }
- }
-
- if(centroids[i] != anL[minFIndex]){
- change[i] = 1;
- centroids[i] = anL[minFIndex];
- }
- }
- else if(centroids[i] != -1){
- change[i] = 1;
- centroids[i] = -1;
- }
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "calcCentroids");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int seqNoise::checkCentroids(vector<double>& weights, vector<int> centroids){
- try {
- int numOTUs = centroids.size();
- vector<int> unique(numOTUs, 1);
-
- double minWeight = MIN_WEIGHT;
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
- if(weights[i] < minWeight){ unique[i] = -1; }
- }
-
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
- if(unique[i] == 1){
- for(int j=i+1; j<numOTUs;j++){
- if(unique[j] == 1){
- if(centroids[i] == centroids[j]){
- unique[j] = 0;
- weights[i] += weights[j];
- weights[j] = 0;
- }
- }
- }
- }
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "checkCentroids");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int seqNoise::setUpOTUData(vector<int>& otuData, vector<double>& percentage, vector<int> cumCount, vector<double> tau, vector<int> otuFreq, vector<int> anP, vector<int> anI){
- try {
-
- int numOTUs = cumCount.size();
- int numSeqs = otuData.size();
-
- vector<double> bestTau(numSeqs, 0);
- vector<double> bestIndex(numSeqs, -1);
-
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
- for(int j=0;j<otuFreq[i];j++){
-
- int index1 = cumCount[i] + j;
- double thisTau = tau[anP[index1]];
- int index2 = anI[index1];
-
- if(thisTau > bestTau[index2]){
- bestTau[index2] = thisTau;
- bestIndex[index2] = i;
- }
- }
- }
-
- for(int i=0;i<numSeqs;i++){
- if (m->control_pressed) { return 0; }
- otuData[i] = bestIndex[i];
- percentage[i] = 1 - bestTau[i];
- }
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "setUpOTUData");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-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){
- try {
- int numSeqs = otuData.size();
- int numOTUs = otuFreq.size();
- int total = numSeqs;
-
- otuFreq.assign(numOTUs, 0);
- tau.assign(numSeqs, 1);
- anP.assign(numSeqs, 0);
- anI.assign(numSeqs, 0);
-
- for(int i=0;i<numSeqs;i++){
- if (m->control_pressed) { return 0; }
- int otu = otuData[i];
- total++;
-
- otuBySeqLookUp[otu][otuFreq[otu]] = i;
- aanI[otu][otuFreq[otu]] = i;
- otuFreq[otu]++;
- }
- updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "finishOTUData");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int seqNoise::getLastMatch(char direction, vector<vector<char> >& alignMoves, int i, int j, vector<int>& seqA, vector<int>& seqB){
- try{
- char nullReturn = -1;
-
- while(i>=1 && j>=1){
- if (m->control_pressed) { return nullReturn; }
- if(direction == 'd'){
- if(seqA[i-1] == seqB[j-1]) { return seqA[i-1]; }
- else { return nullReturn; }
- }
-
- else if(direction == 'l') { j--; }
- else { i--; }
-
- direction = alignMoves[i][j];
- }
-
- return nullReturn;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "getLastMatch");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int seqNoise::countDiffs(vector<int> query, vector<int> ref){
- try {
- //double MATCH = 5.0;
- //double MISMATCH = -2.0;
- //double GAP = -2.0;
-
- vector<vector<double> > correctMatrix(4);
- for(int i=0;i<4;i++){ correctMatrix[i].resize(4); }
-
- correctMatrix[0][0] = 0.000000; //AA
- correctMatrix[1][0] = 11.619259; //CA
- correctMatrix[2][0] = 11.694004; //TA
- correctMatrix[3][0] = 7.748623; //GA
-
- correctMatrix[1][1] = 0.000000; //CC
- correctMatrix[2][1] = 7.619657; //TC
- correctMatrix[3][1] = 12.852562; //GC
-
- correctMatrix[2][2] = 0.000000; //TT
- correctMatrix[3][2] = 10.964048; //TG
-
- correctMatrix[3][3] = 0.000000; //GG
-
- for(int i=0;i<4;i++){
- for(int j=0;j<i;j++){
- correctMatrix[j][i] = correctMatrix[i][j];
- }
- }
-
- int queryLength = query.size();
- int refLength = ref.size();
-
- vector<vector<double> > alignMatrix(queryLength + 1);
- vector<vector<char> > alignMoves(queryLength + 1);
-
- for(int i=0;i<=queryLength;i++){
- if (m->control_pressed) { return 0; }
- alignMatrix[i].resize(refLength + 1, 0);
- alignMoves[i].resize(refLength + 1, 'x');
- }
-
- for(int i=0;i<=queryLength;i++){
- if (m->control_pressed) { return 0; }
- alignMatrix[i][0] = 15.0 * i;
- alignMoves[i][0] = 'u';
- }
-
- for(int i=0;i<=refLength;i++){
- if (m->control_pressed) { return 0; }
- alignMatrix[0][i] = 15.0 * i;
- alignMoves[0][i] = 'l';
- }
-
- for(int i=1;i<=queryLength;i++){
- if (m->control_pressed) { return 0; }
- for(int j=1;j<=refLength;j++){
-
- double nogap;
- nogap = alignMatrix[i-1][j-1] + correctMatrix[query[i-1]][ref[j-1]];
-
-
- double gap;
- double left;
- if(i == queryLength){ //terminal gap
- left = alignMatrix[i][j-1];
- }
- else{
- if(ref[j-1] == getLastMatch('l', alignMoves, i, j, query, ref)){
- gap = 4.0;
- }
- else{
- gap = 15.0;
- }
-
- left = alignMatrix[i][j-1] + gap;
- }
-
-
- double up;
- if(j == refLength){ //terminal gap
- up = alignMatrix[i-1][j];
- }
- else{
-
- if(query[i-1] == getLastMatch('u', alignMoves, i, j, query, ref)){
- gap = 4.0;
- }
- else{
- gap = 15.0;
- }
-
- up = alignMatrix[i-1][j] + gap;
- }
-
-
-
- if(nogap < left){
- if(nogap < up){
- alignMoves[i][j] = 'd';
- alignMatrix[i][j] = nogap;
- }
- else{
- alignMoves[i][j] = 'u';
- alignMatrix[i][j] = up;
- }
- }
- else{
- if(left < up){
- alignMoves[i][j] = 'l';
- alignMatrix[i][j] = left;
- }
- else{
- alignMoves[i][j] = 'u';
- alignMatrix[i][j] = up;
- }
- }
- }
- }
-
- int i = queryLength;
- int j = refLength;
- int diffs = 0;
-
- // string alignA = "";
- // string alignB = "";
- // string bases = "ACTG";
-
- while(i > 0 && j > 0){
- if (m->control_pressed) { return 0; }
- if(alignMoves[i][j] == 'd'){
- // alignA = bases[query[i-1]] + alignA;
- // alignB = bases[ref[j-1]] + alignB;
-
- if(query[i-1] != ref[j-1]) { diffs++; }
-
- i--;
- j--;
- }
- else if(alignMoves[i][j] == 'u'){
- if(j != refLength){
- // alignA = bases[query[i-1]] + alignA;
- // alignB = '-' + alignB;
-
- diffs++;
- }
-
- i--;
- }
- else if(alignMoves[i][j] == 'l'){
- if(i != queryLength){
- // alignA = '-' + alignA;
- // alignB = bases[ref[j-1]] + alignB;
-
- diffs++;
- }
-
- j--;
- }
- }
-
- // cout << diffs << endl;
- // cout << alignA << endl;
- // cout << alignB << endl;
- // cout << endl;
-
- return diffs;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "countDiffs");
- exit(1);
- }
-
-}
-
-/**************************************************************************************************/
-
-vector<int> seqNoise::convertSeq(string bases){
- try {
- vector<int> numbers(bases.length(), -1);
-
- for(int i=0;i<bases.length();i++){
- if (m->control_pressed) { return numbers; }
-
- char b = bases[i];
-
- if(b == 'A') { numbers[i] = 0; }
- else if(b=='C') { numbers[i] = 1; }
- else if(b=='T') { numbers[i] = 2; }
- else if(b=='G') { numbers[i] = 3; }
- else { numbers[i] = 0; }
- }
-
- return numbers;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "convertSeq");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-string seqNoise::degapSeq(string aligned){
- try {
- string unaligned = "";
-
- for(int i=0;i<aligned.length();i++){
-
- if (m->control_pressed) { return ""; }
-
- if(aligned[i] != '-' && aligned[i] != '.'){
- unaligned += aligned[i];
- }
- }
-
- return unaligned;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "degapSeq");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-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){
- try {
- int numOTUs = finalTau.size();
- int numSeqs = uniqueNames.size();
-
- ofstream fastaFile(fastaFileName.c_str());
- ofstream namesFile(namesFileName.c_str());
- ofstream uMapFile(uMapFileName.c_str());
-
- vector<int> maxSequenceAbund(numOTUs, 0);
- vector<int> maxSequenceIndex(numOTUs, 0);
-
- for(int i=0;i<numSeqs;i++){
- if (m->control_pressed) { return 0; }
- if(maxSequenceAbund[otuData[i]] < seqFreq[i]){
- maxSequenceAbund[otuData[i]] = seqFreq[i];
- maxSequenceIndex[otuData[i]] = i;
- }
- }
-
- int count = 1;
-
- for(int i=0;i<numOTUs;i++){
- if (m->control_pressed) { return 0; }
-
- if(finalTau[i] > 0){
-
- if(maxSequenceIndex[i] != centroids[i] && distances[maxSequenceIndex[i]*numSeqs + centroids[i]] == 0){
- // cout << uniqueNames[centroids[i]] << '\t' << uniqueNames[maxSequenceIndex[i]] << '\t' << count << endl;
- centroids[i] = maxSequenceIndex[i];
- }
-
- int index = centroids[i];
-
- fastaFile << '>' << uniqueNames[index] << endl << sequences[index] << endl;
- namesFile << uniqueNames[index] << '\t';
-
- string refSeq = sequences[index];
- string redundantSeqs = redundantNames[index];;
-
-
- vector<freqData> frequencyData;
-
- for(int j=0;j<numSeqs;j++){
- if(otuData[j] == i && j != index){
- frequencyData.push_back(freqData(j, seqFreq[j]));
- }
- }
- sort(frequencyData.rbegin(), frequencyData.rend());
-
- string refDegap = degapSeq(refSeq);
- vector<int> rUnalign = convertSeq(refDegap);
-
- uMapFile << "ideal_seq_" << count << '\t' << finalTau[i] << endl;
- uMapFile << uniqueNames[index] << '\t' << seqFreq[index] << "\t0\t" << refDegap << endl;
-
-
- for(int j=0;j<frequencyData.size();j++){
- if (m->control_pressed) { return 0; }
- redundantSeqs += ',' + redundantNames[frequencyData[j].index];
-
- uMapFile << uniqueNames[frequencyData[j].index] << '\t' << seqFreq[frequencyData[j].index] << '\t';
-
- string querySeq = sequences[frequencyData[j].index];
-
- string queryDegap = degapSeq(querySeq);
- vector<int> qUnalign = convertSeq(queryDegap);
-
- int udiffs = countDiffs(qUnalign, rUnalign);
- uMapFile << udiffs << '\t' << queryDegap << endl;
-
- }
-
- uMapFile << endl;
- namesFile << redundantSeqs << endl;
- count++;
-
- }
- }
- fastaFile.close();
- namesFile.close();
- uMapFile.close();
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "seqNoise", "writeOutput");
- exit(1);
- }
-}
-
-/**************************************************************************************************
-
-int main(int argc, char *argv[]){
-
- double sigma = 100;
- sigma = atof(argv[5]);
-
- double cutOff = 0.08;
- int minIter = 10;
- int maxIter = 1000;
- double minDelta = 1e-6;
-
- string sequenceFileName = argv[1];
- string fileNameStub = sequenceFileName.substr(0,sequenceFileName.find_last_of('.')) + ".shhh";
-
- vector<string> sequences;
- getSequenceData(sequenceFileName, sequences);
-
- int numSeqs = sequences.size();
-
- vector<string> uniqueNames(numSeqs);
- vector<string> redundantNames(numSeqs);
- vector<int> seqFreq(numSeqs);
-
- string namesFileName = argv[4];
- getRedundantNames(namesFileName, uniqueNames, redundantNames, seqFreq);
-
- string distFileName = argv[2];
- vector<double> distances(numSeqs * numSeqs);
- getDistanceData(distFileName, distances);
-
- string listFileName = argv[3];
- vector<int> otuData(numSeqs);
- vector<int> otuFreq;
- vector<vector<int> > otuBySeqLookUp;
-
- getListData(listFileName, cutOff, otuData, otuFreq, otuBySeqLookUp);
-
- int numOTUs = otuFreq.size();
-
- vector<double> weights(numOTUs, 0);
- vector<int> change(numOTUs, 1);
- vector<int> centroids(numOTUs, -1);
- vector<int> cumCount(numOTUs, 0);
-
- vector<double> tau(numSeqs, 1);
- vector<int> anP(numSeqs, 0);
- vector<int> anI(numSeqs, 0);
- vector<int> anN(numSeqs, 0);
- vector<vector<int> > aanI = otuBySeqLookUp;
-
- int numIters = 0;
- double maxDelta = 1e6;
-
- while(numIters < minIter || ((maxDelta > minDelta) && (numIters < maxIter))){
-
- updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
- maxDelta = calcNewWeights(weights, seqFreq, anI, cumCount, anP, otuFreq, tau);
-
- calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau);
- checkCentroids(weights, centroids);
-
- otuFreq.assign(numOTUs, 0);
-
- int total = 0;
-
- for(int i=0;i<numSeqs;i++){
- double offset = 1e6;
- double norm = 0.0000;
- double minWeight = MIN_WEIGHT;
- vector<double> currentTau(numOTUs);
-
- for(int j=0;j<numOTUs;j++){
- if(weights[j] > minWeight && distances[i * numSeqs+centroids[j]] < offset){
- offset = distances[i * numSeqs+centroids[j]];
- }
- }
-
- for(int j=0;j<numOTUs;j++){
- if(weights[j] > minWeight){
- currentTau[j] = exp(sigma * (-distances[(i * numSeqs + centroids[j])] + offset)) * weights[j];
- norm += currentTau[j];
- }
- else{
- currentTau[j] = 0.0000;
- }
- }
-
- for(int j=0;j<numOTUs;j++){
- currentTau[j] /= norm;
- }
-
- for(int j=0;j<numOTUs;j++){
-
- if(currentTau[j] > MIN_TAU){
- int oldTotal = total;
- total++;
-
- tau.resize(oldTotal+1);
- tau[oldTotal] = currentTau[j];
- otuBySeqLookUp[j][otuFreq[j]] = oldTotal;
- aanI[j][otuFreq[j]] = i;
- otuFreq[j]++;
-
- }
- }
-
- anP.resize(total);
- anI.resize(total);
- }
-
- numIters++;
- }
-
- updateOTUCountData(otuFreq, otuBySeqLookUp, aanI, anP, anI, cumCount);
-
- vector<double> percentage(numSeqs);
- setUpOTUData(otuData, percentage, cumCount, tau, otuFreq, anP, anI);
- finishOTUData(otuData, otuFreq, anP, anI, cumCount, otuBySeqLookUp, aanI, tau);
-
- change.assign(numOTUs, 1);
- calcCentroids(anI, anP, change, centroids, cumCount, distances, seqFreq, otuFreq, tau);
-
-
- vector<int> finalTau(numOTUs, 0);
- for(int i=0;i<numSeqs;i++){
- finalTau[otuData[i]] += int(seqFreq[i]);
- }
-
- writeOutput(fileNameStub, finalTau, centroids, otuData, sequences, uniqueNames, redundantNames, seqFreq, distances);
-
- return 0;
-}
-
-**************************************************************************************************/
projects / mothur.git / blobdiff