In [11]:
import deconvolution_methods as DM
import numpy as np
from sklearn.metrics import adjusted_mutual_info_score as AMI
from scipy.stats import spearmanr,pearsonr
import pandas as pd
def bin_variable(var1): # bin with normalization
var1=np.array(var1).astype(np.float)
if abs(np.std(var1))>0.01:
var1 = (var1 - np.mean(var1))/np.std(var1)
else:
var1 = (var1 - np.mean(var1))
val1 = np.digitize(var1, np.histogram(var1, bins='fd')[1])
#print(type(val1))
#print((val1).shape())
return val1
ajd_mi_bin=lambda x,y : AMI(bin_variable(x),bin_variable(y))
In [12]:
# load data
data_in_silico=pd.read_csv('net1_expression_data.tsv',sep='\t')
print(data_in_silico)
In [13]:
P_corr=np.corrcoef(data_in_silico.transpose())
In [4]:
print(P_corr.shape)
print(len(list(data_in_silico['G1'])))
In [25]:
S_corr=np.ones((1643,1643))
AMI_fd=np.ones((1643,1643))
for idx_i,i in zip(range(len(data_in_silico.columns)),data_in_silico.columns):
for idx_j,j in zip(range(len(data_in_silico.columns)),data_in_silico.columns):
if i!=j:
S_corr[idx_i,idx_j]=spearmanr(list(data_in_silico[i]),list(data_in_silico[j]))[0]
AMI_fd[idx_i,idx_j]=ajd_mi_bin(list(data_in_silico[i]),list(data_in_silico[j]))
In [35]:
print(S_corr)
In [34]:
print(AMI_fd)
In [32]:
S_corr=np.loadtxt('S_corrM.nptxt')
np.savetxt('AMI_fd.csv',AMI_fd,delimiter=",")
In [20]:
import deconvolution_methods as DM
??DM
In [16]:
# %load ND.py
def ND(mat,beta=0.90,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. M?©dard 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 Médard and Manolis Kellis
Nature Biotechnology
--------------------------------------------------------------------------
ND.m: network deconvolution
--------------------------------------------------------------------------
DESCRIPTION:
USAGE:
mat_nd = ND(mat)
mat_nd = ND(mat,beta)
mat_nd = ND(mat,beta,alpha,control)
INPUT ARGUMENTS:
mat Input matrix, if it is a square matrix, the program assumes
it is a relevance matrix where mat(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.
alpha fraction of edges of the observed dependency matrix to be kept in
deconvolution process.
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
'''
import scipy.stats.mstats as stat
from numpy import linalg as LA
if beta>=1 or beta<=0:
print 'error: beta should be in (0,1)'
if alpha>1 or alpha<=0:
print 'error: alpha should be in (0,1)';
'''
***********************************
Processing the inut matrix
diagonal values are filtered
'''
n = mat.shape[0]
np.fill_diagonal(mat, 0)
'''
Thresholding the input matrix
'''
y =stat.mquantiles(mat[:],prob=[1-alpha])
th = mat>=y
mat_th=mat*th;
'''
making the matrix symetric if already not
'''
mat_th = (mat_th+mat_th.T)/2
'''
***********************************
eigen decomposition
'''
print 'Decomposition and deconvolution...'
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)))
'''
***********************************
displying direct weights
'''
if control==0:
ind_edges = (mat_th>0)*1.0;
ind_nonedges = (mat_th==0)*1.0;
m1 = np.max(np.max(mat*ind_nonedges));
m2 = np.min(np.min(mat_new1));
mat_new2 = (mat_new1+np.max(m1-m2,0))*ind_edges+(mat*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
In [17]:
ND_S= ND(S_corr)
ND_P= ND(P_corr)
ND_AMI=ND(AMI_fd)
In [18]:
np.savetxt('ND_S.csv',ND_S,delimiter=",")
np.savetxt('ND_P.csv',ND_P,delimiter=",")
np.savetxt('ND_AMI.csv',ND_AMI,delimiter=",")
In [ ]:
PartCorr_S= DM.deconvolution_methods(3,S_corr)
PartCorr_P= DM.deconvolution_methods(3,P_corr)
PartCorr_AMI=DM.deconvolution_methods(3,AMI_fd)
In [25]:
#Reprocessing
PartCorr_S=np.absolute(np.nan_to_num(PartCorr_S))
PartCorr_P=np.absolute(np.nan_to_num(PartCorr_P))
PartCorr_AMI=np.absolute(np.nan_to_num(PartCorr_AMI))
print(PartCorr_S)
In [26]:
np.savetxt('PartCorr_S.csv',PartCorr_S,delimiter=",")
np.savetxt('PartCorr_P.csv',PartCorr_P,delimiter=",")
np.savetxt('PartCorr_AMI.csv',PartCorr_AMI,delimiter=",")
In [28]:
data_in_silico.transpose().to_csv('data_in_silico.csv',sep=',',header=False,index=False)
In [ ]: