Context Navigation

source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWatermanRepeated.java @ 1587

Last change on this file since 1587 was 1587, checked in by rkrimmel, 10 years ago
Adding needleman wunsch
File size: 8.6 KB

Rev	Line
[1558]	1	package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
	2
	3	import java.util.ArrayList;
[1575]	4	import java.util.logging.Level;
[1558]	5
[1572]	6	import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
[1578]	7	import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
[1575]	8	import de.ugoe.cs.util.console.Console;
[1558]	9
[1586]	10	public class SmithWatermanRepeated implements AlignmentAlgorithm {
[1558]	11
	12	/**
	13	* The first input
	14	*/
	15	private int[] input1;
	16
	17	/**
	18	* The second input String
	19	*/
	20	private int[] input2;
	21
	22	/**
	23	* The lengths of the input
	24	*/
	25	private int length1, length2;
	26
	27	/**
	28	* The score matrix. The true scores should be divided by the normalization
	29	* factor.
	30	*/
	31	private MatrixEntry[][] matrix;
	32
[1572]	33
[1585]	34	private ArrayList<NumberSequence> alignment;
[1558]	35
	36	private float scoreThreshold;
	37
	38	/**
	39	* Substitution matrix to calculate scores
	40	*/
	41	private SubstitutionMatrix submat;
	42
	43	public SmithWatermanRepeated(int[] input1, int[] input2, SubstitutionMatrix submat,float threshold) {
	44	this.input1 = input1;
	45	this.input2 = input2;
	46	length1 = input1.length;
	47	length2 = input2.length;
	48	this.submat = submat;
	49
	50	//System.out.println("Starting SmithWaterman algorithm with a "
	51	// + submat.getClass() + " Substitution Matrix: " + submat.getClass().getCanonicalName());
	52	this.scoreThreshold = threshold;
	53
	54	matrix = new MatrixEntry[length1+2][length2+1];
[1572]	55	alignment = new ArrayList<NumberSequence>();
[1558]	56
[1559]	57	for (int i = 0; i <= length1+1; i++) {
[1558]	58	for(int j = 0; j< length2; j++) {
	59	matrix[i][j] = new MatrixEntry();
	60	}
	61	}
	62
	63
	64	buildMatrix();
[1574]	65	traceback();
[1558]	66	}
	67
	68	/**
	69	* Compute the similarity score of substitution The position of the first
	70	* character is 1. A position of 0 represents a gap.
	71	*
	72	* @param i
	73	* Position of the character in str1
	74	* @param j
	75	* Position of the character in str2
	76	* @return Cost of substitution of the character in str1 by the one in str2
	77	*/
[1568]	78	private double similarity(int i, int j) {
[1578]	79	return submat.getScore(input1[i - 1], input2[j - 1]);
[1558]	80	}
	81
	82	/**
[1559]	83	* Build the score matrix using dynamic programming.
[1558]	84	*/
	85	private void buildMatrix() {
	86	if (submat.getGapPenalty() >= 0) {
	87	throw new Error("Indel score must be negative");
	88	}
	89
[1559]	90	// it's a gap
[1558]	91	matrix[0][0].setScore(0);
	92	matrix[0][0].setPrevious(null); // starting point
	93
	94	// the first column
	95	for (int j = 1; j < length2; j++) {
	96	matrix[0][j].setScore(0);
[1575]	97	//We don't need to go back to [0][0] if we reached matrix[0][x], so just end here
	98	//matrix[0][j].setPrevious(matrix[0][j-1]);
	99	matrix[0][j].setPrevious(null);
[1558]	100	}
	101
	102
	103
[1587]	104	for (int i = 1; i < length1 + 2; i++) {
[1558]	105
	106	// Formula for first row:
	107	// F(i,0) = max { F(i-1,0), F(i-1,j)-T j=1,...,m
	108
[1568]	109	double firstRowLeftScore = matrix[i-1][0].getScore();
[1559]	110	//for sequences of length 1
[1568]	111	double tempMax;
[1559]	112	int maxRowIndex;
	113	if(length2 == 1) {
	114	tempMax = matrix[i-1][0].getScore();
	115	maxRowIndex = 0;
	116	} else {
	117	tempMax = matrix[i-1][1].getScore();
	118	maxRowIndex = 1;
	119	//position of the maximal score of the previous row
	120
	121	for(int j = 2; j < length2;j++) {
	122	if(matrix[i-1][j].getScore() > tempMax) {
	123	tempMax = matrix[i-1][j].getScore();
	124	maxRowIndex = j;
	125	}
[1558]	126	}
[1559]	127
[1558]	128	}
[1559]	129
[1558]	130	tempMax -= scoreThreshold;
	131	matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
[1559]	132	if(tempMax ==matrix[i][0].getScore()){
	133	matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
	134	}
[1558]	135
[1559]	136	if(firstRowLeftScore == matrix[i][0].getScore()) {
	137	matrix[i][0].setPrevious(matrix[i-1][0]);
	138	}
[1558]	139
[1559]	140	//The last additional score is not related to a character in the input sequence, it's the total score. Therefore we don't need to save something for it
	141	if(i<length1+1)
	142	{
	143	matrix[i][0].setXvalue(input1[i-1]);
[1578]	144	matrix[i][0].setYvalue(Constants.UNMATCHED_SYMBOL);
[1559]	145	}
	146	else {
	147	//End after we calculated final score
	148	return;
	149	}
	150
	151
[1558]	152	for (int j = 1; j < length2; j++) {
[1568]	153	double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
	154	double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
	155	double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
[1558]	156
	157	matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
	158
	159	// find the directions that give the maximum scores.
[1578]	160	// TODO: Multiple directions are ignored, we choose the first maximum score
[1559]	161	//True if we had a match
[1558]	162	if (diagScore == matrix[i][j].getScore()) {
	163	matrix[i][j].setPrevious(matrix[i-1][j-1]);
[1559]	164	matrix[i][j].setXvalue(input1[i-1]);
	165	matrix[i][j].setYvalue(input2[j-1]);
[1558]	166	}
[1559]	167	//true if we took an event from sequence x and not from y
[1558]	168	if (leftScore == matrix[i][j].getScore()) {
[1559]	169	matrix[i][j].setXvalue(input1[i-1]);
[1578]	170	matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
[1558]	171	matrix[i][j].setPrevious(matrix[i-1][j]);
	172	}
[1559]	173	//true if we took an event from sequence y and not from x
[1558]	174	if (upScore == matrix[i][j].getScore()) {
[1578]	175	matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
[1559]	176	matrix[i][j].setYvalue(input2[j-1]);
[1558]	177	matrix[i][j].setPrevious(matrix[i][j-1]);
	178	}
[1559]	179	//true if we ended a matching region
[1558]	180	if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
[1559]	181	matrix[i][j].setPrevious(matrix[i][0]);
	182	matrix[i][j].setXvalue(input1[i-1]);
[1578]	183	matrix[i][j].setYvalue(Constants.UNMATCHED_SYMBOL);
[1558]	184	}
	185	}
[1559]	186
	187	//Set the complete score cell
	188
[1558]	189	}
	190	}
	191
	192	/**
	193	* Get the maximum value in the score matrix.
	194	*/
	195	public double getMaxScore() {
	196	double maxScore = 0;
	197
	198	// skip the first row and column
	199	for (int i = 1; i <= length1; i++) {
	200	for (int j = 1; j < length2; j++) {
	201	if (matrix[i][j].getScore() > maxScore) {
	202	maxScore = matrix[i][j].getScore();
	203	}
	204	}
	205	}
	206
	207	return maxScore;
	208	}
	209
[1586]	210	/* (non-Javadoc)
	211	* @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignmentScore()
[1558]	212	*/
[1586]	213	@Override
[1568]	214	public double getAlignmentScore() {
[1559]	215	return matrix[length1+1][0].getScore();
[1558]	216	}
	217
	218
	219
[1572]	220	public void traceback() {
[1559]	221	MatrixEntry tmp = matrix[length1+1][0];
[1587]	222	int aligned1[] = new int[length1+1+length2+2];
	223	int aligned2[] = new int[length1+1+length2+2];
[1572]	224	int count = 0;
	225	do
	226	{
	227	if(count != 0)
	228	{
[1587]	229	if (length1+1+length2+2 == count) {
[1581]	230	Console.traceln(Level.WARNING, "Traceback longer than both sequences summed up!");
	231	break;
	232	}
[1572]	233	aligned1[count] = tmp.getXvalue();
	234	aligned2[count] = tmp.getYvalue();
	235	}
	236
	237	tmp = tmp.getPrevious();
	238	count++;
[1581]	239
[1572]	240	} while(tmp != null);
[1575]	241	count--;
	242	//reverse order of the alignment
[1572]	243	int reversed1[] = new int[count];
	244	int reversed2[] = new int[count];
	245
[1574]	246
	247	for(int i = count-1; i > 0; i--) {
	248	reversed1[reversed1.length-i]= aligned1[i];
	249	reversed2[reversed2.length-i]= aligned2[i];
[1572]	250	}
	251
	252	NumberSequence ns1 = new NumberSequence(reversed1.length);
	253	NumberSequence ns2 = new NumberSequence(reversed2.length);
	254	ns1.setSequence(reversed1);
	255	ns2.setSequence(reversed2);
	256
	257	alignment.add(ns1);
	258	alignment.add(ns2);
	259	}
	260
	261	public void printAlignment() {
	262	MatrixEntry tmp = matrix[length1+1][0];
[1559]	263	String aligned1 = "";
	264	String aligned2 = "";
	265	int count = 0;
	266	do
	267	{
	268	String append1="";
	269	String append2="";
	270
[1578]	271	if(tmp.getXvalue() == Constants.GAP_SYMBOL) {
[1559]	272	append1 = " ___";
	273	}
[1578]	274	else if(tmp.getXvalue() == Constants.UNMATCHED_SYMBOL) {
[1559]	275	append1 = " ...";
	276	}
	277	else {
	278	append1 = String.format("%5d", tmp.getXvalue());
	279	}
[1558]	280
[1578]	281	if(tmp.getYvalue() == Constants.GAP_SYMBOL) {
[1559]	282	append2 = " ___";
	283	}
[1578]	284	else if(tmp.getYvalue() == Constants.UNMATCHED_SYMBOL) {
[1559]	285	append2 = " ...";
	286	}
	287	else {
	288	append2 = String.format("%5d", tmp.getYvalue());
	289	}
	290	if(count != 0)
	291	{
	292	aligned1 = append1 + aligned1;
	293	aligned2 = append2 + aligned2;
	294	}
	295
	296	tmp = tmp.getPrevious();
	297	count++;
	298
	299	} while(tmp != null);
	300	System.out.println(aligned1);
	301	System.out.println(aligned2);
	302	}
[1558]	303
[1559]	304
	305
[1558]	306	/**
	307	* print the dynmaic programming matrix
	308	*/
	309	public void printDPMatrix() {
	310	System.out.print(" ");
	311	for (int i = 1; i <= length1; i++)
	312	System.out.format("%5d", input1[i - 1]);
	313	System.out.println();
[1587]	314	for (int j = 0; j <= length2; j++) {
[1558]	315	if (j > 0)
	316	System.out.format("%5d ",input2[j - 1]);
	317	else{
	318	System.out.print(" ");
	319	}
[1559]	320	for (int i = 0; i <= length1 + 1; i++) {
	321	if((i<length1+1) \|\| (i==length1+1 && j==0)) {
	322	System.out.format("%4.1f ",matrix[i][j].getScore());
	323	}
	324
	325	}
[1558]	326	System.out.println();
	327	}
	328	}
	329
	330
[1586]	331	/* (non-Javadoc)
	332	* @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignment()
	333	*/
	334	@Override
[1585]	335	public ArrayList<NumberSequence> getAlignment() {
[1572]	336	return alignment;
	337	}
	338
[1585]	339	public void setAlignment(ArrayList<NumberSequence> alignment) {
[1572]	340	this.alignment = alignment;
	341	}
	342
[1558]	343	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: