--- /dev/null
+/*
+ Fast, but crude, pairwise structural Alignments of RNA sequences
+
+ Possible structures of each RNA are encoded in a linear
+ "probability profile", by computing for each base the probability
+ of being unpaired, or paired upstream or downstream. These profiles
+ can be aligned using standard string alignment.
+
+ The is an extension of the old method in ProfileDist.c with the
+ following changes:
+ - use sequence as well as structure profile for scoring
+ - use similarity alignment instead of distance (maybe add local alinment)
+ - use affine gap costs
+
+ C Ivo L Hofacker, Vienna RNA Package
+*/
+/* Last changed Time-stamp: <2004-08-02 10:11:13 ivo> */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+#include <float.h>
+#include "dist_vars.h"
+#include "fold_vars.h"
+#include "part_func.h"
+#include "utils.h"
+#include "profiledist.h"
+#include "ProfileAln.h"
+
+
+/*@unused@*/
+static char rcsid[] = "$Id: ProfileAln.c,v 1.5 2006/01/18 13:00:30 ivo Exp $";
+
+#define EQUAL(x,y) (fabs((x)-(y)) <= fabs(x)*2*FLT_EPSILON)
+
+PRIVATE int *alignment[2];
+
+PRIVATE void sprint_aligned_bppm(const float *T1, const char *seq1,
+ const float *T2, const char *seq2);
+PRIVATE double PrfEditScore(const float *p1, const float *p2,
+ char c1, char c2);
+PRIVATE double average(double x, double y);
+
+PRIVATE double open=-1.5, ext=-0.666; /* defaults from clustalw */
+PRIVATE double seqw=0.5;
+PRIVATE int free_ends=1; /* whether to use free end gaps */
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE float **newmat(int l1, int l2) {
+ float **a;
+ int i;
+ a = (float **) space((l1+1)*sizeof(float *));
+ for (i=0; i<=l1; i++) a[i] = (float *) space((l2+1)*sizeof(float));
+ return a;
+}
+
+PUBLIC float profile_aln(const float *T1, const char *seq1,
+ const float *T2, const char *seq2)
+{
+ /* align the 2 probability profiles T1, T2 */
+ /* This is like a Needleman-Wunsch alignment, with affine gap-costs
+ ala Gotoh. The score looks at both seq and pair profile */
+
+ float **S, **E, **F, tot_score;
+ int i, j, length1, length2;
+
+ length1 = strlen(seq1);
+ length2 = strlen(seq2);
+ S = newmat(length1, length2);
+ E = newmat(length1, length2);
+ F = newmat(length1, length2);
+
+ E[0][0] = F[0][0] = open - ext;
+ S[0][0] = 0;
+ for (i=1; i<=length1; i++) F[i][0] = -9999; /* impossible */
+ for (j=1; j<=length2; j++) E[0][j] = -9999; /* impossible */
+ if (!free_ends) {
+ for (i=1; i<=length1; i++) S[i][0] = E[i][0] = E[i-1][0] +ext;
+ for (j=1; j<=length2; j++) S[0][j] = F[0][j] = F[0][j-1] +ext;
+ }
+
+ for (i=1; i<=length1; i++) {
+ for (j=1; j<=length2; j++) {
+ float M;
+ E[i][j] = MAX2(E[i-1][j]+ext, S[i-1][j]+open);
+ F[i][j] = MAX2(F[i][j-1]+ext, S[i][j-1]+open);
+ M = S[i-1][j-1] + PrfEditScore(T1+3*i,T2+3*j, seq1[i-1], seq2[j-1]);
+ S[i][j] = MAX3(M, E[i][j], F[i][j]);
+ }
+ }
+
+ if (edit_backtrack) {
+ double score=0;
+ char state = 'S';
+ int pos, i,j;
+ alignment[0] = (int *) space((length1+length2+1)*sizeof(int));
+ alignment[1] = (int *) space((length1+length2+1)*sizeof(int));
+
+ pos = length1+length2;
+ i = length1;
+ j = length2;
+
+ tot_score = S[length1][length2];
+
+ if (free_ends) {
+ /* find starting point for backtracking,
+ search for highest entry in last row or column */
+ int imax=0;
+ for (i=1; i<=length1; i++) {
+ if (S[i][length2]>score) {
+ score=S[i][length2];
+ imax=i;
+ }
+ }
+ for (j=1; j<=length2; j++) {
+ if (S[length1][j]>score) {
+ score=S[length1][j];
+ imax=-j;
+ }
+ }
+ if (imax<0) {
+ for (j=length2; j> -imax; j--) {
+ alignment[0][pos] = 0;
+ alignment[1][pos--] = j;
+ }
+ i=length1;
+ } else {
+ for (i=length1; i>imax; i--) {
+ alignment[0][pos] = i;
+ alignment[1][pos--] = 0;
+ }
+ j=length2;
+ }
+ tot_score=score;
+ }
+
+ while (i>0 && j>0) {
+ switch (state) {
+ case 'E':
+ score = E[i][j];
+ alignment[0][pos] = i;
+ alignment[1][pos--] = 0;
+ if (EQUAL(score, S[i-1][j] + open)) state = 'S';
+ i--;
+ break;
+ case 'F':
+ score = F[i][j];
+ alignment[0][pos] = 0;
+ alignment[1][pos--] = j;
+ if (EQUAL(score, S[i][j-1] + open)) state = 'S';
+ j--;
+ break;
+ case 'S':
+ score = S[i][j];
+ if (EQUAL(score, E[i][j])) state = 'E';
+ else if (EQUAL(score, F[i][j])) state = 'F';
+ else if (EQUAL(score, S[i-1][j-1] +
+ PrfEditScore(T1+3*i,T2+3*j, seq1[i-1], seq2[j-1]))) {
+ alignment[0][pos] = i;
+ alignment[1][pos--] = j;
+ i--; j--;
+ }
+ else nrerror("backtrack of alignment failed");
+ break;
+ }
+ }
+
+ for (; j>0; j--) {
+ alignment[0][pos] = 0;
+ alignment[1][pos--] = j;
+ }
+ for (; i>0; i--) {
+ alignment[0][pos] = i;
+ alignment[1][pos--] = 0;
+ }
+
+ for(i=pos+1; i<=length1+length2; i++){
+ alignment[0][i-pos] = alignment[0][i];
+ alignment[1][i-pos] = alignment[1][i];
+ }
+ alignment[0][0] = length1+length2-pos; /* length of alignment */
+
+ sprint_aligned_bppm(T1,seq1, T2,seq2);
+ free(alignment[0]);
+ free(alignment[1]);
+ }
+ for (i=0; i<=length1; i++) {
+ free(S[i]); free(E[i]); free(F[i]);
+ }
+ free(S); free(E); free(F);
+
+ return tot_score;
+}
+
+
+/*---------------------------------------------------------------------------*/
+PRIVATE inline double average(double x, double y) {
+ /*
+ As in Bonhoeffer et al (1993) 'RNA Multi Structure Landscapes',
+ Eur. Biophys. J. 22: 13-24 we have chosen the geometric mean.
+ */
+ return (float) sqrt(x*y);
+}
+
+PRIVATE double PrfEditScore(const float *p1, const float *p2, char c1, char c2)
+{
+ double score;
+ int k;
+
+ for(score=0.,k=0; k<3; k++)
+ score += average(p1[k],p2[k]);
+
+ score *= (1- seqw);
+ if (c1==c2) score += seqw;
+ else if (((c1=='A') && (c2=='G')) ||
+ ((c1=='G') && (c2=='A')) ||
+ ((c1=='C') && (c2=='U')) ||
+ ((c1=='U') && (c2=='C')))
+ score += 0.5*seqw;
+ else score -= 0.9*seqw;
+ return score;
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void sprint_aligned_bppm(const float *T1, const char *seq1,
+ const float *T2, const char *seq2) {
+ int i, length;
+ length = alignment[0][0];
+ for (i=0; i<4; i++) {
+ if (aligned_line[i] != NULL) free(aligned_line[i]);
+ aligned_line[i] = (char *) space((length+1)*sizeof(char));
+ }
+ for(i=1; i<=length; i++){
+ if (alignment[0][i]==0)
+ aligned_line[0][i-1] = aligned_line[2][i-1] = '_';
+ else {
+ aligned_line[0][i-1] = bppm_symbol(T1+alignment[0][i]*3);
+ aligned_line[2][i-1] = seq1[alignment[0][i]-1];
+ }
+ if (alignment[1][i]==0)
+ aligned_line[1][i-1] = aligned_line[3][i-1] = '_';
+ else {
+ aligned_line[1][i-1] = bppm_symbol(T2+alignment[1][i]*3);
+ aligned_line[3][i-1] = seq2[alignment[1][i]-1];
+ }
+ }
+}
+
+PUBLIC int set_paln_params(double gap_open, double gap_ext,
+ double seq_weight, int freeends) {
+ open = (gap_open>0) ? -gap_open : gap_open;
+ ext = (gap_ext>0) ? -gap_ext : gap_ext;
+ if (open > ext) fprintf(stderr, "Gap extension penalty is smaller than "
+ "gap open. Do you realy want this?\n");
+ seqw = seq_weight;
+ if (seqw<0) {
+ seqw = 0;
+ fprintf(stderr, "Sequence weight set to 0 (must be in [0..1])\n");
+ } else
+ if (seqw>1) {
+ seqw = 1;
+ fprintf(stderr, "Sequence weight set to 1 (must be in [0..1])\n");
+ }
+ free_ends = (freeends) ? 1 : 0;
+ return 0;
+}
+
+/*---------------------------------------------------------------------------*/
git://source.jalview.org / jabaws.git / blobdiff