ndconv.py 4.45 KB
Newer Older
Fabrice Allain's avatar
Fabrice Allain committed
1
# coding=utf-8
2
3
4
"""
                            Network deconvolution tool
"""
fabrice's avatar
fabrice committed
5
6
from __future__ import absolute_import, division, print_function

Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
7
import logging
fabrice's avatar
fabrice committed
8
9
import numpy as np

fabrice's avatar
fabrice committed
10

Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
11
12
13
LOG = logging.getLogger(__name__)


14
def net_deconv(npmat, beta=0.99, alpha=1, control=0):
15
16
    """
    This is a python implementation/translation of network deconvolution by
fabrice's avatar
fabrice committed
17
    MIT-KELLIS LAB
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
18
    
19
    LICENSE: MIT-KELLIS LAB
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
20
    
21
    AUTHORS:
fabrice's avatar
fabrice committed
22
23
        Algorithm was programmed by Soheil Feizi.
        Paper authors are S. Feizi, D. Marbach,  M. Medard and M. Kellis
24

fabrice's avatar
fabrice committed
25
    Python implementation: Gideon Rosenthal
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
26
    
fabrice's avatar
fabrice committed
27
    REFERENCES:
28
        For more details, see the following paper:
fabrice's avatar
fabrice committed
29
30
31
32
        Network Deconvolution as a General Method to Distinguish
        Direct Dependencies over Networks
        By: Soheil Feizi, Daniel Marbach,  Muriel Medard and Manolis Kellis
        Nature Biotechnology
33

fabrice's avatar
fabrice committed
34
    DESCRIPTION:
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
35
    
36
37
38
39
    USAGE:
        mat_nd = ND(npmat)
        mat_nd = ND(npmat,beta)
        mat_nd = ND(npmat,beta,alpha,control)
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
40
41
    
    
fabrice's avatar
fabrice committed
42
    INPUT ARGUMENTS:
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

    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
    -------

    
fabrice's avatar
fabrice committed
82
83
84
85
86
    """
    import scipy.stats.mstats as stat
    from numpy import linalg as la

    if beta >= 1 or beta <= 0:
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
87
        LOG.error('error: beta should be in (0,1)')
fabrice's avatar
fabrice committed
88
89

    if alpha > 1 or alpha <= 0:
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
90
        LOG.error('error: alpha should be in (0,1)')
fabrice's avatar
fabrice committed
91
92
93
94
95
96
97

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

98
99
    # n = npmat.shape[0]
    np.fill_diagonal(npmat, 0)
fabrice's avatar
fabrice committed
100
101
102
103

    '''
    Thresholding the input matrix
    '''
104
105
106
    y = stat.mquantiles(npmat[:], prob=[1 - alpha])
    th = npmat >= y
    mat_th = npmat * th
fabrice's avatar
fabrice committed
107
108
109
110
111
112
113
114
115
116

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

    '''
    ***********************************
    eigen decomposition
    '''
Fabrice  ALLAIN's avatar
Fabrice ALLAIN committed
117
    LOG.info('Decomposition and deconvolution...')
fabrice's avatar
fabrice committed
118

Fabrice Allain's avatar
Fabrice Allain committed
119
    # noinspection PyTypeChecker
fabrice's avatar
fabrice committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    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
143
        m1 = np.max(np.max(npmat * ind_nonedges))
fabrice's avatar
fabrice committed
144
145
        m2 = np.min(np.min(mat_new1))
        mat_new2 = (mat_new1 + np.max(m1 - m2, 0)) * ind_edges + (
146
            npmat * ind_nonedges)
fabrice's avatar
fabrice committed
147
148
149
150
151
152
153
154
155
156
157
158
159
    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