5 * Created by Pat Schloss on 12/16/08.
6 * Copyright 2008 Patrick D. Schloss. All rights reserved.
8 * This is a basic alignment method that gets the blast program to do the heavy lifting. In the future, we should
9 * probably incorporate NCBI's library so that we don't have to call on a user-supplied executable. This is a child
10 * of the Alignment class, which requires a constructor and align method.
15 #include "alignment.hpp"
16 #include "blastalign.hpp"
19 //**************************************************************************************************/
21 BlastAlignment::BlastAlignment(float go, float ge, float m, float mm) :
22 match(m), // This is the score to award for two nucleotides matching (match >= 0)
23 mismatch(mm) // This is the penalty to assess for a mismatch (mismatch <= 0)
25 globaldata = GlobalData::getInstance();
26 path = globaldata->argv;
27 path = path.substr(0, (path.find_last_of('m')));
29 gapOpen = abs(go); // This is the penalty to assess for opening a gap (gapOpen >= 0)
30 gapExtend = abs(ge); // This is the penalty to assess for extending a gap (gapExtend >= 0)
32 int randNumber = rand();
33 candidateFileName = toString(randNumber) + ".candidate";
34 templateFileName = toString(randNumber) + ".template";
35 blastFileName = toString(randNumber) + ".pairwise";
38 //**************************************************************************************************/
40 BlastAlignment::~BlastAlignment(){ // The desctructor should clean up by removing the temporary
41 remove(candidateFileName.c_str()); // files used to run bl2seq
42 remove(templateFileName.c_str());
43 remove(blastFileName.c_str());
46 //**************************************************************************************************/
48 void BlastAlignment::align(string seqA, string seqB){ //Use blastn to align the two sequences
50 ofstream candidateFile(candidateFileName.c_str()); // Write the sequence to be aligned to a temporary candidate seq file
51 candidateFile << ">candidate" << endl << seqA << endl;
52 candidateFile.close();
54 ofstream templateFile(templateFileName.c_str()); // Write the unaligned template sequence to a temporary candidate seq file
55 templateFile << ">template" << endl << seqB << endl;
58 // The blastCommand assumes that we have DNA sequences (blastn) and that they are fairly similar (-e 0.001) and
59 // that we don't want to apply any kind of complexity filtering (-F F)
60 string blastCommand = path + "blast/bin/bl2seq -p blastn -i " + candidateFileName + " -j " + templateFileName + " -e 0.0001 -F F -o " + blastFileName + " -W 11";
61 blastCommand += " -r " + toString(match) + " -q " + toString(mismatch);
62 blastCommand += " -G " + toString(gapOpen) + " -E " + toString(gapExtend);
64 system(blastCommand.c_str()); // Here we assume that "bl2seq" is in the users path or in the same folder as
69 /**************************************************************************************************/
71 void BlastAlignment::setPairwiseSeqs(){ // This method call assigns the blast generated alignment
72 // to the pairwise entry in the Sequence class for the
73 // candidate and template Sequence objects
75 m->openInputFile(blastFileName, blastFile);
80 int candidateLength, templateLength;
83 string candidateName, templateName;
85 while((d=blastFile.get()) != '='){}
86 blastFile >> candidateName; // Get the candidate sequence name from flatfile
88 while((d=blastFile.get()) != '('){}
89 blastFile >> candidateLength; // Get the candidate sequence length from flatfile
91 while((d=blastFile.get())){
93 blastFile >> templateName; // Get the template sequence name from flatfile
96 else if(d == '*'){ // We go here if there is no significant match
105 // while(dummy != "query:"){ m->mothurOut(dummy, ""); m->mothurOutEndLine(); blastFile >> dummy; }
106 // blastFile >> seqBend;
107 // m->mothurOut(toString(seqBend), ""); m->mothurOutEndLine();
108 // for(int i=0;i<seqBend;i++){
112 // pairwiseLength = 0;
117 while((d=blastFile.get()) != '='){}
118 blastFile >> templateLength; // Get the template sequence length from flatfile
120 while((d=blastFile.get()) != 'Q'){} // Suck up everything else until we get to the start of the alignment
121 int queryStart, sbjctStart, queryEnd, sbjctEnd;
122 string queryLabel, sbjctLabel, query, sbjct;
124 blastFile >> queryLabel; queryLabel = 'Q' + queryLabel;
127 while(queryLabel == "Query:"){
128 blastFile >> queryStart >> query >> queryEnd;
130 while((d=blastFile.get()) != 'S'){};
132 blastFile >> sbjctLabel >> sbjctStart >> sbjct >> sbjctEnd;
135 seqAstart = queryStart;
136 seqBstart = sbjctStart;
139 seqAaln += query; // concatenate each line of the sequence to what we already have
140 seqBaln += sbjct; // for the query and template (subject) sequence
142 blastFile >> queryLabel;
146 pairwiseLength = seqAaln.length();
148 for(int i=1;i<seqBstart;i++){ // Since the alignments don't always start at (1, 1), we need to pad
149 seqAaln = 'Z' + seqAaln; // the sequences so that they start at the same point
150 seqBaln = 'X' + seqBaln;
153 for(int i=seqBend+1;i<=templateLength;i++){ // since the sequences don't necessarily end at the same point, we
154 seqAaln += 'Z'; // again need ot pad the sequences so that they extend to the length
155 seqBaln += 'X'; // of the template sequence
160 //**************************************************************************************************/