Source code for matrix_permutation_test

#!/usr/bin/python

# Load required modules
import sys, os
import multi_dendrix as Multi

def parse_args(input_list=None):
    # Parse arguments
    import argparse
    class Args: pass
    args = Args()
    description = 'Calculates network permutation test for given pathway set '\
                  'on the given network.'
    parser = argparse.ArgumentParser(description=description)
    parser.add_argument('-i', '--permuted_matrices_dir', required=True,
                        help='Directory of permuted networks.')
    parser.add_argument('-p', '--pathway_set_file', required=True,
                        help='File containing input pathway set.')
    parser.add_argument('-o', '--output_file', default=None,
                        help='Name of output file.')
    parser.add_argument('-k_min', '--min_pathway_size', required=True, type=int,
                        help='Minimum pathway size.')
    parser.add_argument('-k_max', '--max_pathway_size', required=True, type=int,
                        help='Maximum pathway size.')
    parser.add_argument('-t', '--num_pathways', required=True, type=int,
                        help='Minimum number of pathways.')
    parser.add_argument('-a', '--alpha', type=float, default=1.0,
                        help='Parameter that changes weight function W by '\
                        'weighting the penalty of coverage overlap.')    
    parser.add_argument('--num_overlaps', type=int, default=0,
                        help='Number of overlaps allowed per pathway.')    
    parser.add_argument('--pathways_per_gene', type=int, default=1,
                        help='Number of pathways a gene can be a member of.')
    parser.add_argument('-v', '--verbose', default=False, action='store_true',
                        help='Flag verbose mode.')

    # If called from the command line, parse command line args.
    if input_list: parser.parse_args(input_list, namespace=args)
    else: parser.parse_args(namespace=args)
    
    return args

def load_permuted_matrices(input_dir):
	return [ Multi.load_mutation_data(input_dir + "/" + fh)
	         for fh in os.listdir(input_dir) ]

def load_w_prime(pathway_set_file):
	return sum([int(l.rstrip().split("\t")[0]) for l in open(pathway_set_file)])

def matrix_permutation_test(W_prime, permuted_matrices, t, k_min, k_max,
	                        alpha, num_overlaps, pathways_per_gene):
	count = 0.
	for m, n, genespace, patientspace, G2T, T2G in permuted_matrices:
		multi_args = [ (genespace, T2G, G2T), t, k_min, k_max, alpha,
                       num_overlaps, pathways_per_gene ]
		pathways_w_weights = Multi.multi_dendrix(*multi_args)
		if sum([w for pathway, w in pathways_w_weights]) >= W_prime:
			count += 1.
	return count / float(len(permuted_matrices))

def run(args):
	# Load input
	if args.verbose: print 'Loading matrices...'
	W_prime = load_w_prime(args.pathway_set_file)
	permuted_matrices = load_permuted_matrices(args.permuted_matrices_dir)

	# Evaluate the pathway set
	if args.verbose: print 'Evaluating input pathway set...'
	pval = matrix_permutation_test(W_prime, permuted_matrices, args.num_pathways,
		                           args.min_pathway_size, args.max_pathway_size,
		                           args.alpha, args.num_overlaps,
		                           args.pathways_per_gene)

	# Output and return results
	test_name = 'Matrix permutation pval: '
	if args.output_file:
		open(args.output_file, 'w').write(test_name + str(pval))
	
	if args.verbose: print test_name + str(pval)

if __name__ == "__main__": run(parse_args())