Smith Waterman für C++



  • Nachdem ich dem Code aus github gefolgt bin (https://github.com/ding-lab/pindel2/blob/master/src/smith_waterman_alignment.cpp) habe ich Folgendes implementiert:

    int BinaryAlignment::watermanAlign(std::string seqSynth, std::string seqReal, bool viableAlignment) {
    	const int nSynth = seqSynth.length();                                               //length of sequences
    	const int nReal = seqReal.length();
    
    	//H[nSynth + 1][nReal + 1]
    	std::vector<std::vector<double>> H(nReal + 1, std::vector<double>(nSynth + 1));
    	/*double **H = new double*[nReal + 1];
    	for (int i = 0; i < nReal + 1; i++)
    		H[i] = new double[nSynth + 1];*/
    
    	for (int m = 0; m <= nSynth; m++)
    		for (int n = 0; n <= nReal; n++)
    			H[m][n] = 0;
    
    	double temp[4];
    	std::vector<std::vector<double>> Ii(nReal + 1, std::vector<double>(nSynth + 1));
    	/*int **Ii = new int*[nReal + 1];
    	for (int i = 0; i < nReal + 1; i++)
    		Ii[i] = new int[nSynth + 1];*/
    
    	std::vector<std::vector<double>> Ij(nReal + 1, std::vector<double>(nSynth + 1));
    	/*int **Ij = new int*[nReal + 1];
    	for (int i = 0; i < nReal + 1; i++)
    		Ij[i] = new int[nSynth + 1];
    */
    	for (int i = 1; i <= nSynth; i++) {
    		for (int j = 1; j <= nReal; j++) {
    			temp[0] = H[i - 1][j - 1] + similarityScore(seqSynth[i - 1], seqReal[j - 1]);
    			temp[1] = H[i - 1][j] - delta;
    			temp[2] = H[i][j - 1] - delta;
    			temp[3] = 0;
    			H[i][j] = findArrayMax(temp, 4);
    
    			switch (ind) {
    			case 0:                                  // score in (i,j) stems from a match/mismatch
    				Ii[i][j] = i - 1;
    				Ij[i][j] = j - 1;
    				break;
    			case 1:                                  // score in (i,j) stems from a deletion in sequence A
    				Ii[i][j] = i - 1;
    				Ij[i][j] = j;
    				break;
    			case 2:                                  // score in (i,j) stems from a deletion in sequence B
    				Ii[i][j] = i;
    				Ij[i][j] = j - 1;
    				break;
    			case 3:                                  // (i,j) is the beginning of a subsequence
    				Ii[i][j] = i;
    				Ij[i][j] = j;
    				break;
    			}
    		}
    	}
    
    	//Print matrix H to console 
    	std::cout << "**********************************************" << std::endl;
    	std::cout << "The scoring matrix is given by  " << std::endl << std::endl;
    	for (int i = 1; i <= nSynth; i++) {
    		for (int j = 1; j <= nReal; j++) {
    			std::cout << H[i][j] << " ";
    		}
    		std::cout << std::endl;
    	}
    
    	//search H for the moaximal score
    	double Hmax = 0;
    	int imax = 0, jmax = 0;
    
    	for (int i = 1; i <= nSynth; i++) {
    		for (int j = 1; j <= nReal; j++) {
    			if (H[i][j] > Hmax) {
    				Hmax = H[i][j];
    				imax = i;
    				jmax = j;
    			}
    		}
    	}
    
    	std::cout << Hmax << endl;
    	std::cout << nSynth << ", " << nReal << ", " << imax << ", " << jmax << std::endl;
    	std::cout << "max score: " << Hmax << std::endl;
    	std::cout << "alignment index: " << (imax - jmax) << std::endl;
    
    	//Backtracing from Hmax
    	int icurrent = imax, jcurrent = jmax;
    	int inext = Ii[icurrent][jcurrent];
    	int jnext = Ij[icurrent][jcurrent];
    	int tick = 0;
    	char *consensusSynth = new char[nSynth + nReal + 2];
    	char *consensusReal = new char[nSynth + nReal + 2];
    
    	while (((icurrent != inext) || (jcurrent != jnext)) && (jnext >= 0) && (inext >= 0)) {
    
    		if (inext == icurrent)
    			consensusSynth[tick] = '-';                         //deletion in A
    		else
    			consensusSynth[tick] = seqSynth[icurrent - 1];      //match / mismatch in A
    
    		if (jnext == jcurrent)
    			consensusReal[tick] = '-';                          //deletion in B
    		else
    			consensusReal[tick] = seqReal[jcurrent - 1];        //match/mismatch in B
    
    																//fix for adding first character of the alignment.
    		if (inext == 0)
    			inext = -1;
    		else if (jnext == 0)
    			jnext = -1;
    		else
    			icurrent = inext;
    		jcurrent = jnext;
    
    		inext = Ii[icurrent][jcurrent];
    		jnext = Ij[icurrent][jcurrent];
    
    		tick++;
    	}
    
    	// Output of the consensus motif to the console
    	std::cout << std::endl << "***********************************************" << std::endl;
    	std::cout << "The alignment of the sequences" << std::endl << std::endl;
    	for (int i = 0; i < nSynth; i++) {
    		std::cout << seqSynth[i];
    	};
    	std::cout << "  and" << std::endl;
    	for (int i = 0; i < nReal; i++) {
    		std::cout << seqReal[i];
    	};
    	std::cout << std::endl << std::endl;
    	std::cout << "is for the parameters  mu = " << mu << " and delta = " << delta << " given by" << std::endl << std::endl;
    	for (int i = tick - 1; i >= 0; i--)
    		std::cout << consensusSynth[i];
    	std::cout << std::endl;
    	for (int j = tick - 1; j >= 0; j--)
    		std::cout << consensusReal[j];
    	std::cout << std::endl;
    
    	int numMismatches = 0;
    	for (int i = tick - 1; i >= 0; i--) {
    		if (consensusSynth[i] != consensusReal[i]) {
    			numMismatches++;
    		}
    	}
    
    	viableAlignment = numMismatches <= maxMismatch;
    	return imax - jmax;
    }
    

    Nun probiere ich daran seit einer Stunde rum und ich kriege immer wieder Exceptions an verschiedenen Stellen, diese sind einmal

    inext = Ii[icurrent][jcurrent];
    

    und

    consensusReal[tick] = '-';
    

    . Die Exception sagt mir folgendes: Ausnahme ausgelöst bei 0x00007FF637292246 in Waterman.exe: 0xC0000005: Zugriffsverletzung beim Schreiben an Position 0x000002006A17E000.

    Als Parameter habe ich auch bereits jegliche verschiedene Strings benutzt (für mein Vorhaben sind es immer binäre String-Sequenzen). Einmal hat's geklappt, danach wurde rein gar nichts an dem Code geändert und nun läuft's gar nciht mehr. Was mache ich falsch? 😕



  • new char - Warum nicht std::vector?
    Benutze einen Debugger.
    Ersetzte Ii[icurrent][jcurrent] durch Ii.at(icurrent).at(jcurrent) (auch an anderen Stellen) und warte auf die Exception.



  • Debuggen hab ich bereits getan, kommt auf's gleiche raus. Und mit dem Ersetzen wird es leider auch nicht besser, da kommt die gleiche Exception.



  • vickal93 schrieb:

    Debuggen hab ich bereits getan, kommt auf's gleiche raus. Und mit dem Ersetzen wird es leider auch nicht besser, da kommt die gleiche Exception.

    Die Exception kommt daher das du den Bereich verlässt. Sprich du greifst auf einen Index des Vectors zu der nicht vorhanden ist weil icurrent warscheinlich größer wird als die Größe deines Vectors.



  • axels. schrieb:

    vickal93 schrieb:

    Debuggen hab ich bereits getan, kommt auf's gleiche raus. Und mit dem Ersetzen wird es leider auch nicht besser, da kommt die gleiche Exception.

    Die Exception kommt daher das du den Bereich verlässt. Sprich du greifst auf einen Index des Vectors zu der nicht vorhanden ist weil icurrent warscheinlich größer wird als die Größe deines Vectors.

    Aber ich initialisiere doch icurrent = imax? Bzw. was müsste ich dann ändern?



  • Wenn du alle Stellen ersetzt, bei denen auf einen vector zugegriffen wird und auch das new char entfernst, kommt eine andere Exception.



  • manni66 schrieb:

    Wenn du alle Stellen ersetzt, bei denen auf einen vector zugegriffen wird und auch das new char entfernst, kommt eine andere Exception.

    Nein, da kommt genau die gleiche Exception.



  • Dann hast du irgendwo außerhalb der Funktion schon etwas kaputt gemacht.



  • vickal93 schrieb:

    Debuggen hab ich bereits getan, kommt auf's gleiche raus.

    Was soll das heissen "kommt auf's gleiche raus". Klar passiert das selbe im Debugger. Nur siehst du dann wo es knallt. Also genau die Zeile. Und du kannst gucken was deine ganzen Variablen alle für Werte haben zu dem Zeitpunkt.


Anmelden zum Antworten