/*
centroid structure prediction
Ivo L Hofacker + Ronny Lorenz
Vienna RNA package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/gquad.h"
#include "ViennaRNA/centroid.h"
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC char *
vrna_centroid_from_plist( int length,
double *dist,
vrna_plist_t *pl){
/* compute the centroid structure of the ensemble, i.e. the strutcure
with the minimal average distance to all other structures
<d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
Thus, the centroid is simply the structure containing all pairs with
p_ij>0.5 */
int i;
char *centroid;
if (pl==NULL)
vrna_message_error("vrna_centroid_from_plist@centroid.c: pl==NULL!");
*dist = 0.;
centroid = (char *) vrna_alloc((length+1)*sizeof(char));
for (i=0; i<length; i++) centroid[i]='.';
for (i=0; pl[i].i>0; i++){
if ((pl[i].p)>0.5) {
centroid[pl[i].i-1] = '(';
centroid[pl[i].j-1] = ')';
*dist += (1-pl[i].p);
} else
*dist += pl[i].p;
}
centroid[length] = '\0';
return centroid;
}
PUBLIC char *
vrna_centroid_from_probs( int length,
double *dist,
FLT_OR_DBL *probs){
/* compute the centroid structure of the ensemble, i.e. the strutcure
with the minimal average distance to all other structures
<d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
Thus, the centroid is simply the structure containing all pairs with
p_ij>0.5 */
int i,j;
FLT_OR_DBL p;
char *centroid;
int *index = vrna_idx_row_wise(length);
if (probs == NULL)
vrna_message_error("vrna_centroid_from_probs@centroid.c: probs==NULL!");
*dist = 0.;
centroid = (char *) vrna_alloc((length+1)*sizeof(char));
for (i=0; i<length; i++) centroid[i]='.';
for (i=1; i<=length; i++)
for (j=i+TURN+1; j<=length; j++) {
if ((p=probs[index[i]-j])>0.5) {
centroid[i-1] = '(';
centroid[j-1] = ')';
*dist += (1-p);
} else
*dist += p;
}
free(index);
centroid[length] = '\0';
return centroid;
}
PUBLIC char *
vrna_centroid(vrna_fold_compound_t *vc,
double *dist){
/* compute the centroid structure of the ensemble, i.e. the strutcure
with the minimal average distance to all other structures
<d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
Thus, the centroid is simply the structure containing all pairs with
p_ij>0.5 */
int i,j, k, length;
FLT_OR_DBL p;
char *centroid;
short *S;
vrna_mx_pf_t *matrices;
FLT_OR_DBL *probs;
int *my_iindx;
vrna_exp_param_t *pf_params;
if(!vc){
vrna_message_error("vrna_centroid@centroid.c: run vrna_pf_fold first!");
} else if( !vc->exp_matrices->probs){
vrna_message_error("vrna_centroid@centroid.c: probs==NULL!");
}
length = vc->length;
pf_params = vc->exp_params;
S = vc->sequence_encoding2;
my_iindx = vc->iindx;
matrices = vc->exp_matrices;
probs = matrices->probs;
*dist = 0.;
centroid = (char *) vrna_alloc((length+1)*sizeof(char));
for (i=0; i<length; i++) centroid[i]='.';
for (i=1; i<=length; i++)
for (j=i+TURN+1; j<=length; j++) {
if ((p=probs[my_iindx[i]-j])>0.5) {
if(pf_params->model_details.gquad){
/* check for presence of gquadruplex */
if((S[i] == 3) && (S[j] == 3)){
int L, l[3];
get_gquad_pattern_pf(S, i, j, pf_params, &L, l);
for(k=0;k<L;k++){
centroid[i+k-1]\
= centroid[i+k+L+l[0]-1]\
= centroid[i+k+2*L+l[0]+l[1]-1]\
= centroid[i+k+3*L+l[0]+l[1]+l[2]-1]\
= '+';
}
/* skip everything within the gquad */
i = j; j = j+TURN+1;
*dist += (1-p); /* right? */
break;
}
}
/* regular base pair */
centroid[i-1] = '(';
centroid[j-1] = ')';
*dist += (1-p);
} else
*dist += p;
}
centroid[length] = '\0';
return centroid;
}
/*###########################################*/
/*# deprecated functions below #*/
/*###########################################*/
/* this function is a threadsafe replacement for centroid() */
PUBLIC char *
get_centroid_struct_pl( int length,
double *dist,
vrna_plist_t *pl){
return vrna_centroid_from_plist(length, dist, pl);
}
/* this function is a threadsafe replacement for centroid() */
PUBLIC char *
get_centroid_struct_pr( int length,
double *dist,
FLT_OR_DBL *probs){
return vrna_centroid_from_probs(length, dist, probs);
}