/*
alphabet.c
Code for handling nucleotide and base pair alphabet
Part of the ViennaRNA Package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/alphabet.h"
/*
For now, we neglect all non-standard nucleotides in an input sequence, i.e. only
ACGTUN is allowed.
However, the standard nucleotide ambiguity code table would allow for many more:
A = Adenylic acid
C = Cytidylic acid
G = Guanylic acid
T = Thymidylic acid
U = Uridylic acid
I = Inosylic acid
R = A or G = puRine
Y = C or T = pYrimidine
K = G or T = Keto
M = A or C = aMino
S = G or C = Strong base pair
W = A or T = Weak base pair
B = not A (G or C or T)
D = not C (A or G or T)
H = not G (A or C or T)
V = not T/U (A or C or G)
N = aNy base (by convention, X is used for unknown amino acids, N for unknown nucleotides)
For the future, we aim to accept all of the above codes.
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
PRIVATE const char Law_and_Order[] = "_ACGUTXKI";
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE char *wrap_get_ptypes(const short *S, vrna_md_t *md); /* provides backward compatibility for old ptypes array in pf computations */
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC unsigned int
vrna_sequence_length_max(unsigned int options){
if(options & VRNA_OPTION_WINDOW)
return (unsigned int)INT_MAX;
/*
return (unsigned int)sqrt((double)INT_MAX);
*/
/*
many functions in RNAlib still rely on the sequence length
at pos 0 in the integer encoded sequence array S. Since this
encoding is stored in a short * array, the maximum length
of any sequence is SHRT_MAX
*/
return (unsigned int)SHRT_MAX;
}
PUBLIC int
vrna_nucleotide_IUPAC_identity( char nt,
char mask){
char n1,n2,*p;
p = NULL;
n1 = toupper(nt);
n2 = toupper(mask);
switch(n1){
case 'A': p = strchr("ARMWDHVN", n2);
break;
case 'C': p = strchr("CYMSBHVN", n2);
break;
case 'G': p = strchr("GRKSBDVN", n2);
break;
case 'T': p = strchr("TYKWBDHN", n2);
break;
case 'U': p = strchr("UYKWBDHN", n2);
break;
case 'I': p = strchr("IN", n2);
break;
case 'R': p = strchr("AGR", n2);
break;
case 'Y': p = strchr("CTUY", n2);
break;
case 'K': p = strchr("GTUK", n2);
break;
case 'M': p = strchr("ACM", n2);
break;
case 'S': p = strchr("GCS", n2);
break;
case 'W': p = strchr("ATUW", n2);
break;
case 'B': p = strchr("GCTBU", n2);
break;
case 'D': p = strchr("AGTUD", n2);
break;
case 'H': p = strchr("ACTUH", n2);
break;
case 'V': p = strchr("ACGV", n2);
break;
case 'N': p = strchr("ACGTUN", n2);
break;
}
return (p) ? 1 : 0;
}
PUBLIC char *
vrna_ptypes(const short *S,
vrna_md_t *md){
char *ptype;
int n,i,j,k,l,*idx;
int min_loop_size = md->min_loop_size;
n = S[0];
if((unsigned int)n > vrna_sequence_length_max(VRNA_OPTION_DEFAULT)){
vrna_message_warning("vrna_ptypes@alphabet.c: sequence length of %d exceeds addressable range", n);
return NULL;
}
ptype = (char *)vrna_alloc(sizeof(char)*((n*(n+1))/2+2));
idx = vrna_idx_col_wise(n);
for (k=1; k<n-min_loop_size; k++)
for (l=1; l<=2; l++) {
int type,ntype=0,otype=0;
i=k; j = i+min_loop_size+l; if (j>n) continue;
type = md->pair[S[i]][S[j]];
while ((i>=1)&&(j<=n)) {
if ((i>1)&&(j<n)) ntype = md->pair[S[i-1]][S[j+1]];
if (md->noLP && (!otype) && (!ntype))
type = 0; /* i.j can only form isolated pairs */
ptype[idx[j]+i] = (char) type;
otype = type;
type = ntype;
i--; j++;
}
}
free(idx);
return ptype;
}
PUBLIC short *
vrna_seq_encode(const char *sequence,
vrna_md_t *md){
unsigned int i, l;
short *S = NULL;
if(sequence && md){
S = vrna_seq_encode_simple(sequence, md);
l = (unsigned int)strlen(sequence);
for(i=1; i<=l; i++)
S[i] = md->alias[S[i]];
S[l+1] = S[1];
S[0] = S[l];
}
return S;
}
PUBLIC short *
vrna_seq_encode_simple( const char *sequence,
vrna_md_t *md){
unsigned int i, l;
short *S = NULL;
if(sequence && md){
l = (unsigned int)strlen(sequence);
S = (short *) vrna_alloc(sizeof(short)*(l+2));
for(i=1; i<=l; i++) /* make numerical encoding of sequence */
S[i]= (short) vrna_nucleotide_encode(toupper(sequence[i-1]), md);
S[l+1] = S[1];
S[0] = (short) l;
}
return S;
}
PUBLIC int
vrna_nucleotide_encode( char c,
vrna_md_t *md){
/* return numerical representation of nucleotide used e.g. in vrna_md_t.pair[][] */
int code = -1;
if(md){
if (md->energy_set>0) code = (int) (c-'A')+1;
else {
const char *pos;
pos = strchr(Law_and_Order, c);
if (pos==NULL) code=0;
else code = (int) (pos-Law_and_Order);
if (code>5) code = 0;
if (code>4) code--; /* make T and U equivalent */
}
}
return code;
}
PUBLIC char
vrna_nucleotide_decode( int enc,
vrna_md_t *md){
if(md){
if(md->energy_set > 0)
return (char)enc + 'A' - 1;
else
return (char)Law_and_Order[enc];
} else {
return (char)0;
}
}
PUBLIC void
vrna_aln_encode(const char *sequence,
short **S_p,
short **s5_p,
short **s3_p,
char **ss_p,
unsigned short **as_p,
vrna_md_t *md){
unsigned int i,l;
unsigned short p;
l = strlen(sequence);
(*s5_p) = (short *) vrna_alloc((l + 2) * sizeof(short));
(*s3_p) = (short *) vrna_alloc((l + 2) * sizeof(short));
(*as_p) = (unsigned short *)vrna_alloc((l + 2) * sizeof(unsigned short));
(*ss_p) = (char *) vrna_alloc((l + 2) * sizeof(char));
/* make numerical encoding of sequence */
(*S_p) = vrna_seq_encode_simple(sequence, md);
(*s5_p)[0] = (*s5_p)[1] = 0;
if(md->oldAliEn){
/* use alignment sequences in all energy evaluations */
(*ss_p)[0]=sequence[0];
for(i=1; i<l; i++){
(*s5_p)[i] = (*S_p)[i-1];
(*s3_p)[i] = (*S_p)[i+1];
(*ss_p)[i] = sequence[i];
(*as_p)[i] = i;
}
(*ss_p)[l] = sequence[l];
(*as_p)[l] = l;
(*s5_p)[l] = (*S_p)[l-1];
(*s3_p)[l] = 0;
(*S_p)[l+1] = (*S_p)[1];
(*s5_p)[1] = 0;
if(md->circ){
(*s5_p)[1] = (*S_p)[l];
(*s3_p)[l] = (*S_p)[1];
(*ss_p)[l+1] = (*S_p)[1];
}
}
else{
if(md->circ){
for(i=l; i>0; i--){
char c5;
c5 = sequence[i-1];
if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
(*s5_p)[1] = (*S_p)[i];
break;
}
for (i=1; i<=l; i++) {
char c3;
c3 = sequence[i-1];
if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
(*s3_p)[l] = (*S_p)[i];
break;
}
}
else (*s5_p)[1]=(*s3_p)[l]=0;
for(i=1,p=0; i<=l; i++){
char c5;
c5 = sequence[i-1];
if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
(*s5_p)[i+1]=(*s5_p)[i];
else { /* no gap */
(*ss_p)[p++]=sequence[i-1]; /*start at 0!!*/
(*s5_p)[i+1]=(*S_p)[i];
}
(*as_p)[i]=p;
}
for (i=l; i>=1; i--) {
char c3;
c3 = sequence[i-1];
if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
(*s3_p)[i-1]=(*s3_p)[i];
else
(*s3_p)[i-1]=(*S_p)[i];
}
}
}
PRIVATE char *
wrap_get_ptypes(const short *S,
vrna_md_t *md){
char *ptype;
int n,i,j,k,l,*idx;
n = S[0];
ptype = (char *)vrna_alloc(sizeof(char)*((n*(n+1))/2+2));
idx = vrna_idx_row_wise(n);
int min_loop_size = md->min_loop_size;
for (k=1; k<n-min_loop_size; k++)
for (l=1; l<=2; l++) {
int type,ntype=0,otype=0;
i=k; j = i+min_loop_size+l; if (j>n) continue;
type = md->pair[S[i]][S[j]];
while ((i>=1)&&(j<=n)) {
if ((i>1)&&(j<n)) ntype = md->pair[S[i-1]][S[j+1]];
if (md->noLP && (!otype) && (!ntype))
type = 0; /* i.j can only form isolated pairs */
ptype[idx[i]-j] = (char) type;
otype = type;
type = ntype;
i--; j++;
}
}
free(idx);
return ptype;
}
#ifdef VRNA_BACKWARD_COMPAT
/*###########################################*/
/*# deprecated functions below #*/
/*###########################################*/
PUBLIC char *
get_ptypes( const short *S,
vrna_md_t *md,
unsigned int idx_type){
if(S){
if((unsigned int)S[0] > vrna_sequence_length_max(VRNA_OPTION_DEFAULT)){
vrna_message_warning("get_ptypes@alphabet.c: sequence length of %d exceeds addressable range", (int)S[0]);
return NULL;
}
if(idx_type)
return wrap_get_ptypes(S, md);
else
return vrna_ptypes(S, md);
} else {
return NULL;
}
}
#endif