--- /dev/null
+/*
+ Distances of Secondary Structures
+ Walter Fontana, Ivo L Hofacker, Peter F Stadler
+ Vienna RNA Package
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <unistd.h>
+#include <string.h>
+#include "dist_vars.h"
+#include "RNAstruct.h"
+#include "treedist.h"
+#include "stringdist.h"
+#include "utils.h"
+#include "data_structures.h"
+#include "RNAdistance_cmdl.h"
+
+#define MAXNUM 1000 /* max number of structs for distance matrix */
+
+#define PUBLIC
+#define PRIVATE static
+/*@unused@*/
+static char rcsid[] = "$Id: RNAdistance.c,v 1.8 2005/07/24 08:35:15 ivo Exp $";
+PRIVATE void command_line(int argc, char *argv[]);
+PRIVATE int parse_input(char *line);
+PRIVATE int check_tree(char *line, char alpha[]);
+PRIVATE int check_brackets(char *line);
+PRIVATE void print_aligned_lines(FILE *somewhere);
+
+PRIVATE char ruler[] ="....,....1....,....2....,....3....,....4"
+"....,....5....,....6....,....7....,....8";
+PRIVATE int types=1;
+PRIVATE int task;
+PRIVATE int taxa_list;
+PRIVATE char outfile[FILENAME_MAX_LENGTH], *list_title;
+
+PRIVATE char ttype[10]="f";
+PRIVATE int n=0;
+
+int main(int argc, char *argv[])
+{
+ char *line=NULL, *xstruc, *cc;
+ Tree *T[10][MAXNUM];
+ int tree_types = 0, ttree;
+ swString *S[10][MAXNUM];
+ char *P[MAXNUM]; /* structures for base pair distances */
+ int string_types = 0, tstr;
+ int i,j, tt, istty, type;
+ int it, is;
+ FILE *somewhere=NULL;
+
+ command_line(argc, argv);
+
+ if((outfile[0]=='\0')&&(task==2)&&(edit_backtrack))
+ strcpy(outfile,"backtrack.file");
+ if(outfile[0]!='\0') somewhere = fopen(outfile,"w");
+ if (somewhere==NULL) somewhere = stdout;
+ istty = isatty(fileno(stdin))&&isatty(fileno(stdout));
+
+ do {
+ if ((istty)&&(n==0)) {
+ printf("\nInput structure; @ to quit\n");
+ printf("%s\n", ruler);
+ }
+ do {
+ if (line!=NULL) free(line);
+ line=get_line(stdin);
+ } while ((type=parse_input(line))==0);
+
+ if (((type==999)||(type==888))&&(task==2)) { /* do matrices */
+ if (taxa_list)
+ printf("* END of taxa list\n");
+
+ ttree = 0; tstr = 0;
+ for (tt=0; tt< types; tt++) {
+ printf("> %c %d\n", ttype[tt], n);
+ if(islower(ttype[tt])) {
+ for (i=1; i<n; i++) {
+ for (j=0; j<i; j++) {
+ printf("%g ",tree_edit_distance(T[ttree][i], T[ttree][j]));
+ if(edit_backtrack) {
+ fprintf(somewhere,"%d %d",i+1,j+1);
+ if (ttype[tt]=='f') unexpand_aligned_F(aligned_line);
+ print_aligned_lines(somewhere);
+ }
+ }
+ printf("\n");
+ }
+ printf("\n");
+ for (i=0; i<n; i++) free_tree(T[ttree][i]);
+ ttree++;
+ }
+ if (ttype[tt]=='P') {
+ for (i=1; i<n; i++) {
+ for (j=0; j<i; j++)
+ printf("%g ", (float) bp_distance(P[i], P[j]));
+ printf("\n");
+ }
+ printf("\n");
+ for (i=0; i<n; i++) free(P[i]);
+ }
+ else if(isupper(ttype[tt])) {
+ for (i=1; i<n; i++) {
+ for (j=0; j<i; j++) {
+ printf("%g ",string_edit_distance(S[tstr][i], S[tstr][j]));
+ if (edit_backtrack) {
+ fprintf(somewhere,"%d %d",i+1,j+1);
+ if (ttype[tt]=='F') unexpand_aligned_F(aligned_line);
+ print_aligned_lines(somewhere);
+ }
+ }
+ printf("\n");
+ }
+ printf("\n");
+ for (i=0; i<n; i++) free(S[tstr][i]);
+ tstr++;
+ }
+ }
+ fflush(stdout);
+ if (type==888) { /* do another distance matrix */
+ n = 0;
+ printf("%s\n", list_title);
+ free(list_title);
+ continue;
+ }
+ } /* END do matrices */
+ if (type==999) { /* finito */
+ if (outfile[0]!='\0') fclose(somewhere);
+ return 0;
+ }
+
+ if (type<0) {
+ xstruc = add_root(line);
+ free(line);
+ line = xstruc;
+ type = -type;
+ }
+ if (type==2) {
+ xstruc = unexpand_Full(line);
+ free(line);
+ line = xstruc;
+ type=1;
+ }
+ tree_types = 0;
+ string_types = 0;
+ for(tt=0; tt < types; tt++) {
+ switch(ttype[tt]){
+ case 'f' :
+ case 'F' :
+ if (type!=1) nrerror("Can't convert back to full structure");
+ xstruc = expand_Full(line);
+ if(islower(ttype[tt])) { /* tree_edit */
+ T[tree_types++][n] = make_tree(xstruc);
+ }
+ if(isupper(ttype[tt])) { /* string edit */
+ S[string_types++][n] = Make_swString(xstruc);
+ }
+ free(xstruc);
+ break;
+ case 'P':
+ if (type!=1) nrerror("Can't convert back to full structure");
+ P[n] = strdup(line);
+ break;
+ case 'h' :
+ case 'H' :
+ switch (type) {
+ case 1:
+ xstruc = b2HIT(line);
+ if(islower(ttype[tt])) {
+ T[tree_types++][n] = make_tree(xstruc);
+ }
+ if(isupper(ttype[tt])) {
+ S[string_types++][n] = Make_swString(xstruc);
+ }
+ free(xstruc);
+ break;
+ default:
+ nrerror("Can't convert to HIT structure");
+ }
+ break;
+ case 'c' :
+ case 'C' :
+ switch (type) {
+ case 1:
+ cc = b2C(line);
+ xstruc = expand_Shapiro(cc);
+ free(cc);
+ break;
+ case 4:
+ cc = expand_Shapiro(line);
+ xstruc = unweight(cc);
+ free(cc);
+ break;
+ case 3:
+ xstruc = unweight(line);
+ break;
+ default:
+ nrerror("Unknown structure representation");
+ exit(0);
+ }
+ if(islower(ttype[tt])) {
+ T[tree_types++][n] = make_tree(xstruc);
+ }
+ if(isupper(ttype[tt])) {
+ S[string_types++][n] = Make_swString(xstruc);
+ }
+ free(xstruc);
+ break;
+ case 'w' :
+ case 'W' :
+ if (type==1) {
+ xstruc = b2Shapiro(line);
+ if(islower(ttype[tt])) {
+ T[tree_types++][n] = make_tree(xstruc);
+ }
+ if(isupper(ttype[tt])) {
+ S[string_types++][n] = Make_swString(xstruc);
+ }
+ free(xstruc);
+ }
+ else {
+ if(islower(ttype[tt])) {
+ T[tree_types++][n] = make_tree(line);
+ }
+ if(isupper(ttype[tt])) {
+ S[string_types++][n] = Make_swString(line);
+ }
+ }
+ break;
+ default:
+ nrerror("Unknown distance type");
+ }
+ }
+ n++;
+ switch (task) {
+ float dist;
+ case 1:
+ if (n==2) {
+ for (it=0, is=0, i=0; i<types; i++) {
+ if(islower(ttype[i])) {
+ dist = tree_edit_distance(T[it][0], T[it][1]);
+ free_tree(T[it][0]); free_tree(T[it][1]);
+ it++;
+ }
+ else if (ttype[i]=='P') {
+ dist = (float) bp_distance(P[0], P[1]);
+ free(P[0]); free(P[1]);
+ }
+ else /* isupper(ttype[i]) */ {
+ dist = string_edit_distance(S[is][0], S[is][1]);
+ free(S[is][0]); free(S[is][1]);
+ is++;
+ }
+ printf("%c: %g ", ttype[i], dist);
+ if ((edit_backtrack)&&(ttype[i]!='P')) {
+ if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
+ print_aligned_lines(somewhere);
+ }
+ }
+ printf("\n");
+ n=0;
+ }
+ break;
+ case 3:
+ if (n>1) {
+ for (it=0, is=0, i=0; i<types; i++) {
+ if(islower(ttype[i])) {
+ dist = tree_edit_distance(T[it][1], T[it][0]);
+ free_tree(T[it][1]);
+ it++;
+ }
+ else if (ttype[i]=='P') {
+ dist = (float) bp_distance(P[0], P[1]);
+ free(P[1]);
+ }
+ else /* if(isupper(ttype[i])) */ {
+ dist = string_edit_distance(S[is][0], S[is][1]);
+ free(S[is][1]);
+ is++;
+ }
+ printf("%c: %g ", ttype[i], dist);
+ if ((edit_backtrack)&&(ttype[i]!='P')) {
+ if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
+ print_aligned_lines(somewhere);
+ }
+ }
+ printf("\n");
+ n=1;
+ }
+ break;
+ case 4:
+ if (n>1) {
+ for (it=0, is=0, i=0; i<types; i++) {
+ if(islower(ttype[i])) {
+ dist = tree_edit_distance(T[it][1], T[it][0]);
+ free_tree(T[it][0]);
+ T[it][0] = T[it][1];
+ it++;
+ }
+ else if (ttype[i]=='P') {
+ dist = (float) bp_distance(P[0], P[1]);
+ free(P[0]); P[0] = P[1];
+ }
+ else /* if(isupper(ttype[i])) */ {
+ dist = string_edit_distance(S[is][0], S[is][1]);
+ free(S[is][0]);
+ S[is][0] = S[is][1];
+ is++;
+ }
+ printf("%c: %g ", ttype[i], dist);
+ if ((edit_backtrack)&&(ttype[i]!='P')) {
+ if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
+ print_aligned_lines(somewhere);
+ }
+ }
+ printf("\n");
+ n=1;
+ }
+ break;
+ }
+ fflush(stdout);
+ } while(type!=999);
+ return 0;
+}
+
+/*--------------------------------------------------------------------------*/
+
+PRIVATE int parse_input(char *line)
+{
+ int type, rooted=0, i, xx;
+ char *cp;
+
+ if (line==NULL) return 999;
+ if (line[0]=='*') {
+ if (taxa_list==0) {
+ if (task==2) taxa_list=1;
+ printf("%s\n", line);
+ return 0;
+ } else {
+ list_title = strdup(line);
+ return 888;
+ }
+ }
+ if (line[0]=='>') {
+ if (taxa_list)
+ printf("%d :%s\n", n+1, line+1);
+ else printf("%s\n", line);
+ return 0;
+ }
+
+ cp = strchr(line,' ');
+ if (cp) *cp='\0'; /* get rid of junk at the end of line */
+
+ switch (line[0]) {
+ case '.' :
+ type = 1;
+ break;
+
+ case '(' : /* it's a tree */
+ i=1;
+ xx = 0;
+ type = 4; /* coarse */
+ rooted = 0;
+ while (line[i]) {
+ if (line[i]=='.'){
+ type = 1; /* Full */
+ break;
+ }
+ if ( (line[i]=='U')||(line[i]=='P') ) {
+ type = 2; /* FFull */
+ xx = 1;
+ break;
+ }
+ if (line[i]=='S') {
+ type = 3;
+ xx = 1;
+ break; /* Shapiro tree */
+ }
+ if ((line[i]!='(')&&(line[i]!=')')) xx = 1;
+ i++;
+ }
+ if(!xx) type =1;
+
+ rooted = (line[strlen(line)-2]=='R');
+ break;
+ case '@' :
+ return 999;
+
+ default:
+ return 0;
+ }
+
+ switch (type) {
+ case 1:
+ if(check_brackets(line))
+ return 1;
+ break;
+ case 2:
+ if(check_tree(line,"UP") ) {
+ if(rooted) return 2;
+ else return -2;
+ }
+ break;
+ case 3:
+ if(check_tree(line,"HBIMSE") ){
+ if(rooted) return -3;
+ else return -3;
+ }
+ break;
+ case 4:
+ if(check_tree(line,"HBIM") ){
+ if(rooted) return 4;
+ else return -4;
+ }
+ break;
+ }
+ return 0;
+}
+
+/*--------------------------------------------------------------------------*/
+
+
+PRIVATE int check_tree(char *line, char alpha[])
+{
+ int n, i, o;
+ char *pos;
+
+ n = (int) strlen(line);
+
+ if( line[0] != '(' ) return 0;
+ i=o=1;
+ while( line[i] ){
+ while( line[i] == '(' ){
+ o++;
+ i++;
+ }
+ pos=strchr(alpha, (int)line[i]);
+ if (pos) {
+ while (isdigit((int) line[++i]));
+ if (line[i]!=')') return 0;
+ }
+ if (line[i]=='R') {
+ i++;
+ if ((i!=n-1)||(line[i]!=')')) return 0;
+ }
+ if (line[i]==')') {
+ o--;
+ if(o< 0) return 0;
+ if( (i<n)&&(line[i+1]==')') ) return 0;
+ i++;
+ }
+ else return 0;
+ }
+ if(o>0) return 0;
+ return 1;
+}
+
+/*--------------------------------------------------------------------------*/
+
+
+PRIVATE int check_brackets(char *line)
+{
+ int i,o;
+
+ i=o=0;
+ while( line[i] ){
+ switch(line[i]) {
+ case '(' :
+ o++;
+ i++;
+ break;
+ case '.' :
+ i++;
+ break;
+ case ')' :
+ i++;
+ o--;
+ if(o<0) return 0;
+ break;
+ default:
+ return 0;
+ }
+ }
+ if (o>0) return 0;
+ return 1;
+}
+
+/*--------------------------------------------------------------------------*/
+
+
+PRIVATE void command_line(int argc, char *argv[])
+{
+ int i;
+ struct RNAdistance_args_info args_info;
+
+ /* default values */
+ edit_backtrack = 0;
+ types=1; ttype[0]='f'; task=1;
+
+ /*
+ #############################################
+ # check the command line parameters
+ #############################################
+ */
+ if(RNAdistance_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
+
+ /* use specified distance representations */
+ if(args_info.distance_given){
+ strncpy(ttype, args_info.distance_arg,9);
+ types = (int)strlen(ttype);
+ }
+
+ if(args_info.compare_given){
+ switch(args_info.compare_arg[0]){
+ case 'p': task=1; break;
+ case 'm': task=2; break;
+ case 'f': task=3; break;
+ case 'c': task=4; break;
+ default: RNAdistance_cmdline_parser_print_help();
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ if(args_info.shapiro_given)
+ cost_matrix = 1;
+
+ if(args_info.backtrack_given){
+ if(strcmp(args_info.backtrack_arg, "none")){
+ strncpy(outfile, args_info.backtrack_arg, FILENAME_MAX_LENGTH-1);
+ }
+ edit_backtrack = 1;
+ }
+
+ /* free allocated memory of command line data structure */
+ RNAdistance_cmdline_parser_free (&args_info);
+}
+
+/*--------------------------------------------------------------------------*/
+
+PRIVATE void print_aligned_lines(FILE *somewhere)
+{
+ if (edit_backtrack) {
+ fprintf(somewhere, "\n%s\n%s\n", aligned_line[0], aligned_line[1]);
+ fflush(somewhere);
+ }
+}
git://source.jalview.org / jabaws.git / blobdiff