binaries/src/ViennaRNA/lib/aln_util.c

   1 /*
   2                                aln_util.c
   3                Helper functions frelated to alignments
   4 */
   5 /* Last changed Time-stamp: <2006-01-16 11:42:38 ivo> */
   6
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <errno.h>
  10 #include <time.h>
  11 #include <string.h>
  12 #include <ctype.h>
  13 #include <config.h>
  14 #include "utils.h"
  15 #include "fold_vars.h"
  16 #include "pair_mat.h"
  17 /*@unused@*/
  18 static char rcsid[] = "$Id: aln_util.c,v 1.4 2007/02/02 15:18:13 ivo Exp $";
  19
  20 #define MAX_NUM_NAMES    500
  21 int read_clustal(FILE *clust, char *AlignedSeqs[], char *names[]) {
  22    char *line, name[100]="", *seq;
  23    int  n, nn=0, num_seq = 0, i;
  24
  25    if ((line=get_line(clust)) == NULL) {
  26      fprintf(stderr, "Empty CLUSTAL file\n"); return 0;
  27    }
  28
  29    if ((strncmp(line,"CLUSTAL", 7) !=0) && (!strstr(line,"STOCKHOLM"))) {
  30      fprintf(stderr, "This doesn't look like a CLUSTAL/STOCKHOLM file, sorry\n");
  31      free(line); return 0;
  32    }
  33    free(line);
  34    line = get_line(clust);
  35
  36    while (line!=NULL) {
  37     if(strncmp(line, "//", 2) == 0){
  38       free(line);
  39       break;
  40     }
  41
  42     if (((n=strlen(line))<4) || isspace((int)line[0])) {
  43       /* skip non-sequence line */
  44       free(line); line = get_line(clust);
  45       nn=0; /* reset seqence number */
  46       continue;
  47     }
  48     /* skip comments */
  49     if(line[0] == '#'){
  50       free(line);
  51       line = get_line(clust);
  52       continue;
  53     }
  54
  55      seq = (char *) space( (n+1)*sizeof(char) );
  56      sscanf(line,"%99s %s", name, seq);
  57
  58     for(i=0;i<strlen(seq);i++){
  59       if(seq[i] == '.') seq[i] = '-'; /* replace '.' gaps by '-' */
  60       /* comment the next line and think about something more difficult to deal with
  61          lowercase sequence letters if you really want to */
  62       seq[i] = toupper(seq[i]);
  63     }
  64
  65      if (nn == num_seq) { /* first time */
  66        names[nn] = strdup(name);
  67        AlignedSeqs[nn] = strdup(seq);
  68      }
  69      else {
  70        if (strcmp(name, names[nn])!=0) {
  71          /* name doesn't match */
  72          fprintf(stderr,
  73                  "Sorry, your file is fucked up (inconsitent seq-names)\n");
  74          free(line); free(seq);
  75          return 0;
  76        }
  77        AlignedSeqs[nn] = (char *)
  78          xrealloc(AlignedSeqs[nn], strlen(seq)+strlen(AlignedSeqs[nn])+1);
  79        strcat(AlignedSeqs[nn], seq);
  80      }
  81      nn++;
  82      if (nn>num_seq) num_seq = nn;
  83      free(seq);
  84      free(line);
  85      if (num_seq>=MAX_NUM_NAMES) {
  86        fprintf(stderr, "Too many sequences in CLUSTAL/STOCKHOLM file");
  87        return 0;
  88      }
  89
  90      line = get_line(clust);
  91    }
  92
  93    AlignedSeqs[num_seq] = NULL;
  94    names[num_seq] = NULL;
  95    if (num_seq == 0) {
  96      fprintf(stderr, "No sequences found in CLUSTAL/STOCKHOLM file\n");
  97      return 0;
  98    }
  99    n = strlen(AlignedSeqs[0]);
 100    for (nn=1; nn<num_seq; nn++) {
 101      if (strlen(AlignedSeqs[nn])!=n) {
 102        fprintf(stderr, "Sorry, your file is fucked up.\n"
 103                "Unequal lengths!\n\n");
 104        return 0;
 105      }
 106    }
 107
 108    fprintf(stderr, "%d sequences; length of alignment %d.\n", nn, n);
 109    return num_seq;
 110 }
 111
 112 char *consensus(const char *AS[]) {
 113   /* simple consensus sequence (most frequent character) */
 114   char *string;
 115   int i,n;
 116   n = strlen(AS[0]);
 117   string = (char *) space((n+1)*sizeof(char));
 118   for (i=0; i<n; i++) {
 119     int s,c,fm, freq[8] = {0,0,0,0,0,0,0,0};
 120     for (s=0; AS[s]!=NULL; s++)
 121       freq[encode_char(AS[s][i])]++;
 122     for (s=c=fm=0; s<8; s++) /* find the most frequent char */
 123       if (freq[s]>fm) {c=s, fm=freq[c];}
 124     if (s>4) s++; /* skip T */
 125     string[i]=Law_and_Order[c];
 126   }
 127   return string;
 128 }
 129
 130 /* IUP nucleotide classes indexed by a bit string of the present bases */
 131 /* A C AC G AG CG ACG U AU CU ACU GU AGU CGU ACGU */
 132 static char IUP[17] = "-ACMGRSVUWYHKDBN";
 133 char *consens_mis(const char*AS[]) {
 134   /* MIS displays the 'most informative sequence' (Freyhult et al 2004),
 135      elements in columns with frequency greater than the background
 136      frequency are projected into iupac notation. Columns where gaps are
 137      over-represented are in lower case. */
 138
 139   char *cons;
 140   int i, s, n, N, c;
 141   int bgfreq[8] = {0,0,0,0,0,0,0,0};
 142
 143   n = strlen(AS[0]);
 144   for (N=0; AS[N]!=NULL; N++);
 145   cons = (char *) space((n+1)*sizeof(char));
 146
 147   for (i=0; i<n; i++)
 148     for (s=0; s<N; s++) {
 149       c = encode_char(AS[s][i]);
 150       if (c>4) c=5;
 151       bgfreq[c]++;
 152     }
 153
 154   for (i=0; i<n; i++) {
 155     int freq[8] = {0,0,0,0,0,0,0,0};
 156     int code = 0;
 157     for (s=0; s<N; s++) {
 158       c = encode_char(AS[s][i]);
 159       if (c>4) c=5;
 160       freq[c]++;
 161     }
 162     for (c=4; c>0; c--) {
 163       code <<=1;
 164       if (freq[c]*n>=bgfreq[c]) code++;
 165     }
 166     cons[i] = IUP[code];
 167     if (freq[0]*n>bgfreq[0])
 168       cons[i] = tolower(IUP[code]);
 169   }
 170   return cons;
 171 }