|
2
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import argparse
|
|
|
4
|
|
|
5 import numpy
|
|
|
6 import scipy.stats
|
|
|
7
|
|
|
8 def main():
|
|
|
9 parser = generate_parser()
|
|
|
10 args = parser.parse_args()
|
|
4
|
11 data, names = load_data(args)
|
|
2
|
12 corr = find_correlations(data, args)
|
|
|
13 save_data(corr, names, args)
|
|
|
14
|
|
|
15 def load_data(args):
|
|
|
16 infile = open(args.input)
|
|
|
17 names = []
|
|
3
|
18 cnames = None
|
|
2
|
19 data = []
|
|
|
20 if args.column:
|
|
3
|
21 cnames = infile.readline().split()
|
|
4
|
22 if args.row:
|
|
3
|
23 cnames = cnames[1:]
|
|
2
|
24 temp = infile.readline()
|
|
|
25 if args.int:
|
|
|
26 dtype = int
|
|
|
27 else:
|
|
|
28 dtype = float
|
|
|
29 while temp:
|
|
|
30 temp = temp.split()
|
|
|
31 if args.row:
|
|
|
32 names.append(temp[0])
|
|
|
33 temp = temp[1:]
|
|
|
34 data.append([])
|
|
|
35 for i in range(len(temp)):
|
|
|
36 data[-1].append(dtype(temp[i]))
|
|
|
37 temp = infile.readline()
|
|
|
38 if len(names) == 0:
|
|
|
39 names = None
|
|
|
40 data = numpy.array(data)
|
|
|
41 if args.features:
|
|
|
42 data = data.T
|
|
3
|
43 names = cnames
|
|
4
|
44 return data, names
|
|
2
|
45
|
|
|
46 def find_correlations(data, args):
|
|
|
47 corr = numpy.ones((data.shape[0], data.shape[0]), dtype=numpy.float32)
|
|
|
48 if args.test == 'pearson':
|
|
|
49 findcorr = scipy.stats.pearsonr
|
|
|
50 elif args.test == 'spearman':
|
|
|
51 findcorr = scipy.stats.spearmanr
|
|
|
52 else:
|
|
|
53 findcorr = scipy.stats.kendalltau
|
|
|
54 for i in range(data.shape[0] - 1):
|
|
|
55 for j in range(i + 1, data.shape[0]):
|
|
|
56 corr[i, j] = findcorr(data[i, :], data[j, :])[0]
|
|
|
57 corr[j, i] = corr[i, j]
|
|
|
58 return corr
|
|
|
59
|
|
|
60 def save_data(data, names, args):
|
|
|
61 output = open(args.output, 'w')
|
|
|
62 if names is not None:
|
|
|
63 output.write("%s\n" % '\t'.join(['sample'] + names))
|
|
|
64 for i in range(data.shape[0]):
|
|
|
65 if names is not None:
|
|
|
66 temp = [names[i]]
|
|
|
67 else:
|
|
|
68 temp = []
|
|
|
69 for j in range(data.shape[1]):
|
|
|
70 temp.append("%0.6f" % data[i, j])
|
|
|
71 output.write("%s\n" % '\t'.join(temp))
|
|
|
72 output.close()
|
|
|
73
|
|
|
74 def generate_parser():
|
|
|
75 """Generate an argument parser."""
|
|
|
76 description = "%(prog)s -- Create a raw file of paired aligned reads for a HiC experiment from bam files"
|
|
|
77 parser = argparse.ArgumentParser(description=description)
|
|
|
78 parser.add_argument('-f', dest="features", action='store_true', help="Rows represent features.")
|
|
|
79 parser.add_argument('-i', dest='int', action='store_true', help="Data is of type int.")
|
|
|
80 parser.add_argument('-t', dest='test', action='store', default='pearson',
|
|
|
81 choices=['spearman', 'pearson', 'kendall'], help="Type of correlation to perform.")
|
|
|
82 parser.add_argument('-r', dest='row', action='store_true', help="Row names present.")
|
|
|
83 parser.add_argument('-c', dest='column', action='store_true', help="Column names present.")
|
|
|
84 parser.add_argument(dest="input", type=str, action='store', help="Text files conatining table to be correlated.")
|
|
|
85 parser.add_argument(dest="output", type=str, action='store', help="Output destination.")
|
|
|
86 return parser
|
|
|
87
|
|
|
88 if __name__ == "__main__":
|
|
|
89 main() |
