Reputation: 45
EDIT: I found the error! The program wasn't actually getting stuck in the loop in the first place. Really it was getting stuck in a loop further down in the code, but because there was no newline at the end of my confirmation cout at the end of this loop, it never executed. Thank you for all the help!
I'm working on a program that will align two nucleotide sequences. It's working well up until this point. The loop executed up until 661, the size of seq1length being 663, and then doesn't leave the loop. I put a fail safe into the loop so that it won't go on forever, but it doesn't appear to work. I feel like I'm missing something silly and I would really appreciate another pair of eyes.
//calculate the matrix using the model. Figure out whether a up, left or diagonal move would be best
int CalcMatrix [seq1length][seq2length]; //this will store the cumulative alignment scores at each spot
int DirectionalMatrix [seq1length][seq2length] ; //this will store the direction of the alignment //0 is diagonal, 1 is up, 2 is left
int diag;
int up;
int left;
int incase = 0; //set up a fail safe in case it gets too big
while(incase < 2*(seq1length+seq2length)) // to make sure the loop doesn't go on forever. doesn't work
{
for(int i = 0; i < seq1length-1; i++) //for the length of the first sequence
{
CalcMatrix[i][0] = i*gappenalty; //first row is a cumulative gap penalty
DirectionalMatrix[i][0] = 2; //to the left
for(int j = 0; j < seq2length-1; j++)
{
CalcMatrix[0][j] = j*gappenalty; //first column is a cumulative gap penalty
DirectionalMatrix[0][j] = 1; // to the up
if(i > 0 && j > 0) //if i and j are greater than 0, that is we aren't in the first row or column, make the rest of the matrix
{
diag = CalcMatrix[i-1][j-1] + ScoreMatrix[i][j];
up = CalcMatrix[i-1][j] + gappenalty;
left = CalcMatrix[i][j-1] + gappenalty;
CalcMatrix[i][j] = max(diag, max(up,left));
if(diag == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 0;}
else if(up == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 1;}
else if(left == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 2;}
incase++;
}
}
cout << "Calculating Alignment... Please Hold... " << i << " Nucleotides Aligned \n";
}
}
cout << "Yay! You got out of the loop!";
EDIT
I made a self contained compilable version. However this one does leave the loop successfully. Which means it's either an issue with the size of the values I'm working with in the original code, or because of how the values are initialized higher up in the original code
#include <iostream>
#include <string>
using namespace std;
int main(){
//build matrix of scores
string seq1seq = "ACTGACTGACTGACTGACTACTG";
string seq2seq = "ACTGTGTCTTCTGATCTCGATCCA";
int seq1length = seq1seq.length();
int seq2length = seq2seq.length();
int gappenalty = -40;
cout << "Sequence 1 length is: " << seq1length << "\n";
cout << "Sequence 2 length is: " << seq2length << "\n";
int incase = 0;
int a [21] = {0};
int ScoreMatrix[seq1length][seq2length];
for (int i = 0; i < seq1length; i++)
{
for(int j = 0; j < seq2length; j++)
{
if(seq1seq[i] == 'A')
{
if(seq2seq[j] == 'A')
{
ScoreMatrix[0][0] = a[5];//the substitution score array at [0][0] ;
}
else if(seq2seq[j] == 'C')
{
ScoreMatrix[0][1] = a[6];//the substitution score array at [0][1] that is, An A in sequence 1 and a C in sequence 2
}
else if(seq2seq[j] == 'G')
{
ScoreMatrix[0][2] = a[7];//the substitution score array at [0][2]
}
else if(seq2seq[j] == 'T')
{
ScoreMatrix[0][3] = a[8];//the substitution score array at [0][3]
}
}
if(seq1seq[i] == 'C')
{
if(seq2seq[j] == 'A')
{
ScoreMatrix[1][0] = a[9];//the substitution score array at [1][0]
}
else if(seq2seq[j] == 'C')
{
ScoreMatrix[1][1] = a[10];//the substitution score array at [1][1] that is, An C in sequence 1 and a C in sequence 2
}
else if(seq2seq[j] == 'G')
{
ScoreMatrix[1][2] = a[11];//the substitution score array at [1][2]
}
else if(seq2seq[j] == 'T')
{
ScoreMatrix[1][3] = a[12];//the substitution score array at [1][3
}
}
if(seq1seq[i] == 'G')
{
if(seq2seq[j] == 'A')
{
ScoreMatrix[2][0] = a[13];//the substitution score array at [2][0]
}
else if(seq2seq[j] == 'C')
{
ScoreMatrix[2][1] = a[14];//the substitution score array at [2][1] that is, An G in sequence 1 and a C in sequence 2
}
else if(seq2seq[j] == 'G')
{
ScoreMatrix[2][2] = a[15];//the substitution score array at [2][2]
}
else if(seq2seq[j] == 'T')
{
ScoreMatrix[2][3] = a[16];//the substitution score array at [2][3]
}
}
if(seq1seq[i] == 'A')
{
if(seq2seq[j] == 'A')
{
ScoreMatrix[3][0] = a[17];//the substitution score array at [3][0]
}
else if(seq2seq[j] == 'C')
{
ScoreMatrix[3][1] = a[18];//the substitution score array at [3][1] that is, An T in sequence 1 and a C in sequence 2
}
else if(seq2seq[j] == 'G')
{
ScoreMatrix[3][2] = a[19];//the substitution score array at [3][2]
}
else if(seq2seq[j] == 'T')
{
ScoreMatrix[3][3] = a[20];//the substitution score array at [3][3]
}
}
}
cout << "Assembling Score Matrix.... Please Hold.... " << i << " Nucleotides Covered \n" ;
}
cout << "YAY! Got out of the scoring assembly. Aren't you proud of me?" ;
//calculate the matrix using the silly named model. Figure out whether a up, left or diagonal move would be best
int CalcMatrix [seq1length][seq2length]; //this will store the cumulative alignment scores at each spot
int DirectionalMatrix [seq1length][seq2length] ; //this will store the direction of the alignment //0 is diagonal, 1 is up, 2 is left
int diag;
int up;
int left;
incase = 0; //set up a fail safe incase it gets too big
for(int i = 0; i < seq1length; i++) //for the length of the first sequence
{
CalcMatrix[i][0] = i*gappenalty; //first row is a cumulative gap penalty
DirectionalMatrix[i][0] = 2; //to the left
for(int j = 0; j < seq2length; j++)
{
CalcMatrix[0][j] = j*gappenalty; //first column is a cumulative gap penalty
DirectionalMatrix[0][j] = 1; // to the up
if(i > 0 && j > 0) //if i and j are greater than 0, that is we aren't in the first row or column, make the rest of the matrix
{
diag = CalcMatrix[i-1][j-1] + ScoreMatrix[i][j];
up = CalcMatrix[i-1][j] + gappenalty;
left = CalcMatrix[i][j-1] + gappenalty;
CalcMatrix[i][j] = max(diag, max(up,left));
if(diag == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 0;}
else if(up == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 1;}
else if(left == CalcMatrix[i][j]) { DirectionalMatrix[i][j] = 2;}
}
}
//cout << "Calculating Alignment... Please Hold... " << i << " Nucleotides Aligned \n";
}
cout << "Yay! You got out of the loop!";
return 0;}
Upvotes: 0
Views: 157
Reputation: 1434
Because you put your incase++ inside the if statment of your nested for loop here:
for(int i = 0; i < seq1length-1; i++){
//Your Code
for(int j = 0; j < seq2length-1; j++){
//Your Code
incase++;
}
}
Your innermost nested for loop only loops until less than the seq2length-1, and then your outer for loop only loops until less than the seq1length-1, it can't reach your 2*(seq1length+seq2length) limit in your while loop, even if your if statement always got executed.
Let's say your innermost for loop loops 1 times, and your outer for loop loops 663 times, in total your incase value will be incremented to 663*1 = 663. But your while loop condition is 2*(664+2), which is 1332. This will never get your incase to hit the limit of the while loop.
In addition, you put your incase++ inside your if statement which doesn't execute until both of your for loops i and j are greater than 0, which might make your incase increments to only less than 663 times in total (which explained why your incase has been incremented to only 661) and it stops right there, without reaching your while loop limit. For example, when either your i or j is 0, your for loops even though are looping, but your incase won't get incremented because i and j are not both greater than 0 (the if statement condition).
You can either put your incase++ outside of both the for loops but in the while loop and set the condition correctly (depending on what your program is supposed to do), or your while loop condition needs to be set correctly, taking into account that your if statement doesn't make incase increment until both your i and j are greater than 0. Or put it just outside of your if statement, or change your if statement condition. Depend on what you want to do.
Upvotes: 1
Reputation: 5291
I ran this reduced test application from your code, and it seems to exit fine:
#include <iostream>
using namespace std;
#define seq1length 3
#define seq2length 3
int main(int argc, char *argv[])
{
int incase = 0; //set up a fail safe incase it gets too big
while (incase < 2 * (seq1length + seq2length)) // to make sure the loop doesn't go on forever. doesn't work
{
for (int i = 0; i < seq1length - 1; i++) //for the length of the first sequence
{
for (int j = 0; j < seq2length - 1; j++)
{
if (i > 0 && j > 0) //if i and j are greater than 0, that is we aren't in the first row or column, make the rest of the matrix
{
incase++;
}
}
cout << "Calculating Alignment... Please Hold... " << i << " Nucleotides Aligned \n";
}
}
cout << "Yay! You got out of the loop!";
getchar();
return(0);
}
I have removed the code that is not dependant on the loop conditin check. Can you tell a bit more about your problem...Like what values are you passing for seq1length and seq2length...
Here is the output i saw:
After refering to http://ideone.com/Zrvzpa from your comment, i see that it is exiting at 661 because your condition is i < seq1length-1;
and j < seq2length-1;
.
If you change it to i < seq1length;
and j < seq2length;
. It will loop till 662.
Note that it should not loop till 663 as the array is indexed from 0 - 662.
Upvotes: 1