5 * Created by westcott on 9/9/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "sequenceparser.h"
13 /************************************************************/
14 SequenceParser::SequenceParser(string groupFile, string fastaFile, string nameFile) {
17 m = MothurOut::getInstance();
21 groupMap = new GroupMap(groupFile);
22 error = groupMap->readMap();
24 if (error == 1) { m->control_pressed = true; }
27 vector<string> namesOfGroups = groupMap->getNamesOfGroups();
28 for (int i = 0; i < namesOfGroups.size(); i++) {
29 vector<Sequence> temp;
30 map<string, string> tempMap;
31 seqs[namesOfGroups[i]] = temp;
32 nameMapPerGroup[namesOfGroups[i]] = tempMap;
35 //read fasta file making sure each sequence is in the group file
37 m->openInputFile(fastaFile, in);
39 map<string, string> seqName; //stores name -> sequence string so we can make new "unique" sequences when we parse the name file
42 if (m->control_pressed) { break; }
44 Sequence seq(in); m->gobble(in);
46 if (seq.getName() != "") {
48 string group = groupMap->getGroup(seq.getName());
49 if (group == "not found") { error = 1; m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your groupfile, please correct."); m->mothurOutEndLine(); }
51 seqs[group].push_back(seq);
52 seqName[seq.getName()] = seq.getAligned();
58 if (error == 1) { m->control_pressed = true; }
62 m->openInputFile(nameFile, inName);
66 set<string> thisnames1;
67 while(!inName.eof()) {
69 if (m->control_pressed) { break; }
71 inName >> first; m->gobble(inName);
72 inName >> second; m->gobble(inName);
75 m->splitAtChar(second, names, ',');
77 //get aligned string for these seqs from the fasta file
78 string alignedString = "";
79 map<string, string>::iterator itAligned = seqName.find(names[0]);
80 if (itAligned == seqName.end()) {
81 error = 1; m->mothurOut("[ERROR]: " + names[0] + " is in your name file and not in your fasta file, please correct."); m->mothurOutEndLine();
83 alignedString = itAligned->second;
86 //separate by group - parse one line in name file
87 map<string, string> splitMap; //group -> name1,name2,...
88 map<string, string>::iterator it;
89 for (int i = 0; i < names.size(); i++) {
91 string group = groupMap->getGroup(names[i]);
92 if (group == "not found") { error = 1; m->mothurOut("[ERROR]: " + names[i] + " is in your name file and not in your groupfile, please correct."); m->mothurOutEndLine(); }
95 it = splitMap.find(group);
96 if (it != splitMap.end()) { //adding seqs to this group
97 (it->second) += "," + names[i];
98 thisnames1.insert(names[i]);
100 }else { //first sighting of this group
101 splitMap[group] = names[i];
103 thisnames1.insert(names[i]);
105 //is this seq in the fasta file?
106 if (i != 0) { //if not then we need to add a duplicate sequence to the seqs for this group so the new "fasta" and "name" files will match
107 Sequence tempSeq(names[i], alignedString); //get the first guys sequence string since he's in the fasta file.
108 seqs[group].push_back(tempSeq);
113 allSeqsMap[names[i]] = names[0];
117 //fill nameMapPerGroup - holds all lines in namefile separated by group
118 for (it = splitMap.begin(); it != splitMap.end(); it++) {
120 string firstName = "";
121 for(int i = 0; i < (it->second).length(); i++) {
122 if (((it->second)[i]) != ',') {
123 firstName += ((it->second)[i]);
127 //group1 -> seq1 -> seq1,seq2,seq3
128 nameMapPerGroup[it->first][firstName] = it->second;
134 if (error == 1) { m->control_pressed = true; }
136 if (countName != (groupMap->getNumSeqs())) {
137 vector<string> groupseqsnames = groupMap->getNamesSeqs();
138 for (int i = 0; i < groupseqsnames.size(); i++) {
139 set<string>::iterator itnamesfile = thisnames1.find(groupseqsnames[i]);
140 if (itnamesfile == thisnames1.end()){
141 cout << "missing name " + groupseqsnames[i] << '\t' << allSeqsMap[groupseqsnames[i]] << endl;
144 m->mothurOutEndLine();
145 m->mothurOut("[ERROR]: Your name file contains " + toString(countName) + " valid sequences, and your groupfile contains " + toString(groupMap->getNumSeqs()) + ", please correct.");
146 m->mothurOutEndLine();
147 m->control_pressed = true;
151 catch(exception& e) {
152 m->errorOut(e, "SequenceParser", "SequenceParser");
156 /************************************************************/
157 SequenceParser::SequenceParser(string groupFile, string fastaFile) {
160 m = MothurOut::getInstance();
164 groupMap = new GroupMap(groupFile);
165 error = groupMap->readMap();
167 if (error == 1) { m->control_pressed = true; }
170 vector<string> namesOfGroups = groupMap->getNamesOfGroups();
171 for (int i = 0; i < namesOfGroups.size(); i++) {
172 vector<Sequence> temp;
173 seqs[namesOfGroups[i]] = temp;
176 //read fasta file making sure each sequence is in the group file
178 m->openInputFile(fastaFile, in);
183 if (m->control_pressed) { break; }
185 Sequence seq(in); m->gobble(in);
187 if (seq.getName() != "") {
189 string group = groupMap->getGroup(seq.getName());
190 if (group == "not found") { error = 1; m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your groupfile, please correct."); m->mothurOutEndLine(); }
191 else { seqs[group].push_back(seq); count++; }
196 if (error == 1) { m->control_pressed = true; }
198 if (count != (groupMap->getNumSeqs())) {
199 m->mothurOutEndLine();
200 m->mothurOut("[ERROR]: Your fasta file contains " + toString(count) + " valid sequences, and your groupfile contains " + toString(groupMap->getNumSeqs()) + ", please correct.");
201 if (count < (groupMap->getNumSeqs())) { m->mothurOut(" Did you forget to include the name file?"); }
202 m->mothurOutEndLine();
203 m->control_pressed = true;
207 catch(exception& e) {
208 m->errorOut(e, "SequenceParser", "SequenceParser");
212 /************************************************************/
213 SequenceParser::~SequenceParser(){ delete groupMap; }
214 /************************************************************/
215 int SequenceParser::getNumGroups(){ return groupMap->getNumGroups(); }
216 /************************************************************/
217 vector<string> SequenceParser::getNamesOfGroups(){ return groupMap->getNamesOfGroups(); }
218 /************************************************************/
219 bool SequenceParser::isValidGroup(string g){ return groupMap->isValidGroup(g); }
220 /************************************************************/
221 string SequenceParser::getGroup(string g){ return groupMap->getGroup(g); }
222 /************************************************************/
223 int SequenceParser::getNumSeqs(string g){
225 map<string, vector<Sequence> >::iterator it;
229 if(it == seqs.end()) {
230 m->mothurOut("[ERROR]: " + g + " is not a valid group, please correct."); m->mothurOutEndLine();
232 num = (it->second).size();
237 catch(exception& e) {
238 m->errorOut(e, "SequenceParser", "getNumSeqs");
242 /************************************************************/
243 vector<Sequence> SequenceParser::getSeqs(string g){
245 map<string, vector<Sequence> >::iterator it;
246 vector<Sequence> seqForThisGroup;
249 if(it == seqs.end()) {
250 m->mothurOut("[ERROR]: No sequences available for group " + g + ", please correct."); m->mothurOutEndLine();
252 seqForThisGroup = it->second;
255 return seqForThisGroup;
257 catch(exception& e) {
258 m->errorOut(e, "SequenceParser", "getSeqs");
262 /************************************************************/
263 int SequenceParser::getSeqs(string g, string filename, bool uchimeFormat=false){
265 map<string, vector<Sequence> >::iterator it;
266 vector<Sequence> seqForThisGroup;
267 vector<seqPriorityNode> nameVector;
270 if(it == seqs.end()) {
271 m->mothurOut("[ERROR]: No sequences available for group " + g + ", please correct."); m->mothurOutEndLine();
275 m->openOutputFile(filename, out);
277 seqForThisGroup = it->second;
280 // format should look like
281 //>seqName /ab=numRedundantSeqs/
284 map<string, string> nameMapForThisGroup = getNameMap(g);
285 map<string, string>::iterator itNameMap;
288 for (int i = 0; i < seqForThisGroup.size(); i++) {
289 itNameMap = nameMapForThisGroup.find(seqForThisGroup[i].getName());
291 if (itNameMap == nameMapForThisGroup.end()){
293 m->mothurOut("[ERROR]: " + seqForThisGroup[i].getName() + " is in your fastafile, but is not in your namesfile, please correct."); m->mothurOutEndLine();
295 int num = m->getNumNames(itNameMap->second);
297 seqPriorityNode temp(num, seqForThisGroup[i].getAligned(), seqForThisGroup[i].getName());
298 nameVector.push_back(temp);
302 if (error == 1) { out.close(); m->mothurRemove(filename); return 1; }
304 //sort by num represented
305 sort(nameVector.begin(), nameVector.end(), compareSeqPriorityNodes);
307 //print new file in order of
308 for (int i = 0; i < nameVector.size(); i++) {
310 if(m->control_pressed) { out.close(); m->mothurRemove(filename); return 1; }
312 out << ">" << nameVector[i].name << "/ab=" << nameVector[i].numIdentical << "/" << endl << nameVector[i].seq << endl;
316 for (int i = 0; i < seqForThisGroup.size(); i++) {
318 if(m->control_pressed) { out.close(); m->mothurRemove(filename); return 1; }
320 seqForThisGroup[i].printSequence(out);
328 catch(exception& e) {
329 m->errorOut(e, "SequenceParser", "getSeqs");
334 /************************************************************/
335 map<string, string> SequenceParser::getNameMap(string g){
337 map<string, map<string, string> >::iterator it;
338 map<string, string> nameMapForThisGroup;
340 it = nameMapPerGroup.find(g);
341 if(it == nameMapPerGroup.end()) {
342 m->mothurOut("[ERROR]: No nameMap available for group " + g + ", please correct."); m->mothurOutEndLine();
344 nameMapForThisGroup = it->second;
347 return nameMapForThisGroup;
349 catch(exception& e) {
350 m->errorOut(e, "SequenceParser", "getNameMap");
354 /************************************************************/
355 int SequenceParser::getNameMap(string g, string filename){
357 map<string, map<string, string> >::iterator it;
358 map<string, string> nameMapForThisGroup;
360 it = nameMapPerGroup.find(g);
361 if(it == nameMapPerGroup.end()) {
362 m->mothurOut("[ERROR]: No nameMap available for group " + g + ", please correct."); m->mothurOutEndLine();
364 nameMapForThisGroup = it->second;
367 m->openOutputFile(filename, out);
369 for (map<string, string>::iterator itFile = nameMapForThisGroup.begin(); itFile != nameMapForThisGroup.end(); itFile++) {
371 if(m->control_pressed) { out.close(); m->mothurRemove(filename); return 1; }
373 out << itFile->first << '\t' << itFile->second << endl;
381 catch(exception& e) {
382 m->errorOut(e, "SequenceParser", "getNameMap");
386 /************************************************************/