47790b284b066c70471a34a5a70d0277d88ee6ba tdreszer Fri May 6 10:47:51 2011 -0700 Another insignificant font color change that I forgot to check in yesterday. diff --git src/xenoAli/smithWat/smithWat.c src/xenoAli/smithWat/smithWat.c index 7e53112..720fdaf 100644 --- src/xenoAli/smithWat/smithWat.c +++ src/xenoAli/smithWat/smithWat.c @@ -1,268 +1,268 @@ /* dynAlign.c - align using dynamic programming */ #include "common.h" #include "cheapcgi.h" #include "htmshell.h" boolean isOnWeb; struct dynCell { struct dynCell *parent; int score; int qSlipScore; int tSlipScore; boolean visited; }; #define gapOpenCost (4*1024) #define matchBonus (1024) #define mismatchCost (1024) struct aliPair { struct aliPair *next; char querySym; char targetSym; }; void dumpMatrix(struct dynCell *matrix, int qDim, int tDim, struct dynCell *end) { int qIx, tIx; struct dynCell *cell, *parent; printf("<A HREF=#end>shortcut to end token</A>\n"); for (qIx = 0; qIx < qDim; ++qIx) { for (tIx = 0; tIx < tDim; ++tIx) { cell = &matrix[tIx*qDim + qIx]; parent = cell->parent; if (cell == end) - uglyf("<A NAME=end></A><FONT COLOR=#FF0000>"); + uglyf("<A NAME=end></A><span style='color:#FF0000;'>"); if (parent == NULL) uglyf("."); else if (parent == cell-1) uglyf("^"); else if (parent == cell-qDim) uglyf("<"); else uglyf("\\"); if (cell == end) - uglyf("</FONT>"); + uglyf("</span>"); } uglyf("\n"); } } void dynAlign(char *query, char *target, int matchMatrix[256][256], int gapFunc[], int gapFuncLen, FILE *f) /* Do dynamic programming to make alignment between query and target, and print it to stdout. */ { int querySize = strlen(query); int targetSize = strlen(target); int qDim = querySize+1; int tDim = targetSize+1; int matrixCellCount = qDim * tDim; size_t matrixBytes = matrixCellCount * sizeof(struct dynCell); struct dynCell *matrix = needLargeMem(matrixBytes); int matrixIx; int qIx, tIx, diag; int diagMax = querySize + targetSize + 1; int badScore = -0x6fffffff; int tSlipScore, qSlipScore, diagScore; int score; int endScore = -1; struct dynCell *parent, *cell, *newCell, *end = NULL, *start; struct aliPair *pairList = NULL, *pair; int qStart, tStart; int dif; /* Set first cell to zero */ zeroBytes(matrix, sizeof(matrix[0])); /* Do Smith-Waterman dynamic programming to find best local alignment. */ /* This combination of diag, qStart, tStart, qIx, and tIx loop control variables * is a little complex looking. All it's doing is letting us go over the rectangle * from 0 to querySize by 0 to target diagonally. If querySize was 5 and targetSize 2 * we'd traverse the matrix like so: * 0 1 2 3 4 5 6 * 1 2 3 4 5 6 7 * 2 3 4 5 6 7 8 */ for (diag = 1; diag < diagMax; ++diag) { dif = diag - querySize; if (dif > 0) { qStart = querySize; tStart = dif; } else { qStart = diag; tStart = 0; } for (qIx = qStart, tIx = tStart; qIx >= 0 && tIx < tDim; --qIx, ++tIx) { newCell = matrix + tIx*qDim + qIx; if (qIx > 0 && tIx > 0) { /* Initially set things to go as if query and target match at this position. */ parent = cell = newCell - qDim - 1; score = diagScore = cell->score + matchMatrix[query[qIx-1]][target[tIx-1]]; /* Look what happens if insert a gap in query. */ cell = newCell - qDim; qSlipScore = cell->score - gapOpenCost; if (qSlipScore > score) { score = qSlipScore; parent = cell; } /* Look what happens if insert a gap in target. */ cell = newCell-1; tSlipScore = cell->score - gapOpenCost; if (tSlipScore > score) { score = tSlipScore; parent = cell; } /* If can't raise positive score start over. */ if (score <= 0) { parent = NULL; score = 0; } /* Keep track of best score. */ if (score > endScore) { endScore = score; end = newCell; } newCell->parent = parent; newCell->score = score; } } } start = matrix; /* Find best scoring position in matrix to use as start. */ parent = end; for (;;) { cell = parent; if ((parent = cell->parent) == NULL) break; AllocVar(pair); matrixIx = cell-matrix; tIx = matrixIx/qDim; qIx = matrixIx%qDim; if (parent == cell - qDim - 1) /* Diagonal step */ { pair->querySym = query[qIx-1]; pair->targetSym = target[tIx-1]; } else if (parent == cell-1) { pair->querySym = query[qIx-1]; pair->targetSym = '-'; } else { assert(parent == cell-qDim); pair->querySym = '-'; pair->targetSym = target[tIx-1]; } slAddHead(&pairList, pair); } printf("%s\n%s\n------------------------------------------------\n", query, target); for (pair = pairList; pair != NULL; pair = pair->next) putc(pair->querySym, f); putc('\n', f); for (pair = pairList; pair != NULL; pair = pair->next) { char c = (pair->targetSym == pair->querySym ? '|' : ' '); putc(c, f); } putc('\n', f); for (pair = pairList; pair != NULL; pair = pair->next) putc(pair->targetSym, f); putc('\n', f); printf("------------------------------------------------\n"); dumpMatrix(matrix, qDim, tDim, end); } static int matchMatrix[256][256]; int gapFuncTable[1024]; void makeGapFuncTable() { int i; int decay = mismatchCost; int gapCost = gapOpenCost; gapFuncTable[0] = 0; gapFuncTable[1] = gapOpenCost; for (i=2; i<ArraySize(gapFuncTable); ++i) { gapFuncTable[i] = decay; if ((decay >>= 1) == 0) decay = 1; } } void dustData(char *in,char *out) /* Clean up data by getting rid non-alphabet stuff. */ { char c; while ((c = *in++) != 0) { if (isalpha(c)) *out++ = toupper(c); } *out++ = 0; } void doMiddle() { char *query = cgiString("query"); char *target = cgiString("target"); dustData(query, query); dustData(target, target); printf("<PRE><TT>"); dynAlign(query, target, matchMatrix, gapFuncTable, ArraySize(gapFuncTable), stdout); } int main(int argc, char *argv[]) { char *query, *target; int i,j; for (i = 0; i<256; ++i) for (j=0; j<256; ++j) { matchMatrix[i][j] = -mismatchCost; } for (i=0; i<256; ++i) matchMatrix[i][i] = matchBonus; makeGapFuncTable(); if ((isOnWeb = cgiIsOnWeb()) == TRUE) { htmShell("DynaAlign Results", doMiddle, NULL); } else { if (argc != 3) { errAbort("dynAlign - used dynamic programming to find alignment between two strings\n" "usage:\n" " dynAlign string1 string2"); } query = argv[1]; target = argv[2]; dynAlign(query, target, matchMatrix, gapFuncTable, ArraySize(gapFuncTable), stdout); } return 0; }