--- /dev/null
+/*
+ * readblast.cpp
+ * Mothur
+ *
+ * Created by westcott on 12/10/09.
+ * Copyright 2009 Schloss Lab. All rights reserved.
+ *
+ */
+
+#include "readblast.h"
+#include "progress.hpp"
+
+//********************************************************************************************************************
+//sorts lowest to highest
+inline bool compareOverlap(seqDist left, seqDist right){
+ return (left.dist < right.dist);
+}
+/*********************************************************************************************/
+ReadBlast::ReadBlast(string file, float c, float p, int l, bool ms, bool h) : blastfile(file), cutoff(c), penalty(p), length(l), minWanted(ms), hclusterWanted(h) {
+ try {
+ m = MothurOut::getInstance();
+ matrix = NULL;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ReadBlast", "ReadBlast");
+ exit(1);
+ }
+}
+/*********************************************************************************************/
+//assumptions about the blast file:
+//1. if duplicate lines occur the first line is always best and is chosen
+//2. blast scores are grouped together, ie. a a .... score, a b .... score, a c ....score...
+int ReadBlast::read(NameAssignment* nameMap) {
+ try {
+
+ //if the user has not given a names file read names from blastfile
+ if (nameMap->size() == 0) { readNames(nameMap); }
+ int nseqs = nameMap->size();
+
+ if (m->control_pressed) { return 0; }
+
+ ifstream fileHandle;
+ m->openInputFile(blastfile, fileHandle);
+
+ string firstName, secondName, eScore, currentRow;
+ string repeatName = "";
+ int count = 1;
+ float distance, thisoverlap, refScore;
+ float percentId;
+ float numBases, mismatch, gap, startQuery, endQuery, startRef, endRef, score, lengthThisSeq;
+
+ ofstream outDist;
+ ofstream outOverlap;
+
+ //create objects needed for read
+ if (!hclusterWanted) {
+ matrix = new SparseMatrix();
+ }else{
+ overlapFile = m->getRootName(blastfile) + "overlap.dist";
+ distFile = m->getRootName(blastfile) + "hclusterDists.dist";
+
+ m->openOutputFile(overlapFile, outOverlap);
+ m->openOutputFile(distFile, outDist);
+ }
+
+ if (m->control_pressed) {
+ fileHandle.close();
+ if (!hclusterWanted) { delete matrix; }
+ else { outOverlap.close(); m->mothurRemove(overlapFile); outDist.close(); m->mothurRemove(distFile); }
+ return 0;
+ }
+
+ Progress* reading = new Progress("Reading blast: ", nseqs * nseqs);
+
+ //this is used to quickly find if we already have a distance for this combo
+ vector< map<int,float> > dists; dists.resize(nseqs); //dists[0][1] = distance from seq0 to seq1
+ map<int, float> thisRowsBlastScores;
+
+ if (!fileHandle.eof()) {
+ //read in line from file
+ fileHandle >> firstName >> secondName >> percentId >> numBases >> mismatch >> gap >> startQuery >> endQuery >> startRef >> endRef >> eScore >> score;
+ m->gobble(fileHandle);
+
+ currentRow = firstName;
+ lengthThisSeq = numBases;
+ repeatName = firstName + secondName;
+
+ if (firstName == secondName) { refScore = score; }
+ else{
+ //convert name to number
+ map<string,int>::iterator itA = nameMap->find(firstName);
+ map<string,int>::iterator itB = nameMap->find(secondName);
+ if(itA == nameMap->end()){ m->mothurOut("AAError: Sequence '" + firstName + "' was not found in the names file, please correct\n"); exit(1); }
+ if(itB == nameMap->end()){ m->mothurOut("ABError: Sequence '" + secondName + "' was not found in the names file, please correct\n"); exit(1); }
+
+ thisRowsBlastScores[itB->second] = score;
+
+ //calc overlap score
+ thisoverlap = 1.0 - (percentId * (lengthThisSeq - startQuery) / endRef / 100.0 - penalty);
+
+ //if there is a valid overlap, add it
+ if ((startRef <= length) && ((endQuery+length) >= lengthThisSeq) && (thisoverlap < cutoff)) {
+ if (!hclusterWanted) {
+ seqDist overlapValue(itA->second, itB->second, thisoverlap);
+ overlap.push_back(overlapValue);
+ }else {
+ outOverlap << itA->first << '\t' << itB->first << '\t' << thisoverlap << endl;
+ }
+ }
+ }
+ }else { m->mothurOut("Error in your blast file, cannot read."); m->mothurOutEndLine(); exit(1); }
+
+
+ //read file
+ while(!fileHandle.eof()){
+
+ if (m->control_pressed) {
+ fileHandle.close();
+ if (!hclusterWanted) { delete matrix; }
+ else { outOverlap.close(); m->mothurRemove(overlapFile); outDist.close(); m->mothurRemove(distFile); }
+ delete reading;
+ return 0;
+ }
+
+ //read in line from file
+ fileHandle >> firstName >> secondName >> percentId >> numBases >> mismatch >> gap >> startQuery >> endQuery >> startRef >> endRef >> eScore >> score;
+ //cout << firstName << '\t' << secondName << '\t' << percentId << '\t' << numBases << '\t' << mismatch << '\t' << gap << '\t' << startQuery << '\t' << endQuery << '\t' << startRef << '\t' << endRef << '\t' << eScore << '\t' << score << endl;
+ m->gobble(fileHandle);
+
+ string temp = firstName + secondName; //to check if this file has repeat lines, ie. is this a blast instead of a blscreen file
+
+ //if this is a new pairing
+ if (temp != repeatName) {
+ repeatName = temp;
+
+ if (currentRow == firstName) {
+ //cout << "first = " << firstName << " second = " << secondName << endl;
+ if (firstName == secondName) {
+ refScore = score;
+ reading->update((count + nseqs));
+ count++;
+ }else{
+ //convert name to number
+ map<string,int>::iterator itA = nameMap->find(firstName);
+ map<string,int>::iterator itB = nameMap->find(secondName);
+ if(itA == nameMap->end()){ m->mothurOut("AAError: Sequence '" + firstName + "' was not found in the names file, please correct\n"); exit(1); }
+ if(itB == nameMap->end()){ m->mothurOut("ABError: Sequence '" + secondName + "' was not found in the names file, please correct\n"); exit(1); }
+
+ //save score
+ thisRowsBlastScores[itB->second] = score;
+
+ //calc overlap score
+ thisoverlap = 1.0 - (percentId * (lengthThisSeq - startQuery) / endRef / 100.0 - penalty);
+
+ //if there is a valid overlap, add it
+ if ((startRef <= length) && ((endQuery+length) >= lengthThisSeq) && (thisoverlap < cutoff)) {
+ if (!hclusterWanted) {
+ seqDist overlapValue(itA->second, itB->second, thisoverlap);
+ //cout << "overlap = " << itA->second << '\t' << itB->second << '\t' << thisoverlap << endl;
+ overlap.push_back(overlapValue);
+ }else {
+ outOverlap << itA->first << '\t' << itB->first << '\t' << thisoverlap << endl;
+ }
+ }
+ } //end else
+ }else { //end row
+ //convert blast scores to distance and add cell to sparse matrix if we can
+ map<int, float>::iterator it;
+ map<int, float>::iterator itDist;
+ for(it=thisRowsBlastScores.begin(); it!=thisRowsBlastScores.end(); it++) {
+ distance = 1.0 - (it->second / refScore);
+
+
+ //do we already have the distance calculated for b->a
+ map<string,int>::iterator itA = nameMap->find(currentRow);
+ itDist = dists[it->first].find(itA->second);
+
+ //if we have it then compare
+ if (itDist != dists[it->first].end()) {
+
+ //if you want the minimum blast score ratio, then pick max distance
+ if(minWanted) { distance = max(itDist->second, distance); }
+ else{ distance = min(itDist->second, distance); }
+
+ //is this distance below cutoff
+ if (distance < cutoff) {
+ if (!hclusterWanted) {
+ PCell value(itA->second, it->first, distance);
+ matrix->addCell(value);
+ }else{
+ outDist << itA->first << '\t' << nameMap->get(it->first) << '\t' << distance << endl;
+ }
+ }
+ //not going to need this again
+ dists[it->first].erase(itDist);
+ }else { //save this value until we get the other ratio
+ dists[itA->second][it->first] = distance;
+ }
+ }
+ //clear out last rows info
+ thisRowsBlastScores.clear();
+
+ currentRow = firstName;
+ lengthThisSeq = numBases;
+
+ //add this row to thisRowsBlastScores
+ if (firstName == secondName) { refScore = score; }
+ else{ //add this row to thisRowsBlastScores
+
+ //convert name to number
+ map<string,int>::iterator itA = nameMap->find(firstName);
+ map<string,int>::iterator itB = nameMap->find(secondName);
+ if(itA == nameMap->end()){ m->mothurOut("AAError: Sequence '" + firstName + "' was not found in the names file, please correct\n"); exit(1); }
+ if(itB == nameMap->end()){ m->mothurOut("ABError: Sequence '" + secondName + "' was not found in the names file, please correct\n"); exit(1); }
+
+ thisRowsBlastScores[itB->second] = score;
+
+ //calc overlap score
+ thisoverlap = 1.0 - (percentId * (lengthThisSeq - startQuery) / endRef / 100.0 - penalty);
+
+ //if there is a valid overlap, add it
+ if ((startRef <= length) && ((endQuery+length) >= lengthThisSeq) && (thisoverlap < cutoff)) {
+ if (!hclusterWanted) {
+ seqDist overlapValue(itA->second, itB->second, thisoverlap);
+ overlap.push_back(overlapValue);
+ }else {
+ outOverlap << itA->first << '\t' << itB->first << '\t' << thisoverlap << endl;
+ }
+ }
+ }
+ }//end if current row
+ }//end if repeat
+ }//end while
+
+ //get last rows info stored
+ //convert blast scores to distance and add cell to sparse matrix if we can
+ map<int, float>::iterator it;
+ map<int, float>::iterator itDist;
+ for(it=thisRowsBlastScores.begin(); it!=thisRowsBlastScores.end(); it++) {
+ distance = 1.0 - (it->second / refScore);
+
+ //do we already have the distance calculated for b->a
+ map<string,int>::iterator itA = nameMap->find(currentRow);
+ itDist = dists[it->first].find(itA->second);
+
+ //if we have it then compare
+ if (itDist != dists[it->first].end()) {
+ //if you want the minimum blast score ratio, then pick max distance
+ if(minWanted) { distance = max(itDist->second, distance); }
+ else{ distance = min(itDist->second, distance); }
+
+ //is this distance below cutoff
+ if (distance < cutoff) {
+ if (!hclusterWanted) {
+ PCell value(itA->second, it->first, distance);
+ matrix->addCell(value);
+ }else{
+ outDist << itA->first << '\t' << nameMap->get(it->first) << '\t' << distance << endl;
+ }
+ }
+ //not going to need this again
+ dists[it->first].erase(itDist);
+ }else { //save this value until we get the other ratio
+ dists[itA->second][it->first] = distance;
+ }
+ }
+ //clear out info
+ thisRowsBlastScores.clear();
+ dists.clear();
+
+ if (m->control_pressed) {
+ fileHandle.close();
+ if (!hclusterWanted) { delete matrix; }
+ else { outOverlap.close(); m->mothurRemove(overlapFile); outDist.close(); m->mothurRemove(distFile); }
+ delete reading;
+ return 0;
+ }
+
+ if (!hclusterWanted) {
+ sort(overlap.begin(), overlap.end(), compareOverlap);
+ }else {
+ outDist.close();
+ outOverlap.close();
+ }
+
+ if (m->control_pressed) {
+ fileHandle.close();
+ if (!hclusterWanted) { delete matrix; }
+ else { m->mothurRemove(overlapFile); m->mothurRemove(distFile); }
+ delete reading;
+ return 0;
+ }
+
+ reading->finish();
+ delete reading;
+ fileHandle.close();
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ReadBlast", "read");
+ exit(1);
+ }
+}
+/*********************************************************************************************/
+int ReadBlast::readNames(NameAssignment* nameMap) {
+ try {
+ m->mothurOut("Reading names... "); cout.flush();
+
+ string name, hold, prevName;
+ int num = 1;
+
+ ifstream in;
+ m->openInputFile(blastfile, in);
+
+ //ofstream outName;
+ //m->openOutputFile((blastfile + ".tempOutNames"), outName);
+
+ //read first line
+ in >> prevName;
+
+ for (int i = 0; i < 11; i++) { in >> hold; }
+ m->gobble(in);
+
+ //save name in nameMap
+ nameMap->push_back(prevName);
+
+ while (!in.eof()) {
+ if (m->control_pressed) { in.close(); return 0; }
+
+ //read line
+ in >> name;
+
+ for (int i = 0; i < 11; i++) { in >> hold; }
+ m->gobble(in);
+
+ //is this a new name?
+ if (name != prevName) {
+ prevName = name;
+ nameMap->push_back(name);
+ num++;
+ }
+ }
+
+ in.close();
+
+ //write out names file
+ //string outNames = m->getRootName(blastfile) + "names";
+ //ofstream out;
+ //m->openOutputFile(outNames, out);
+ //nameMap->print(out);
+ //out.close();
+
+ if (m->control_pressed) { return 0; }
+
+ m->mothurOut(toString(num) + " names read."); m->mothurOutEndLine();
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ReadBlast", "readNames");
+ exit(1);
+ }
+}
+/*********************************************************************************************/
+
+
+
projects / mothur.git / blobdiff