view release on metacpan or  search on metacpan

src/wurstsrc/score_fx.c  view on Meta::CPAN

    for (i = 0; i < s->size; i++) {
        for (j = i+4; j < s->size; j++) {
            float d = dist2 ( s->rp_cb[i], s->rp_cb[j]);
            if (d < n_cut) {             /* 8A cutoff CB-CB dist */
                ++cnt[i];
                ++cnt[j];
            }
        }
    }
                                        /* calc psi dihedral angle */
    for (i = 0; i < s->size-1; i++) {
        psi[i] = conv * dihedral(s->rp_n[i], s->rp_ca[i], s->rp_c[i],
                                 s->rp_n[i+1]);
    }

                                        /* run over all overlapping fragments
                                           and determine optimal class */
    for (i = 0; i < s->size - fx->nr_inst; i++) {
        float emax = -1.0e+10;
        int jmax = -1;

        float d;  /* CA-CA end-to-end distance */
        d = sqrt (dist2 (s->rp_ca[i], s->rp_ca[i+fx->nr_inst-1]));

        for (j = 0; j < fx->nr_groups; j++) {
            float tmp;
            ej = 0;
            for (k = 0; k < fx->nr_inst; k++) {
                float tmp2;
                                        /* Gaussian model for psi dihedral */
                dpsi = psi[i+k] - fx->psi[j][k];
                if (dpsi > 180)
                    dpsi -= 360;
                else if (dpsi < -180)
                    dpsi += 360;
                tmp2 = 2 * fx->dpsi[j][k] * fx->dpsi[j][k];
                ppsi = exp(-dpsi*dpsi/tmp2) / sqrt(M_PI*tmp2);

                if (ppsi > FX_TINY)
                    ej += log(ppsi);
                else
                    ej += log(FX_TINY);

                                /* discrete profile for number of neighbours */
                nbr = cnt[i+k];
                if (nbr >= 20)          /* just in case */
                    nbr = 19;
                ej += fx->pna[j][k][nbr];
            }

                                         /* because we use a non-linear scale
                                           we have to find the right bin */
            for ( k = 0; k < fx->nr_dbins; k++)
                if (fx->dbin[k] > d)
                    break;
            if (k)
                k--;

            dpsi = k - fx->pdav[j];
            tmp = 2 * fx->pdsig[j] * fx->pdsig[j];
            ppsi = exp(-dpsi*dpsi/tmp) / sqrt(M_PI*tmp);

            if (ppsi > FX_TINY)
                ej += log(ppsi);
            else
                ej += log(FX_TINY);


            if (ej > emax) {
                emax = ej;
                jmax = j;
            }
        }
        if (jmax != -1)
            class[i] = jmax;
        else {
            err_printf (this_sub, "BUGGER ALL: no class of fragment!!! \n");
            exit(EXIT_FAILURE);
        }
    }
    free(cnt);
    free(psi);
    return EXIT_SUCCESS;
}

/* ---------------- score_fx -------------------------------
 * Score a sequence and structure using a fragment based
 * scoring function.
 * The contents of the param should be whatever was returned
 * by the param reading routine.
 * Note, the matrix we get is has two extra rows and two extra
 * columns. One at the start and end.
 * You *are* allowed to assume that the matrix is zeroed.
 */

int
score_fx (struct score_mat *score_mat, struct seq *s,
          struct coord *c1, struct FXParam *fx)
{
    const char *this_sub = "score_fx";
    size_t   i;
    int      j, jlast, k, kk;
    int      aa1, classi;
    int      *class, middle, middle_plus1;
    float    ei, **scores = score_mat->mat;
    size_t   to_mall;

    /* It is almost impossible to be given obvious garbage by the
     * interpreter, but this is the kind of check and error exit
     * we do
     */
    if (score_mat == NULL || s == NULL || c1 == NULL || fx == NULL) {
        err_printf (this_sub, "null parameter, FIX \n");
        return (EXIT_FAILURE);
    }

    seq_std2thomas (s);          /* Force Thomas style names for amino acids */
    coord_a_2_nm (c1);     /* Force coordinates to nanometers if they were A */
    to_mall = c1->size * sizeof (class[0]);
    class = E_MALLOC (to_mall);              /* Bayes class of each fragment */