ndconv.py

# coding=utf-8
"""
                            Network deconvolution tool
"""
from __future__ import absolute_import, division, print_function

import logging
import numpy as np


LOG = logging.getLogger(__name__)


def net_deconv(npmat, beta=0.99, alpha=1, control=0):
    """
    This is a python implementation/translation of network deconvolution by
    MIT-KELLIS LAB
    
    LICENSE: MIT-KELLIS LAB
    
    AUTHORS:
        Algorithm was programmed by Soheil Feizi.
        Paper authors are S. Feizi, D. Marbach,  M. Medard and M. Kellis

    Python implementation: Gideon Rosenthal
    
    REFERENCES:
        For more details, see the following paper:
        Network Deconvolution as a General Method to Distinguish
        Direct Dependencies over Networks
        By: Soheil Feizi, Daniel Marbach,  Muriel Medard and Manolis Kellis
        Nature Biotechnology

    DESCRIPTION:
    
    USAGE:
        mat_nd = ND(npmat)
        mat_nd = ND(npmat,beta)
        mat_nd = ND(npmat,beta,alpha,control)
    
    
    INPUT ARGUMENTS:

    Parameters
    ----------
    npmat :
        Input matrix, if it is a square matrix, the program assumes
        it is a relevance matrix where npmat(i,j) represents the
        similarity content between nodes i and j. Elements of
        matrix should be non-negative.
        optional parameters:
    beta :
        Scaling parameter, the program maps the largest absolute
        eigenvalue of the direct dependency matrix to beta. It
        should be between 0 and 1. (Default value = 0.99)
    alpha :
        fraction of edges of the observed dependency matrix to be
        kept in deconvolution process. (Default value = 1)
    control :
        if 0, displaying direct weights for observed
        interactions, if 1, displaying direct weights for both
        observed and
        non-observed interactions.
        OUTPUT ARGUMENTS:
        mat_nd        Output deconvolved matrix (direct dependency matrix). Its
        components
        represent direct edge weights of observed interactions.
        Choosing top direct interactions (a cut-off) depends on
        the application and
        is not implemented in this code.
        To apply ND on regulatory networks, follow steps explained in
        Supplementary notes
        1.4.1 and 2.1 and 2.3 of the paper.
        In this implementation, input matrices are made symmetric.
        **************************************************************************
        loading scaling and thresholding parameters (Default value = 0)

    Returns
    -------

    
    """
    import scipy.stats.mstats as stat
    from numpy import linalg as la

    if beta >= 1 or beta <= 0:
        LOG.error('error: beta should be in (0,1)')

    if alpha > 1 or alpha <= 0:
        LOG.error('error: alpha should be in (0,1)')

    '''
    ***********************************
     Processing the input matrix
     diagonal values are filtered
    '''

    # n = npmat.shape[0]
    np.fill_diagonal(npmat, 0)

    '''
    Thresholding the input matrix
    '''
    y = stat.mquantiles(npmat[:], prob=[1 - alpha])
    th = npmat >= y
    mat_th = npmat * th

    '''
    making the matrix symetric if already not
    '''
    mat_th = (mat_th + mat_th.T) / 2

    '''
    ***********************************
    eigen decomposition
    '''
    LOG.info('Decomposition and deconvolution...')

    # noinspection PyTypeChecker
    dv, u = la.eigh(mat_th)
    d = np.diag(dv)
    lam_n = np.abs(np.min(np.min(np.diag(d)), 0))
    lam_p = np.abs(np.max(np.max(np.diag(d)), 0))

    m1 = lam_p * (1 - beta) / beta
    m2 = lam_n * (1 + beta) / beta
    m = max(m1, m2)

    # network deconvolution
    for i in range(d.shape[0]):
        d[i, i] = (d[i, i]) / (m + d[i, i])

    mat_new1 = np.dot(u, np.dot(d, la.inv(u)))

    '''

    ***********************************
     displaying direct weights
    '''
    if control == 0:
        ind_edges = (mat_th > 0) * 1.0
        ind_nonedges = (mat_th == 0) * 1.0
        m1 = np.max(np.max(npmat * ind_nonedges))
        m2 = np.min(np.min(mat_new1))
        mat_new2 = (mat_new1 + np.max(m1 - m2, 0)) * ind_edges + (
            npmat * ind_nonedges)
    else:
        m2 = np.min(np.min(mat_new1))
        mat_new2 = (mat_new1 + np.max(-m2, 0))

    '''
    ***********************************
     linearly mapping the deconvolved matrix to be between 0 and 1
    '''
    m1 = np.min(np.min(mat_new2))
    m2 = np.max(np.max(mat_new2))
    mat_nd = (mat_new2 - m1) / (m2 - m1)

    return mat_nd