function [ ranked_list_genes, scores ] = netics_fun( varargin )

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  

%   INPUTS:

%   filenameMu(varargin{1}):         input file that contains the genetically aberrant genes of each sample. 

%                                    It contains two columns that map every gene (1st column) to the samples that 

%                                    it is genetically aberrant (2nd column)

%   adj(varargin{2}):                adjacency matrix of the directed interaction network

%   restart_prob(varargin{3}):       restart probability for the insulated diffusion. For (Wu et al., 2010) network use 0.4.

%   rank_method_str(varargin{4}):    'MEDIAN' uses the median of the sample-specific ranks

%                                    'RRA' uses the Robust Rank Aggregation method to integrate sample-specific ranked lists

%                                    'SUM' uses the sum of the sample-specific ranks

%   filenameNet(varargin{5}):        input file that contains the list of the genes that are present in the network. 

%                                    They should be in the same order as in the rows of the adjacency matrix adj. 

%                                    An example file is given that contains the gene names of the network described in (Wu et al, 2010).

%   filenameRNA(varargin{6}):        tab delimited file with two columns. First column contains the genes for which differential expression 

%                                    between the tumor and normal samples at the RNA level was measured. Second column contains the p-values 

%                                    of these measurements. This file can be the result of a tool for differential expression analysis such as DESeq2. 

%                                    Each gene in this file should have only one entry.

%   filenamePR(varargin{7}):         tab delimited file with two columns. First column contains the proteins for which differential expression between 

%                                    the tumor and normal samples at the protein level was measured. Second column contains the p-values of these 

%                                    measurements. Each gene in this file should have only one entry.

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  

%   OUTPUTS:

%   ranked_list_genes:  list with the ranked gene names

%   scores:             list with the scores of the ranked gene names

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  

    filenameRNA = '';

    filenamePR = '';

    ranked_list_genes={};

    % Accept arguments    

    if nargin < 5,

        error('Error: Missing input arguments.')

    elseif nargin > 7,

        error('Error: Too many input arguments.')

    else            

        mutation_data = read_mutations( varargin{1} );

        adj = varargin{2};

        restart_prob = varargin{3};

        rank_method_str = varargin{4};

        filenameNet = varargin{5};

        choose_mut_diff = 1;

    end

    if nargin >= 6,            

        filenameRNA = varargin{6};

        choose_mut_diff = 3;        

    end

    if nargin >= 7,

        filenamePR = varargin{7};

    end

    %read network genes

    fidnet = fopen(filenameNet,'r');

    g = textscan(fidnet, '%s', 'delimiter', '\n');

    network_genes = g{1};

    fclose(fidnet);

    %Extract first token

    rank_method_str1 = rank_method_str(1:strfind(rank_method_str,'='));

    %Handle input

    rank_method_str = rank_method_str(strfind(rank_method_str,'=')+1:end);

    %Check input - Token 1

    if ~strcmp(rank_method_str1,'RANK_AGGREG='),

        disp(strcat('Wrong input: ',rank_method_str1));

        disp('Input should be RANK_AGGREG=');       

        return; 

    end

    %Check input - Token 2

    if  ~strcmp(rank_method_str,'SUM') && ~strcmp(rank_method_str,'MEDIAN') && ~strcmp(rank_method_str,'RRA'),

        disp(strcat('Wrong input: ',rank_method_str));

        disp(strcat('Input should be MEDIAN, RRA or SUM'));

        return;     

    end

    %copy data from input mutation_data 

    samples = {};

    for i = 1:length(mutation_data),

        samples.DNA{i} = upper(mutation_data{i});

        samples.DNA_network{i} = upper(intersect(network_genes,samples.DNA{i}));

    end

    clear mutation_data;

    %read differentially expressed genes

    diff_expr_genes = {};

    if choose_mut_diff~=1,

        diff_expr_genes = read_diff_expr( filenameRNA, filenamePR );

        if isempty(diff_expr_genes),

            choose_mut_diff=1;

        end

    else

        disp('No Input files were provided for differentially expressed genes.');        

    end

    for i = 1:length(samples.DNA),

        samples.RNA{i} = diff_expr_genes;

        samples.RNA_network{i} = upper(intersect(network_genes,samples.RNA{i}));

    end

    %compute diffused matrix

    disp('Computing diffused matrix...');

    F = insulated_diff(norm_adj(adj), restart_prob);

    F_opp = insulated_diff(norm_adj(adj'), restart_prob);    

    disp('Running NetICS...');

    [ranked_list_genes, scores] = prioritization( samples, F, F_opp, network_genes, choose_mut_diff, rank_method_str);

end

function [ mutation_data ] = read_mutations( filename )

    fid = fopen(filename,'r');

    g = textscan(fid, '%s%s', 'delimiter', '\t');

    fclose(fid);

    unique_samples = unique(g{2});

    mutation_data = cell(length(unique_samples),1);

    for i = 1:length(unique_samples),

        mutation_data{i} = g{1}(ismember(g{2},unique_samples{i}));

    end

end

function [ diff_expr_genes ] = read_diff_expr( filenameDE, filenamePR )

    diff_expr_genes ={}; RNA_names={}; rppa_names={};

    if ~isempty(filenamePR),

        [rppa_names, rppa_pval] = read_file(filenamePR);

    else

        disp('No Input file was provided for differentially expressed genes at the proteome level.');        

    end

    if ~isempty(filenameDE),

        [RNA_names, RNA_pval] = read_file(filenameDE);        

    else

        disp('No Input file was provided for differentially expressed genes at the RNA level.');        

    end

    %check for duplicate genes in the input files

    if (length(unique(RNA_names)) ~= length(RNA_names)) || (length(unique(rppa_names)) ~= length(rppa_names)),

        disp('Input files for differentially expressed genes should contain only one entry per gene.');

        disp('Input files for differentially expressed genes were ignored.');        

        return;         

    end

    %check if the input files for differential expression are empty    

    if isempty(rppa_names) && isempty(RNA_names),

        disp('Only genetically aberrant genes are used for diffusion.');

        return;

    %only differential expressed genes at the RNA level are provided

    elseif isempty(rppa_names) && ~isempty(RNA_names),

        [~,~,pval_all] = fdr(RNA_pval,0.05);

        diff_expr_genes = RNA_names(pval_all < 0.05);

    %only differential expressed genes at the protein level are provided        

    elseif isempty(RNA_names) && ~isempty(rppa_names),

        [~,~,pval_all] = fdr(rppa_pval,0.05);

        diff_expr_genes = rppa_names(pval_all < 0.05);

    else

        names_intersection = intersect(RNA_names,rppa_names);

        names_setdiff_rna = setdiff(RNA_names,rppa_names);

        names_setdiff_pr = setdiff(rppa_names,RNA_names);

        rna_pval_names_intersection = RNA_pval(ismember(RNA_names,names_intersection));

        rppa_pval_names_intersection = rppa_pval(ismember(rppa_names,names_intersection));

        rna_pval_names_setdiff_rna = RNA_pval(ismember(RNA_names,names_setdiff_rna));

        pr_pval_names_setdiff_pr = rppa_pval(ismember(rppa_names,names_setdiff_pr));

        all_names = [names_intersection' names_setdiff_rna' names_setdiff_pr'];

        %Fishers method

        pval_all = (1-pchisq(-2*(log(rna_pval_names_intersection)+log(rppa_pval_names_intersection)),4));

        pval_all = [pval_all' rna_pval_names_setdiff_rna' pr_pval_names_setdiff_pr'];    

        [~,~,pval_all] = fdr(pval_all,0.05);

        diff_expr_genes = all_names(pval_all < 0.05);

    end

end

function [names, pval] = read_file(filenamePR)

    fid = fopen(filenamePR,'r');

    g = textscan(fid, '%s%s', 'delimiter', '\t');

    fclose(fid);

    names = g{1}; pval = str2double(g{2}); 

end