Merge pull request #100 from netZoo/devel

Devel

Author: marouenbg

README.md

@@ -11,7 +11,7 @@ Macos[![MAC](https://travis-ci-job-status.herokuapp.com/badge/netZoo/netZooM/mas
 netZooM is a MATLAB package of network methods.
 
 ## Zoo animals
-(gpu)PANDA, (gpu)LIONESS, PUMA, SPIDER, and optPANDA.
+(gpu)PANDA, (gpu)LIONESS, PUMA, SPIDER, optPANDA, and OTTER.
 
 ## Quick guide
 `git clone https://github.com/netZoo/netZooM.git` into your local disk. 

netZooM/lioness/lioness_run.m

@@ -112,12 +112,13 @@ function lioness_run(exp_file, motif_file, ppi_file, panda_file, save_dir,...
 if isequal(computing,'gpu')
     parpool(gpuDeviceCount); 
     parfor i = indexes
-        ExpGPU = Exp;
+        ExpGPU = gpuArray(Exp);
         fprintf('Running LIONESS for sample %d:\n', i);
         idx = [1:(i-1), (i+1):NumConditions];  % all samples except i
 
-        disp('Computing coexpresison network:');
+        disp('Computing coexpression network:');
         tic; GeneCoReg = Coexpression(ExpGPU(idx,:)); toc;
+        GeneCoReg = gather(GeneCoReg);
 
         disp('Normalizing Networks:');
         tic; GeneCoReg = NormalizeNetwork(GeneCoReg); toc;
@@ -134,7 +135,7 @@ function lioness_run(exp_file, motif_file, ppi_file, panda_file, save_dir,...
         fprintf('Running LIONESS for sample %d:\n', i);
         idx = [1:(i-1), (i+1):NumConditions];  % all samples except i
 
-        disp('Computing coexpresison network:');
+        disp('Computing coexpression network:');
         tic; GeneCoReg = Coexpression(Exp(idx,:)); toc;
 
         disp('Normalizing Networks:');

netZooM/otter/otter.m

@@ -0,0 +1,71 @@
+function W = otter(W,P,C,lambda,gamma,Imax,eta)
+% Description:
+%               OTTER infers gene regulatory networks using TF DNA binding
+%               motif (W), TF PPI (P), and gene coexpression (C) through 
+%               minimzing the following objective:
+%                                  min f(W) 
+%               with f(W) = (1-lambda)*||WW' - P||^2 + lambda*||W'W - C||^2 + (gamma/2)*||W||^2
+%
+% Inputs:
+%               W     : TF-gene regulatory network based on TF motifs as a
+%                       matrix of size (t,g), g=number of genes, t=number of TFs
+%               P     : TF-TF protein interaction network as a matrix of size (t,t)
+%               C     : gene coexpression as a matrix of size (g,g) 
+%               lambda: tuning parameter in [0,1] (higher gives more weight to C)
+%               gamma : regularization parameter
+%               Imax  : number of iterations of the algorithm
+%               eta   : learning rate
+%               bexp  : exponent influencing learning rate (higher means smaller)
+%
+%
+% Outputs:
+%               W  : Predicted TF-gene complete regulatory network as an adjacency matrix of size (t,g).
+%
+% Authors: 
+%               Rebekka Burkholz 4/2020
+
+%global parameters
+if nargin<4
+    lambda = 0.0035;
+end
+if nargin<5
+    gamma = 0.335;
+end
+if nargin<6
+    Imax = 300;
+end
+if nargin<7
+    eta = 0.00001;
+end
+%ADAM parameters
+b1 = 0.9;
+b2 = 0.999;
+eps = 0.00000001;
+%initial transformation
+C = C/trace(C);
+P = P/trace(P) + 0.0013;
+W = P*W;
+W = W/sqrt(trace(W*W'));
+
+[t, g] = size(W);
+m = zeros(t, g);
+v = m;
+b1t = b1;
+b2t = b2;
+%To save computations:
+P = -P*(1-lambda) + gamma*eye(t);
+C = -C*lambda;
+%ADAM gradient descent                  
+for i = 1:Imax
+    %gradient
+    grad = (P*W + W*C + W*W'*W);
+    m = b1*m + ((1-b1)*4)*grad;
+    v = b2*v + ((1-b2)*16)*grad.^2; 
+    b1t = b1t*b1;
+    b2t = b2t*b2;
+    alpha = sqrt(1-b2t)/(1-b1t)*eta;
+    epst = eps*sqrt((1-b2t));
+    %update of gene ragulatory matrix
+    W = W - alpha*(m./(epst+sqrt(v)));
+end
+end

netZooM/panda/panda_run.m

@@ -66,7 +66,7 @@
 %                              1 saves memory on the GPU (slower)
 %                          
 % Outputs:
-%               AgNet     : Predicted TF-gene gene complete regulatory network using PANDA as a matrix of size (t,g).
+%               AgNet     : Predicted TF-gene complete regulatory network using PANDA as a matrix of size (t,g).
 %
 % Authors: 
 %               cychen, marieke, kglass

netZooM/tools/gpuPANDA.m

@@ -209,6 +209,7 @@
     runtime = toc;
     if isequal(computing,'gpu')
         RegNet=gather(RegNet);
+        gpuDevice(1);%Clear GPU device memory 
     end
     fprintf('Running PANDA on %d Genes and %d TFs took %f seconds!\n', NumGenes, NumTFs, runtime);
 end

tests/testOtter.m

@@ -0,0 +1,51 @@
+function test_suite=testOtter()
+	initTestSuite;
+end
+
+function testOtterSimple()
+	% Tell if this is Octave (Unit tests) or Matlab
+        isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
+
+        % Load statistics package from Octave
+        if isOctave
+            %we need the nan package because it has a fast implementation of corrcoeff
+            %pkg load statistics
+            pkg load nan;
+        end
+
+        % Set Program Parameters
+        generatePriors=0;
+        if generatePriors==1
+            exp_file   = 'test_data/expression.txt';
+            motif_file = 'test_data/motifTest.txt';
+            ppi_file   = 'test_data/ppi.txt';
+            [Exp,RegNet,TFCoop,TFNames,GeneNames]=processData(exp_file,motif_file,ppi_file,'intersection');
+            tic; GeneCoReg = Coexpression(Exp); toc;
+            csvwrite('c.csv',GeneCoReg);
+            csvwrite('p.csv',TFCoop);
+            csvwrite('w.csv',RegNet);
+        elseif generatePriors==0
+            C=csvread('tests/test_data/otter/c.csv');
+            P=csvread('tests/test_data/otter/p.csv');
+            W=csvread('tests/test_data/otter/w.csv');
+        end
+
+        % Add path
+        addpath(genpath(fullfile(pwd,'tests')));
+
+        % Call Otter
+        Imax=1;
+        lambda = 0.0035;
+        gamma = 0.335;
+        tic;W = otter(W,P,C,lambda,gamma,Imax);toc;
+        
+        % Load the expected result
+        filename = 'tests/test_data/otter/test_otter.csv';
+        W_test = csvread(filename);
+
+        % Compare the outputs
+        tolMat  =1e-6;
+        deltaMat=max(max(abs(W-W_test)));
+	    assertTrue( deltaMat < tolMat );
+        
+end