📄 sim.c
字号:
/************************/
/* sim.c, 27 Oct 1994 */
/* runtime --> sim.res */
/************************/
/* A PROGRAM FOR LOCAL SIMILARITIES WITH AFFINE WEIGHTS:
Modified 12/6/90 for consistency with the "blast", "convert" and "lav"
programs.
copyright (c) 1990 Xiaoqiu Huang and Webb Miller
For the description of the algorithm, please see the paper
"A Time-Efficient, Linear-Space Local Similarity Algorithm"
(to appear in Advances in Applied Mathematics)
Xiaoqiu Huang and Webb Miller
Department of Computer Science
The Pennsylvania State University
University Park, PA 16802
An early version of this program is presented in the paper
"A Space-Efficient Algorithm for Local Similarities"
(to appear in Computer Applications in the Biosciences)
Xiaoqiu Huang, Ross C. Hardison and Webb Miller
Department of Computer Science
Department of Molecular and Cell Biology
The Pennsylvania State University
University Park, PA 16802
The SIM program finds k best non-intersecting alignments between
two sequences or within one sequence. Using dynamic programming
techniques, SIM is guaranteed to find optimal alignments. The
alignments are reported in order of similarity score, with the
highest scoring alignment first. The k best alignments share no
aligned pairs. SIM requires space proportional to the sum of the
input sequence lengths and the output alignment lengths. Thus
SIM can handle sequences of tens of thousands, even a hundred of
thousands, of base pairs on a workstation.
Users can supply values for the parameters:
M = cost of a matching aligned pair (default = 1)
I = cost of a transition (default is -1)
V = cost of a transversion (default is -1)
O = gap open penalty (default is 6.0)
E = gap extension penalty (default is 0.2)
Thus the score for an N-symbol indel is -(Q + R*N). Values are
specified with a command argument like O=5.5, where the given number
is rounded by SIM to the nearest tenth.
For example, to find 7 best non-intersecting alignments of segments
from two sequences in files A and B, and using the above default
values except that transistions are scored 0, use the command
sim 7 A B I=0
Acknowledgments
The functions diff() and display() are from Gene Myers. We made
the following modifications: similarity weights (integer), instead of
distance weights (float), are used, the aligned pairs already output
are not used in the subsequent computation, and the printed alignments
have been re-formatted in a number of small ways.
*/
#include <stdio.h>
static char *name1, *name2; /* names of sequence files */
#define CMAX 256 /* for EBCDIC */
/*#define round(x) nint(x) /* round to larger magnitude */
#define round(x) ((x)>0.0 ? (x)+0.5 : (x)-0.5)
main(argc, argv) int argc; char *argv[];
{ long M, N, K; /* Sequence lengths and k */
char *A, *B; /* Storing two sequences */
int symbol;
long V[CMAX][CMAX], Q, R; /* Converted integer weights */
double starttime, benchtime, dtime();
double parm_M, parm_I, parm_V, parm_O, parm_E, v;
double atof();
FILE *Bp, *Ap, *Cp, *ckopen();
char *ckalloc(); /* space-allocating function */
if ((Cp = fopen("sim.res","a+")) == NULL)
{
printf("Can not open sim.res\n\n");
exit(1);
}
starttime = dtime();
if ( argc < 4)
fatal("SIM k file1 file2 [M=] [I=] [V=] [O=] [E=]");
/* read k: the number of local alignments to find */
sscanf(argv[1],"%d", &K);
if (K == 0)
fatal("specified 0 alignments");
/* determine the sequence lengths */
Ap = ckopen(argv[2], "r");
for (M = 0; ( symbol = getc(Ap)) != EOF ; )
if ( symbol != '\n' )
++M;
fclose(Ap);
name1 = argv[2];
/* allocate space for A */
A = ( char * ) ckalloc( (M + 1) * sizeof(char));
/* read the first sequence into A */
Ap = ckopen(argv[2], "r");
for (M = 0; ( symbol = getc(Ap)) != EOF ; )
if ( symbol != '\n' )
A[++M] = symbol;
if (strcmp(argv[2],argv[3])) { /* sequences are different */
/* read the second sequence into B */
Bp = ckopen(argv[3], "r");
for (N = 0; ( symbol = getc(Bp)) != EOF ; )
if ( symbol != '\n' )
++N;
fclose(Bp);
name2 = argv[3];
B = ( char * ) ckalloc( (N + 1) * sizeof(char));
Bp = ckopen(argv[3], "r");
for (N = 0; ( symbol = getc(Bp)) != EOF ; )
if ( symbol != '\n' )
B[++N] = symbol;
}
parm_M = 1.0;
parm_I = -1.0;
parm_V = -1.0;
parm_O = 6.0;
parm_E = 0.2;
while (--argc > 3) {
if (argv[argc][1] != '=')
fatalf("argument %d has improper form", argc);
v = atof(argv[argc]+2);
switch (argv[argc][0]) {
case 'M': parm_M = v; break;
case 'I': parm_I = v; break;
case 'V': parm_V = v; break;
case 'O': parm_O = v; break;
case 'E': parm_E = v; break;
default: fatal("options are M, I, V, O and E.");
}
}
printf("\t\tSIM output with parameters:\n");
printf("\t\tM = %g, I = %g, V = %g\n\t\tO = %g, E = %g\n\n",
parm_M, parm_I, parm_V, parm_O, parm_E);
/* set up scoring matrix */
V['A']['A'] = V['C']['C'] =
V['G']['G'] = V['T']['T'] = round(10*parm_M);
V['A']['G'] = V['G']['A'] =
V['C']['T'] = V['T']['C'] = round(10*parm_I);
V['A']['C'] = V['A']['T'] =
V['C']['A'] = V['C']['G'] =
V['G']['C'] = V['G']['T'] =
V['T']['A'] = V['T']['G'] = round(10*parm_V);
Q = round (10*parm_O);
R = round (10*parm_E);
if (strcmp(argv[2], argv[3]))
SIM(A,B,M,N,K,V,Q,R,2L);
else
SIM(A,A,M,M,K,V,Q,R,1L);
benchtime = dtime() - starttime;
printf("\nRun Time (sec) = %9.1lf\n",benchtime);
fprintf(Cp,"\n");
fprintf(Cp," Run Line: sim 8 tob.38-44 liv.42-48\n");
fprintf(Cp," Run Time: %9.1lf (sec)\n",benchtime);
fprintf(Cp,"#######################################################\n");
fclose(Cp);
exit(0);
}
static long (*v)[CMAX]; /* substitution scores */
static long q, r; /* gap penalties */
static long qr; /* qr = q + r */
typedef struct ONE { long COL ; struct ONE *NEXT ;} pair, *pairptr;
pairptr *row, z; /* for saving used aligned pairs */
static short tt;
typedef struct NODE
{ long SCORE;
long STARI;
long STARJ;
long ENDI;
long ENDJ;
long TOP;
long BOT;
long LEFT;
long RIGHT; } vertex, *vertexptr;
vertexptr *LIST; /* an array for saving k best scores */
vertexptr low = 0; /* lowest score node in LIST */
vertexptr most = 0; /* latestly accessed node in LIST */
static long numnode; /* the number of nodes in LIST */
static long *CC, *DD; /* saving matrix scores */
static long *RR, *SS, *EE, *FF; /* saving start-points */
static long *HH, *WW; /* saving matrix scores */
static long *II, *JJ, *XX, *YY; /* saving start-points */
static long m1, mm, n1, nn; /* boundaries of recomputed area */
static long rl, cl; /* left and top boundaries */
static long min; /* minimum score in LIST */
static short flag; /* indicate if recomputation necessary*/
/* DIAG() assigns value to x if (ii,jj) is never used before */
#define DIAG(ii, jj, x, value) \
{ for ( tt = 1, z = row[(ii)]; z != 0; z = z->NEXT ) \
if ( z->COL == (jj) ) \
{ tt = 0; break; } \
if ( tt ) \
x = ( value ); \
}
/* replace (ss1, xx1, yy1) by (ss2, xx2, yy2) if the latter is large */
#define ORDER(ss1, xx1, yy1, ss2, xx2, yy2) \
{ if ( ss1 < ss2 ) \
{ ss1 = ss2; xx1 = xx2; yy1 = yy2; } \
else \
if ( ss1 == ss2 ) \
{ if ( xx1 < xx2 ) \
{ xx1 = xx2; yy1 = yy2; } \
else \
if ( xx1 == xx2 && yy1 < yy2 ) \
yy1 = yy2; \
} \
}
/* The following definitions are for function diff() */
long diff(), display();
static long zero = 0; /* long type zero */
#define gap(k) ((k) <= 0 ? 0 : q+r*(k)) /* k-symbol indel score */
static long *sapp; /* Current script append ptr */
static long last; /* Last script op appended */
static long I, J; /* current positions of A ,B */
static long no_mat; /* number of matches */
static long no_mis; /* number of mismatches */
static long al_len; /* length of alignment */
/* Append "Delete k" op */
#define DEL(k) \
{ I += k; \
al_len += k; \
if (last < 0) \
last = sapp[-1] -= (k); \
else \
last = *sapp++ = -(k); \
}
/* Append "Insert k" op */
#define INS(k) \
{ J += k; \
al_len += k; \
if (last < 0) \
{ sapp[-1] = (k); *sapp++ = last; } \
else \
last = *sapp++ = (k); \
}
/* Append "Replace" op */
#define REP \
{ last = *sapp++ = 0; \
al_len += 1; \
}
/* SIM(A,B,M,N,K,V,Q,R) reports K best non-intersecting alignments of
the segments of A and B in order of similarity scores, where
V[a][b] is the score of aligning a and b, and -(Q+R*i) is the score
of an i-symbol indel. */
SIM(A,B,M,N,K,V,Q,R,nseq) char A[],B[]; long M,N,K,V[][CMAX],Q,R,nseq;
{ long endi, endj, stari, starj; /* endpoint and startpoint */
long score; /* the max score in LIST */
long count; /* maximum size of list */
register long i, j; /* row and column indices */
char *ckalloc(); /* space-allocating function */
long *S; /* saving operations for diff */
vertexptr cur; /* temporary pointer */
vertexptr findmax(); /* return the largest score node */
/* allocate space for all vectors */
j = (N + 1) * sizeof(long);
CC = ( long * ) ckalloc(j);
DD = ( long * ) ckalloc(j);
RR = ( long * ) ckalloc(j);
SS = ( long * ) ckalloc(j);
EE = ( long * ) ckalloc(j);
FF = ( long * ) ckalloc(j);
i = (M + 1) * sizeof(long);
HH = ( long * ) ckalloc(i);
WW = ( long * ) ckalloc(i);
II = ( long * ) ckalloc(i);
JJ = ( long * ) ckalloc(i);
XX = ( long * ) ckalloc(i);
YY = ( long * ) ckalloc(i);
S = ( long * ) ckalloc(i + j);
row = ( pairptr * ) ckalloc( (M + 1) * sizeof(pairptr));
/* set up list for each row */
for ( i = 1; i <= M; i++ )
if ( nseq == 2 )
row[i] = 0;
else
{ row[i] = z = ( pairptr ) ckalloc( (long) sizeof(pair));
z->COL = i;
z->NEXT = 0;
}
v = V;
q = Q;
r = R;
qr = q + r;
LIST = ( vertexptr * ) ckalloc( K * sizeof(vertexptr));
for ( i = 0; i < K ; i++ )
LIST[i] = ( vertexptr ) ckalloc( (long) sizeof(vertex));
printf(" Upper Sequence : %s\n", name1);
printf(" Length : %d\n", M);
printf(" Lower Sequence : %s\n", name2);
printf(" Length : %d\n", N);
numnode = min = 0;
big_pass(A,B,M,N,K,nseq);
/* Report the K best alignments one by one. After each
alignment is output, recompute part of the matrix. First
determine the size of the area to be recomputed, then do the
recomputation */
for ( count = K - 1; count >= 0 ; count-- )
{ if ( numnode == 0 )
fatal("The number of alignments computed is too large");
cur = findmax();/* Return a pointer to a node with max score*/
score = cur->SCORE;
stari = ++cur->STARI;
starj = ++cur->STARJ;
endi = cur->ENDI;
endj = cur->ENDJ;
m1 = cur->TOP;
mm = cur->BOT;
n1 = cur->LEFT;
nn = cur->RIGHT;
rl = endi - stari + 1;
cl = endj - starj + 1;
I = stari - 1;
J = starj - 1;
sapp = S;
last = 0;
al_len = 0;
no_mat = 0;
no_mis = 0;
diff(&A[stari]-1, &B[starj]-1,rl,cl,q,q);
/* Output the best alignment */
printf("\n*********************************************************\n");
printf(" Number %d Local Alignment\n", K - count);
printf(" Similarity Score : %g\n",score/10.0);
printf(" Match Percentage : %d%%\n", (100*no_mat)/al_len);
printf(" Number of Matches : %d%\n", no_mat);
printf(" Number of Mismatches : %d%\n", no_mis);
printf(" Total Length of Gaps : %d%\n", al_len-no_mat-no_mis);
printf(" Begins at (%d, %d) and Ends at (%d, %d)\n",
stari,starj, endi,endj);
display(&A[stari]-1,&B[starj]-1,rl,cl,S,stari,starj);
fflush(stdout);
if ( count )
{ flag = 0;
locate(A,B,nseq);
if ( flag )
small_pass(A,B,count,nseq);
}
}
}
/* A big pass to compute K best classes */
big_pass(A,B,M,N,K,nseq) char A[],B[]; long M,N,K,nseq;
{ register long i, j; /* row and column indices */
register long c; /* best score at current point */
register long f; /* best score ending with insertion */
register long d; /* best score ending with deletion */
register long p; /* best score at (i-1, j-1) */
register long ci, cj; /* end-point associated with c */
register long di, dj; /* end-point associated with d */
register long fi, fj; /* end-point associated with f */
register long pi, pj; /* end-point associated with p */
long *va; /* pointer to v(A[i], B[j]) */
long addnode(); /* function for inserting a node */
/* Compute the matrix and save the top K best scores in LIST
CC : the scores of the current row
RR and EE : the starting point that leads to score CC
DD : the scores of the current row, ending with deletion
SS and FF : the starting point that leads to score DD */
/* Initialize the 0 th row */
for ( j = 1; j <= N ; j++ )
{ CC[j] = 0;
RR[j] = 0;
EE[j] = j;
DD[j] = - (q);
SS[j] = 0;
FF[j] = j;
}
for ( i = 1; i <= M; i++)
{ c = 0; /* Initialize column 0 */
f = - (q);
ci = fi = i;
va = v[A[i]];
if ( nseq == 2 )
{ p = 0;
pi = i - 1;
cj = fj = pj = 0;
}
else
{ p = CC[i];
pi = RR[i];
pj = EE[i];
cj = fj = i;
}
for ( j = (nseq == 2 ? 1 : (i+1)) ; j <= N ; j++ )
{ f = f - r;
c = c - qr;
ORDER(f, fi, fj, c, ci, cj)
c = CC[j] - qr;
ci = RR[j];
cj = EE[j];
d = DD[j] - r;
di = SS[j];
dj = FF[j];
ORDER(d, di, dj, c, ci, cj)
c = 0;
DIAG(i, j, c, p+va[B[j]]) /* diagonal */
if ( c <= 0 )
{ c = 0; ci = i; cj = j; }
else
{ ci = pi; cj = pj; }
ORDER(c, ci, cj, d, di, dj)
ORDER(c, ci, cj, f, fi, fj)
p = CC[j];
CC[j] = c;
pi = RR[j];
pj = EE[j];
RR[j] = ci;
EE[j] = cj;
DD[j] = d;
SS[j] = di;
FF[j] = dj;
if ( c > min ) /* add the score into list */
min = addnode(c, ci, cj, i, j, K, min);
}
}
}
/* Determine the left and top boundaries of the recomputed area */
locate(A,B,nseq) char A[],B[]; long nseq;
{ register long i, j; /* row and column indices */
register long c; /* best score at current point */
register long f; /* best score ending with insertion */
register long d; /* best score ending with deletion */
register long p; /* best score at (i-1, j-1) */
register long ci, cj; /* end-point associated with c */
register long di, dj; /* end-point associated with d */
register long fi, fj; /* end-point associated with f */
register long pi, pj; /* end-point associated with p */
short cflag, rflag; /* for recomputation */
long *va; /* pointer to v(A[i], B[j]) */
long addnode(); /* function for inserting a node */
long limit; /* the bound on j */
/* Reverse pass
rows
CC : the scores on the current row
RR and EE : the endpoints that lead to CC
DD : the deletion scores
SS and FF : the endpoints that lead to DD
columns
HH : the scores on the current columns
II and JJ : the endpoints that lead to HH
WW : the deletion scores
XX and YY : the endpoints that lead to WW
*/
for ( j = nn; j >= n1 ; j-- )
{ CC[j] = 0;
EE[j] = j;
DD[j] = - (q);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -