+++ /dev/null
-/*
- * pds.seqdist.cpp
- *
- *
- * Created by Pat Schloss on 8/12/11.
- * Copyright 2011 Patrick D. Schloss. All rights reserved.
- *
- */
-
-#include "myseqdist.h"
-#include "sequence.hpp"
-
-/**************************************************************************************************/
-correctDist::correctDist(int p) : processors(p) {
- try {
- m = MothurOut::getInstance();
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "correctDist");
- exit(1);
- }
-}
-/**************************************************************************************************/
-correctDist::correctDist(string sequenceFileName, int p) : processors(p) {
- try {
- m = MothurOut::getInstance();
- getSequences(sequenceFileName);
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "correctDist");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int correctDist::addSeq(string seqName, string seqSeq){
- try {
- names.push_back(seqName);
- sequences.push_back(fixSequence(seqSeq));
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "addSeq");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int correctDist::execute(string distanceFileName){
- try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
-#else
- processors = 1;
-#endif
- correctMatrix.resize(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];
- }
- }
-
- numSeqs = names.size();
-
- if(processors == 1){ driver(0, numSeqs, distanceFileName); }
- else{
-
- for(int i=0;i<processors;i++){
- start.push_back(int (sqrt(float(i)/float(processors)) * numSeqs));
- end.push_back(int (sqrt(float(i+1)/float(processors)) * numSeqs));
- }
-
- createProcess(distanceFileName);
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "execute");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int correctDist::getSequences(string sequenceFileName){
- try {
- ifstream sequenceFile;
- m->openInputFile(sequenceFileName, sequenceFile);
- string seqName, seqSeq;
-
- while(!sequenceFile.eof()){
- if (m->control_pressed) { break; }
-
- Sequence temp(sequenceFile); m->gobble(sequenceFile);
-
- if (temp.getName() != "") {
- names.push_back(temp.getName());
- sequences.push_back(fixSequence(temp.getAligned()));
- }
- }
- sequenceFile.close();
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "getSequences");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-vector<int> correctDist::fixSequence(string sequence){
- try {
- int alignLength = sequence.length();
-
- vector<int> seqVector;
-
- for(int i=0;i<alignLength;i++){
- if(sequence[i] == 'A') { seqVector.push_back(0); }
- else if(sequence[i] == 'C') { seqVector.push_back(1); }
- else if(sequence[i] == 'T') { seqVector.push_back(2); }
- else if(sequence[i] == 'G') { seqVector.push_back(3); }
- else if(sequence[i] == 'N') { seqVector.push_back(0); }//hmmmm....
- }
-
- return seqVector;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "fixSequence");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int correctDist::createProcess(string distanceFileName){
- try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
- int process = 1;
- vector<int> processIDs;
-
- while(process != processors){
-
- int pid = fork();
-
- if(pid > 0){
- processIDs.push_back(pid);
- process++;
- }
- else if(pid == 0){
- driver(start[process], end[process], distanceFileName + toString(getpid()) + ".temp");
- exit(0);
- }
- else{
- cout << "Doh!" << endl;
- for (int i=0;i<processIDs.size();i++) { int temp = processIDs[i]; kill(temp, SIGINT); }
- exit(0);
- }
- }
-
- driver(start[0], end[0], distanceFileName);
-
- for (int i=0;i<processIDs.size();i++) {
- int temp = processIDs[i];
- wait(&temp);
- }
-
- for(int i=0;i<processIDs.size();i++){
- m->appendFiles((distanceFileName + toString(processIDs[i]) + ".temp"), distanceFileName);
- remove((distanceFileName + toString(processIDs[i]) + ".temp").c_str());
- }
-#endif
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "createProcess");
- exit(1);
- }
-}
-
-/**************************************************************************************************/
-
-int correctDist::driver(int start, int end, string distFileName){
- try {
- ofstream distFile;
- m->openOutputFile(distFileName, distFile);
- distFile << setprecision(9);
-
- if(start == 0){
- distFile << numSeqs << endl;
- }
-
- int startTime = time(NULL);
-
- m->mothurOut("\nCalculating distances...\n");
-
- for(int i=start;i<end;i++){
-
- distFile << i;
-
- for(int j=0;j<i;j++){
-
- if (m->control_pressed) { distFile.close(); return 0; }
-
- double dist = getDist(sequences[i], sequences[j]);
-
- distFile << ' ' << dist;
- }
- distFile << endl;
-
- if(i % 100 == 0){ m->mothurOut(toString(i) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine(); }
- }
- distFile.close();
-
- if((end-1) % 100 != 0){ m->mothurOut(toString(end-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine(); }
- m->mothurOut("Done.\n");
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "driver");
- exit(1);
- }
-}
-/**************************************************************************************************/
-double correctDist::getDist(vector<int>& seqA, vector<int>& seqB){
- try {
-
- int lengthA = seqA.size();
- int lengthB = seqB.size();
-
- vector<vector<double> > alignMatrix(lengthA+1);
- vector<vector<char> > alignMoves(lengthA+1);
-
- for(int i=0;i<=lengthA;i++){
- alignMatrix[i].resize(lengthB+1, 0);
- alignMoves[i].resize(lengthB+1, 'x');
- }
-
- for(int i=0;i<=lengthA;i++){
- alignMatrix[i][0] = 15.0 * i;
- alignMoves[i][0] = 'u';
- }
- for(int i=0;i<=lengthB;i++){
- alignMatrix[0][i] = 15.0 * i;
- alignMoves[0][i] = 'l';
- }
-
- for(int i=1;i<=lengthA;i++){
- for(int j=1;j<=lengthB;j++){
-
- if (m->control_pressed) { return 0; }
-
- double nogap;
- nogap = alignMatrix[i-1][j-1] + correctMatrix[seqA[i-1]][seqB[j-1]];
-
-
- double gap;
- double left;
- if(i == lengthA){ //terminal gap
- left = alignMatrix[i][j-1];
- }
- else{
- if(seqB[j-1] == getLastMatch('l', alignMoves, i, j, seqA, seqB)){
- gap = 4.0;
- }
- else{
- gap = 15.0;
- }
-
- left = alignMatrix[i][j-1] + gap;
- }
-
-
- double up;
- if(j == lengthB){ //terminal gap
- up = alignMatrix[i-1][j];
- }
- else{
-
- if(seqA[i-1] == getLastMatch('u', alignMoves, i, j, seqA, seqB)){
- 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 = lengthA;
- int j = lengthB;
- int count = 0;
-
-
- // string alignA = "";
- // string alignB = "";
- // string bases = "ACTG";
- //
- // for(int i=0;i<lengthA;i++){
- // cout << bases[seqA[i]];
- // }cout << endl;
- //
- // for(int i=0;i<lengthB;i++){
- // cout << bases[seqB[i]];
- // }cout << endl;
-
- while(i > 0 && j > 0){
-
- if (m->control_pressed) { return 0; }
-
- if(alignMoves[i][j] == 'd'){
- // alignA = bases[seqA[i-1]] + alignA;
- // alignB = bases[seqB[j-1]] + alignB;
-
- count++;
- i--;
- j--;
- }
- else if(alignMoves[i][j] == 'u'){
- if(j != lengthB){
- // alignA = bases[seqA[i-1]] + alignA;
- // alignB = '-' + alignB;
- count++;
- }
-
- i--;
- }
- else if(alignMoves[i][j] == 'l'){
- if(i != lengthA){
- // alignA = '-' + alignA;
- // alignB = bases[seqB[j-1]] + alignB;
- count++;
- }
-
- j--;
- }
- }
-
- // cout << alignA << endl << alignB << endl;
-
- return alignMatrix[lengthA][lengthB] / (double)count;
- }
- catch(exception& e) {
- m->errorOut(e, "correctDist", "getDist");
- exit(1);
- }
-}
-/**************************************************************************************************/
-int correctDist::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, "correctDist", "getLastMatch");
- exit(1);
- }
-}
-/**************************************************************************************************/
-
-
-
projects / mothur.git / blobdiff