mdpbase.c File Reference

#include "../headers/mdpbase.h"

Go to the source code of this file.

Functions
void	store_vertice_positions (s_mdconcat *m, c_lst_pockets *pockets)
void	calculate_md_dens_grid (s_mdgrid *g, c_lst_pockets *pockets, s_mdparams *par)
void	update_md_grid (s_mdgrid *g, s_mdgrid *refg, c_lst_pockets *pockets, s_mdparams *par)
void	reset_grid (s_mdgrid *g)
void	project_grid_on_atoms (s_mdgrid *g, s_pdb *pdb)
void	normalize_grid (s_mdgrid *g, int n)
s_mdconcat *	init_md_concat (void)
s_min_max_pockets *	float_get_min_max_from_pockets (c_lst_pockets *pockets)
s_mdgrid *	init_md_grid (c_lst_pockets *pockets)
void	alloc_first_md_concat (s_mdconcat *m, size_t n)
void	realloc_md_concat (s_mdconcat *m, size_t n)
void	free_mdconcat (s_mdconcat *m)

---

Function Documentation

void alloc_first_md_concat	(	s_mdconcat *	m,
		size_t	n
	)

## FUNCTION: alloc_first_md_concat

## SPECIFICATION: Allocate memory for the md concat structure (first snapshot)

## PARAMETRES: @ s_mdconcat *m: Pointer to the structure to free, @ size_t n: Number of vertices in the snapshot to add to m

## RETURN: void

Definition at line 490 of file mdpbase.c.

Referenced by store_vertice_positions().

                                                    {
    size_t z;
    m->vertpos = (float **) my_malloc(sizeof (float *) * n);
    for (z = 0; z < n; z++) {
        m->vertpos[z] = (float *) my_malloc(sizeof (float) *4);
        m->vertpos[z][0] = 0.0;
        m->vertpos[z][1] = 0.0;
        m->vertpos[z][2] = 0.0;
        m->vertpos[z][3] = 0.0;
    }
    m->n_vertpos = n;
}

void calculate_md_dens_grid	(	s_mdgrid *	g,
		c_lst_pockets *	pockets,
		s_mdparams *	par
	)

## FUNCTION: calculate_md_dens_grid

## SPECIFICATION: Do the actual grid calculation (count number of vertices near a grid point) NOT USED ANYMORE ## PARAMETRES: @ s_mdgrid *g: Grid structure @ c_lst_pockets *pockets : Pocket chained list @ s_mdparams *par : Parameters of mdpocket

## RETURN: void

Definition at line 104 of file mdpbase.c.

References s_desc::drug_score, c_lst_vertices::first, c_lst_pockets::first, s_mdparams::flag_scoring, s_mdgrid::gridvalues, M_MDP_CUBE_SIDE, node_pocket::next, node_vertice::next, s_mdgrid::origin, s_pocket::pdesc, node_pocket::pocket, s_mdgrid::resolution, s_pocket::v_lst, node_vertice::vertice, s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by mdpocket_detect().

                                                                                  {
    int xidx = 0, yidx = 0, zidx = 0; /*direct indices of the positions in the grid*/
    int sxidx, syidx, szidx; /*indices nearby xidx,yidx,zidx that are within M_MDP_CUBE_SIDE of this point*/
    int sxi = 0, syi = 0, szi = 0;
    float vx, vy, vz;
    float incre = 1.0;
    node_pocket *cp = NULL;
    s_pocket *pocket = NULL;
    node_vertice *cnv = NULL;
    s_vvertice *cv = NULL;
    int s = (int) (((float) M_MDP_CUBE_SIDE / 2.0) / g->resolution); /*possible stepsize (resolution dependent) in the discrete grid*/
    /*loop over all known vertices and CALCULATE the grid positions and increment grid values by 1*/
    /*important : no distance calculations are done here, thus this routine is very fast*/
    cp = pockets->first;
    while (cp) { /*loop over pockets*/
        pocket = cp->pocket;
        cnv = pocket->v_lst->first;
        /*if(par->flag_scoring) incre=pocket->score;*/
        if (par->flag_scoring) incre = pocket->pdesc->drug_score; /*if we score by drug score, use it as increment here instead of the 1 by default*/
        while (cnv) { /*loop over all vertices of a pocket*/
            cv = cnv->vertice;
            vx = cv->x; /*tmp the vertice position*/
            vy = cv->y;
            vz = cv->z;
            xidx = (int) roundf((vx - g->origin[0]) / g->resolution); /*calculate the nearest grid point internal coordinates*/
            yidx = (int) roundf((vy - g->origin[1]) / g->resolution);
            zidx = (int) roundf((vz - g->origin[2]) / g->resolution);
            /*we use a cube of M_MDP_CUBE_SIDE**3 volume to check if there are grid points nearby,
             *thus a check in the neighbourhood is performed here...this looks a bit awful and should be changed later on*/
            for (sxi = -s; sxi <= s; sxi++) {
                for (syi = -s; syi <= s; syi++) {
                    for (szi = -s; szi <= s; szi++) {
                        /*check if we are not on the point xidx,yidx, zidx....the square sum is just a little trick*/
                        if ((sxi * sxi + syi * syi + szi * szi) != 0) {
                            /*next we have to check in a discrete manner if we can include a new grid point or not.*/
                            if (sxi < 0) sxidx = (int) ceil(((float) sxi + vx - g->origin[0]) / g->resolution);
                            else if (sxi > 0) sxidx = (int) floor(((float) sxi + vx - g->origin[0]) / g->resolution);
                            else if (sxi == 0) sxidx = xidx;
                            if (syi < 0) syidx = (int) ceil(((float) syi + vy - g->origin[1]) / g->resolution);
                            else if (syi > 0) syidx = (int) floor(((float) syi + vy - g->origin[1]) / g->resolution);
                            else if (syi == 0) syidx = yidx;
                            if (szi < 0) szidx = (int) ceil(((float) szi + vz - g->origin[2]) / g->resolution);
                            else if (szi > 0) szidx = (int) floor(((float) szi + vz - g->origin[2]) / g->resolution);
                            else if (szi == 0) szidx = zidx;
                            /*double check if we are not on the already incremented grid point*/
                            if (((sxidx != xidx) || (syidx != yidx) || (szidx != zidx))
                                    && (sxidx >= 0 && syidx >= 0 && szidx >= 0 && sxidx < g->nx && syidx < g->ny && szidx < g->nz)) {
                                g->gridvalues[sxidx][syidx][szidx] += incre;
                            } /*added condition if grid value already set, do not update, to better density measurements and match drug score*/
                        }
                    }
                }
            }
            cnv = cnv->next;
        }

        if (xidx < g->nx && yidx < g->ny && zidx < g->nz && xidx >= 0 && yidx >= 0 && zidx >= 0) {

            g->gridvalues[xidx][yidx][zidx] += incre; /*increment the already known grid point value*/

        } else fprintf(stderr, "\n\nWarning (oh oh!!) : Your structure is not aligned or is moving a lot. Results might not reflect what you expect. Consider first structural alignemnt\nIf your structure is moving a lot, consider to split up analysis in two distinct parts.\n\n");
        cp = cp->next;
    }
}

s_min_max_pockets* float_get_min_max_from_pockets

(

c_lst_pockets *

pockets

)

## FUNCTION: float_get_min_max_from_pockets

## SPECIFICATION: Get the absolute minimum and maximum point from pockets

## PARAMETRES: @ c_lst_pockets *pockets : Chained list of pockets

## RETURN: @ s_min_max_pockets : Structure containing the minimum & maximum

Definition at line 397 of file mdpbase.c.

References s_min_max_pockets::maxx, s_min_max_pockets::maxy, s_min_max_pockets::maxz, s_min_max_pockets::minx, s_min_max_pockets::miny, s_min_max_pockets::minz, my_malloc(), s_lst_vvertice::nvert, s_lst_vvertice::vertices, c_lst_pockets::vertices, s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by init_md_grid().

                                                                          {
    if (pockets) {
        int z;
        float minx = 50000., maxx = -50000., miny = 50000., maxy = -50000., minz = 50000., maxz = -50000.;
        int n = pockets->vertices->nvert;
        /*if there a no vertices in m before, first allocate some space*/
        for (z = 0; z < n; z++) { /*loop over all vertices*/
            /*store the positions and radius of the vertices*/
            if (minx > pockets->vertices->vertices[z].x) minx = pockets->vertices->vertices[z].x;
            else if (maxx < pockets->vertices->vertices[z].x)maxx = pockets->vertices->vertices[z].x;
            if (miny > pockets->vertices->vertices[z].y) miny = pockets->vertices->vertices[z].y;
            else if (maxy < pockets->vertices->vertices[z].y)maxy = pockets->vertices->vertices[z].y;
            if (minz > pockets->vertices->vertices[z].z) minz = pockets->vertices->vertices[z].z;
            else if (maxz < pockets->vertices->vertices[z].z)maxz = pockets->vertices->vertices[z].z;
        }
        s_min_max_pockets *r = (s_min_max_pockets *) my_malloc(sizeof (s_min_max_pockets));
        r->maxx = maxx;
        r->maxy = maxy;
        r->maxz = maxz;
        r->minx = minx;
        r->miny = miny;
        r->minz = minz;
        return (r);
    }
    return (NULL);
}

void free_mdconcat

(

s_mdconcat *

m

)

## FUNCTION: free_mdconcat

## SPECIFICATION: Free the mdconcat structure

## PARAMETRES: @ s_mdconcat *m: Pointer to the structure to free

## RETURN: void

Definition at line 545 of file mdpbase.c.

                                  {
    size_t i;
    for (i = 0; i < m->n_vertpos; i++) {
        my_free(m->vertpos[i]);
    }
    my_free(m->vertpos);
    my_free(m);
}

s_mdconcat* init_md_concat

(

void

)

## FUNCTION: init_md_concat

## SPECIFICATION: Initialize the md concat (alloc)

## PARAMETRES: void

## RETURN: s_mdconcat * : the md concat structure

Definition at line 374 of file mdpbase.c.

                                 {
    s_mdconcat *m = (s_mdconcat *) my_malloc(sizeof (s_mdconcat));
    m->n_vertpos = 0;
    m->vertpos = NULL;
    m->n_snapshots = 0;
    return m;
}

s_mdgrid* init_md_grid

(

c_lst_pockets *

pockets

)

## FUNCTION: init_md_grid

## SPECIFICATION: Initialize the md grid (allocate + 0)

## PARAMETRES: @ s_mdconcat *mdc: Pointer to the mdconcat structure (vertices),

## RETURN: s_mdgrid * : the grid

Definition at line 438 of file mdpbase.c.

References float_get_min_max_from_pockets(), s_mdgrid::gridvalues, M_MDP_CUBE_SIDE, M_MDP_GRID_RESOLUTION, s_min_max_pockets::maxx, s_min_max_pockets::maxy, s_min_max_pockets::maxz, s_min_max_pockets::minx, s_min_max_pockets::miny, s_min_max_pockets::minz, my_free(), my_malloc(), s_mdgrid::nx, s_mdgrid::ny, s_mdgrid::nz, s_mdgrid::origin, and s_mdgrid::resolution.

Referenced by mdpocket_detect().

                                               {
    s_mdgrid *g = (s_mdgrid *) my_malloc(sizeof (s_mdgrid));
    s_min_max_pockets *mm = float_get_min_max_from_pockets(pockets);
    float xmax = mm->maxx;
    float ymax = mm->maxy;
    float zmax = mm->maxz;
    float xmin = mm->minx;
    float ymin = mm->miny;
    float zmin = mm->minz;
    float initvalue = 0.0;
    my_free(mm);
    int cx, cy, cz;
    float span = (M_MDP_CUBE_SIDE / 2.0) / M_MDP_GRID_RESOLUTION;
    g->resolution = M_MDP_GRID_RESOLUTION;

    g->nx = 1 + (int) (xmax + 30. * span - xmin) / (g->resolution);
    g->ny = 1 + (int) (ymax + 30. * span - ymin) / (g->resolution);
    g->nz = 1 + (int) (zmax + 30. * span - zmin) / (g->resolution);

    g->gridvalues = (float ***) my_malloc(sizeof (float **) * g->nx);
    for (cx = 0; cx < g->nx; cx++) {
        g->gridvalues[cx] = (float **) my_malloc(sizeof (float *) * g->ny);
        for (cy = 0; cy < g->ny; cy++) {
            g->gridvalues[cx][cy] = (float *) my_malloc(sizeof (float) * g->nz);
            for (cz = 0; cz < g->nz; cz++) {
                g->gridvalues[cx][cy][cz] = initvalue;

            }
        }
    }

    g->origin = (float *) my_malloc(sizeof (float) *3);
    g->origin[0] = xmin - 15. * span;
    g->origin[1] = ymin - 15. * span;
    g->origin[2] = zmin - 15. * span;
    return g;
}

void normalize_grid	(	s_mdgrid *	g,
		int	n
	)

## FUNCTION: normalize_grid

## SPECIFICATION: Normalize all grid values by the number of snapshots

## PARAMETRES: @ s_mdgrid *g: Grid structure @ int n : Number of snapshots ## RETURN: void

Definition at line 348 of file mdpbase.c.

References s_mdgrid::gridvalues, s_mdgrid::nx, s_mdgrid::ny, and s_mdgrid::nz.

Referenced by mdpocket_detect().

                                        {
    int x, y, z;
    for (x = 0; x < g->nx; x++) {
        for (y = 0; y < g->ny; y++) {
            for (z = 0; z < g->nz; z++) {
                g->gridvalues[x][y][z] /= (float) n;
            }
        }
    }
}

void project_grid_on_atoms	(	s_mdgrid *	g,
		s_pdb *	pdb
	)

## FUNCTION: project_grid_on_atoms

## SPECIFICATION: Project grid values to atoms to store these into the bfactor column of the pdb

## PARAMETRES: @ s_mdgrid *g: Grid structure ## RETURN: void

Definition at line 295 of file mdpbase.c.

References s_atm::bfactor, s_mdgrid::gridvalues, s_pdb::latoms_p, M_MDP_ATOM_DENSITY_DIST, s_mdgrid::n_snapshots, s_pdb::natoms, s_mdgrid::origin, s_mdgrid::resolution, s_atm::x, s_atm::y, and s_atm::z.

Referenced by mdpocket_detect().

                                                    {
    s_atm *ca = NULL;
    float x, y, z;
    int ix, iy, iz;
    float maxix, maxiy, maxiz;
    int i;
    float density;
    float ngridpoints1D = (float) M_MDP_ATOM_DENSITY_DIST / g->resolution;
    for (i = 0; i < pdb->natoms; i++) {/*for all protein atoms*/
        density = 0.0;
        ca = pdb->latoms_p[i];  /*set up the current atom*/
        /*get the coordinates*/
        x = ca->x;
        y = ca->y;
        z = ca->z;
        /*nearest position of the atom in x in the grid*/
        ix = (int) floor((x - M_MDP_ATOM_DENSITY_DIST - g->origin[0]) / g->resolution);

        maxix = (int) ceil((x + M_MDP_ATOM_DENSITY_DIST - g->origin[0]) / g->resolution);
        maxiy = (int) ceil((y + M_MDP_ATOM_DENSITY_DIST - g->origin[1]) / g->resolution);
        maxiz = (int) ceil((z + M_MDP_ATOM_DENSITY_DIST - g->origin[2]) / g->resolution);
        while (ix <= maxix && g->nx > ix) {
            iy = (int) floor((y - M_MDP_ATOM_DENSITY_DIST - g->origin[1]) / g->resolution);
            while (iy <= maxiy && g->ny > iy) {
                iz = (int) floor((z - M_MDP_ATOM_DENSITY_DIST - g->origin[2]) / g->resolution);
                while (iz <= maxiz && g->nz > iz) {
                    density += g->gridvalues[ix][iy][iz];
                    iz++;
                }
                iy++;
            }
            ix++;
        }
        /*transform the density to something more contrasted for the bfactor column*/
        ca->bfactor = log(1 + density / (ngridpoints1D) / (float) g->n_snapshots);
    }
}

void realloc_md_concat	(	s_mdconcat *	m,
		size_t	n
	)

## FUNCTION: realloc_md_concat

## SPECIFICATION: Reallocate memory for the md concat structure (to add a new snapshot)

## PARAMETRES: @ s_mdconcat *m: Pointer to the structure to free, @ size_t n: Number of vertices in the snapshot to add to m

## RETURN: void

Definition at line 518 of file mdpbase.c.

Referenced by store_vertice_positions().

                                                {
    size_t z;
    m->vertpos = (float **) my_realloc(m->vertpos, sizeof (float *) *(m->n_vertpos + n));
    for (z = 0; z < n; z++) {
        m->vertpos[m->n_vertpos + z] = (float *) my_malloc(sizeof (float) *4);
        m->vertpos[m->n_vertpos + z][0] = 0.0;
        m->vertpos[m->n_vertpos + z][1] = 0.0;
        m->vertpos[m->n_vertpos + z][2] = 0.0;
        m->vertpos[m->n_vertpos + z][3] = 0.0;
    }
    m->n_vertpos += n;
}

void reset_grid

(

s_mdgrid *

g

)

## FUNCTION: reset_grid

## SPECIFICATION: Reset all grid values to 0

## PARAMETRES: @ s_mdgrid *g: Grid structure ## RETURN: void

Definition at line 269 of file mdpbase.c.

References s_mdgrid::gridvalues, s_mdgrid::nx, s_mdgrid::ny, and s_mdgrid::nz.

Referenced by mdpocket_detect().

                             {
    int x, y, z;
    for (x = 0; x < g->nx; x++) {
        for (y = 0; y < g->ny; y++) {
            for (z = 0; z < g->nz; z++) {
                g->gridvalues[x][y][z] = 0.0;
            }
        }
    }
}

void store_vertice_positions	(	s_mdconcat *	m,
		c_lst_pockets *	pockets
	)

## FUNCTION: store_vertice_positions

## SPECIFICATION: Store vertice positions of current pockets in the mdconcat structure

## PARAMETRES: @ s_mdconcat *m: Pointer to the mdconcat structure (vertices), @ c_lst_pockets *pockets : Pointer to all pockets found by fpocket

## RETURN: void

Definition at line 70 of file mdpbase.c.

References alloc_first_md_concat(), s_mdconcat::n_vertpos, s_lst_vvertice::nvert, s_vvertice::ray, realloc_md_concat(), s_lst_vvertice::vertices, c_lst_pockets::vertices, s_mdconcat::vertpos, s_vvertice::x, s_vvertice::y, and s_vvertice::z.

                                                                    {
    if (pockets) {
        int z;
        int n = pockets->vertices->nvert;
        size_t old_n = m->n_vertpos;
        /*if there a no vertices in m before, first allocate some space*/
        if (m->n_vertpos == 0) alloc_first_md_concat(m, n);
        else realloc_md_concat(m, n); /*else, reallocate the object size*/
        for (z = 0; z < n; z++) { /*loop over all vertices*/
            /*store the positions and radius of the vertices*/
            m->vertpos[old_n + z][0] = pockets->vertices->vertices[z].x;
            m->vertpos[old_n + z][1] = pockets->vertices->vertices[z].y;
            m->vertpos[old_n + z][2] = pockets->vertices->vertices[z].z;
            m->vertpos[old_n + z][3] = pockets->vertices->vertices[z].ray;
        }
    }
}

void update_md_grid	(	s_mdgrid *	g,
		s_mdgrid *	refg,
		c_lst_pockets *	pockets,
		s_mdparams *	par
	)

## FUNCTION: update_md_dens_grid

## SPECIFICATION: Do the actual grid calculation (count number of vertices near a grid point)

## PARAMETRES: @ s_mdgrid *g: Grid structure @ c_lst_pockets *pockets : Pocket chained list @ s_mdparams *par : Parameters of mdpocket ## RETURN: void

Definition at line 184 of file mdpbase.c.

References s_desc::drug_score, c_lst_vertices::first, c_lst_pockets::first, s_mdparams::flag_scoring, s_mdgrid::gridvalues, M_MDP_CUBE_SIDE, node_pocket::next, node_vertice::next, s_mdgrid::origin, s_pocket::pdesc, node_pocket::pocket, s_mdgrid::resolution, s_pocket::v_lst, node_vertice::vertice, s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by mdpocket_detect().

                                                                                          {
    int xidx = 0, yidx = 0, zidx = 0; /*direct indices of the positions in the grid*/
    int sxidx = 0, syidx = 0, szidx = 0; /*indices nearby xidx,yidx,zidx that are within M_MDP_CUBE_SIDE of this point*/
    int sxi = 0, syi = 0, szi = 0;
    float vx, vy, vz;
    float incre = 1.0;
    node_pocket *cp = NULL;
    s_pocket *pocket = NULL;
    node_vertice *cnv = NULL;
    s_vvertice *cv = NULL;
    int s = (int) (((float) M_MDP_CUBE_SIDE / 2.0) / g->resolution); /*possible stepsize (resolution dependent) in the discrete grid*/
    /*loop over all known vertices and CALCULATE the grid positions and increment grid values by 1*/
    /*important : no distance calculations are done here, thus this routine is very fast*/
    cp = pockets->first;
    while (cp) { /*loop over pockets*/
        pocket = cp->pocket;
        cnv = pocket->v_lst->first;
        /*if(par->flag_scoring) incre=pocket->score;*/
        if (par->flag_scoring) incre = pocket->pdesc->drug_score;
        while (cnv) { /*loop over vertices*/
            cv = cnv->vertice;
            vx = cv->x; /*tmp the vertice position*/
            vy = cv->y;
            vz = cv->z;
            xidx = (int) roundf((vx - g->origin[0]) / g->resolution); /*calculate the nearest grid point internal coordinates*/
            yidx = (int) roundf((vy - g->origin[1]) / g->resolution);
            zidx = (int) roundf((vz - g->origin[2]) / g->resolution);
            /*we use a cube of M_MDP_CUBE_SIDE**3 volume to check if there are grid points nearby,
             *thus a check in the neighbourhood is performed here...this looks a bit awful and should be changed later on*/

            for (sxi = -s; sxi <= s; sxi++) {
                for (syi = -s; syi <= s; syi++) {
                    for (szi = -s; szi <= s; szi++) {
                        /*check if we are not on the point xidx,yidx, zidx....the square sum is just a little trick*/
                        if ((sxi * sxi + syi * syi + szi * szi) != 0) {
                            /*next we have to check in a discrete manner if we can include a new grid point or not.*/
                            if (sxi < 0) sxidx = (int) ceil(((float) sxi + vx - g->origin[0]) / g->resolution);
                            else if (sxi > 0) sxidx = (int) floor(((float) sxi + vx - g->origin[0]) / g->resolution);
                            else if (sxi == 0) sxidx = xidx;
                            if (syi < 0) syidx = (int) ceil(((float) syi + vy - g->origin[1]) / g->resolution);
                            else if (syi > 0) syidx = (int) floor(((float) syi + vy - g->origin[1]) / g->resolution);
                            else if (syi == 0) syidx = yidx;
                            if (szi < 0) szidx = (int) ceil(((float) szi + vz - g->origin[2]) / g->resolution);
                            else if (szi > 0) szidx = (int) floor(((float) szi + vz - g->origin[2]) / g->resolution);
                            else if (szi == 0) szidx = zidx;
                            /*double check if we are not on the already incremented grid point*/
                            if (((sxidx != xidx) || (syidx != yidx) || (szidx != zidx))
                                    && (sxidx >= 0 && syidx >= 0 && szidx >= 0 && sxidx < g->nx && syidx < g->ny && szidx < g->nz)) {
                                if (refg->gridvalues[sxidx][syidx][szidx] < 1.0) {
                                    g->gridvalues[sxidx][syidx][szidx] += incre;
                                    refg->gridvalues[sxidx][syidx][szidx] = 1.0;

                                }
                            } /*added condition if grid value already set, do not update, to better density measurements and match drug score*/
                        }
                    }
                }
            }
            cnv = cnv->next;
        }

        if (xidx < g->nx && yidx < g->ny && zidx < g->nz && xidx >= 0 && yidx >= 0 && zidx >= 0) {
            if (refg->gridvalues[xidx][yidx][zidx] < 1.0) {
                g->gridvalues[xidx][yidx][zidx] += incre; /*increment the already known grid point value*/
                refg->gridvalues[xidx][yidx][zidx] = 1.0;
            }
        } else fprintf(stderr, "\n\nWarning (oh oh!!) : Your structure is not aligned or is moving a lot. Results might not reflect what you expect. Consider first structural alignemnt\nIf your structure is moving a lot, consider to split up analysis in two distinct parts.\n\n");
        cp = cp->next;
    }
    //return g;
}