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
