|
0
|
1 # ----------------------------------------------------------------------#
|
|
|
2 # Copyright (c) 2011, Richard Lupat & Jason Li.
|
|
|
3 #
|
|
|
4 # > Source License <
|
|
|
5 # This file is part of CONTRA.
|
|
|
6 #
|
|
|
7 # CONTRA is free software: you can redistribute it and/or modify
|
|
|
8 # it under the terms of the GNU General Public License as published by
|
|
|
9 # the Free Software Foundation, either version 3 of the License, or
|
|
|
10 # (at your option) any later version.
|
|
|
11 #
|
|
|
12 # CONTRA is distributed in the hope that it will be useful,
|
|
|
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
15 # GNU General Public License for more details.
|
|
|
16 #
|
|
|
17 # You should have received a copy of the GNU General Public License
|
|
|
18 # along with CONTRA. If not, see <http://www.gnu.org/licenses/>.
|
|
|
19 #
|
|
|
20 #
|
|
|
21 #-----------------------------------------------------------------------#
|
|
|
22 # Last Updated : 12 October 2011 16:43PM
|
|
|
23
|
|
|
24
|
|
|
25 def assignBin(binNumber, srcFile, binFile, targetList, minExons):
|
|
|
26 #def assignBin(binNumber, srcFile, binFile, minExons):
|
|
|
27 from math import log
|
|
|
28
|
|
|
29 src = file.readlines(open(srcFile))
|
|
|
30 #binOut = open(binFile, "w")
|
|
|
31
|
|
|
32 minExons = int(minExons)
|
|
|
33 count = 0
|
|
|
34 logcov_list = []
|
|
|
35
|
|
|
36 # Get the Log Coverage List for the normal sample
|
|
|
37 for exon in src:
|
|
|
38 exon = exon.split()
|
|
|
39 cov1 = float(exon[7])
|
|
|
40 cov2 = float(exon[8])
|
|
|
41 cov = (cov1 + cov2) / 2
|
|
|
42 if (cov > 0):
|
|
|
43 logcov = log(cov, 2)
|
|
|
44 else:
|
|
|
45 logcov = 0
|
|
|
46 logcov_list.append(logcov)
|
|
|
47 #print "Len of logcov_list:", len(logcov_list)
|
|
|
48
|
|
|
49 # Calculate the boundaries of the bins
|
|
|
50 minLog = min(logcov_list)
|
|
|
51 maxLog = max(logcov_list)
|
|
|
52 boundary = (maxLog-minLog)/binNumber
|
|
|
53 #print "Min, Max, Boundary, BinNumber: ", minLog, maxLog, boundary, binNumber
|
|
|
54
|
|
|
55
|
|
|
56 # Split exons to different bins
|
|
|
57 bin_list = []
|
|
|
58 boundary_dict = {}
|
|
|
59 for logcov in logcov_list:
|
|
|
60 i = 1
|
|
|
61 set_boundary = minLog + boundary
|
|
|
62 while (logcov > set_boundary):
|
|
|
63 i += 1
|
|
|
64 set_boundary = minLog + (boundary * i)
|
|
|
65 #boundary_dict[i] = set_boundary
|
|
|
66 bin_list.append(i)
|
|
|
67
|
|
|
68 for i in range(binNumber+2):
|
|
|
69 boundary_dict[i] = minLog + (boundary * i)
|
|
|
70
|
|
|
71
|
|
|
72 # Check the number of exons for each bin
|
|
|
73 # Merge with the adjacent bins if it is too small
|
|
|
74 for z in range(1, binNumber+1):
|
|
|
75 element = bin_list.count(z)
|
|
|
76 #print "Bin", z, "has", element, "elements"
|
|
|
77 if (element < minExons):
|
|
|
78 while (bin_list.count(z) != 0):
|
|
|
79 if (z != binNumber):
|
|
|
80 bin_list[bin_list.index(z)]+=1
|
|
|
81 else:
|
|
|
82 bin_list[bin_list.index(z)]-=1
|
|
|
83
|
|
|
84
|
|
|
85 # Check the number of exons in the last bin
|
|
|
86 last_bin_number = sorted(set(bin_list))[-1]
|
|
|
87 if len(set(bin_list)) > 1:
|
|
|
88 second_last_bin = sorted(set(bin_list))[-2]
|
|
|
89 else:
|
|
|
90 second_last_bin = last_bin_number
|
|
|
91 element = bin_list.count(last_bin_number)
|
|
|
92 if (element < minExons):
|
|
|
93 while (bin_list.count(last_bin_number) != 0):
|
|
|
94 if (last_bin_number != 1):
|
|
|
95 #bin_list[bin_list.index(last_bin_number)] -= 1
|
|
|
96 bin_list[bin_list.index(last_bin_number)] = second_last_bin
|
|
|
97
|
|
|
98 final_number_of_bin = len(set(bin_list))
|
|
|
99
|
|
|
100 # List out the binning boundaries in a txt file
|
|
|
101 boundary_list = [boundary_dict[x] for x in sorted(set(bin_list))]
|
|
|
102 i = 1
|
|
|
103
|
|
|
104 #boundary_file = binFile + str(final_number_of_bin) + ".boundaries.txt"
|
|
|
105 boundary_file = binFile + str(binNumber) + ".boundaries.txt"
|
|
|
106 boOut = open(boundary_file, "w")
|
|
|
107 boOut.write("\t".join([str(0), str(minLog)])+"\n")
|
|
|
108 for z in boundary_list:
|
|
|
109 if (i==final_number_of_bin):
|
|
|
110 z = maxLog
|
|
|
111 boOut.write("\t".join([str(i), str(z)])+"\n")
|
|
|
112 i += 1
|
|
|
113 boOut.close()
|
|
|
114
|
|
|
115 # Re-sort the bin numbers - get rid of gaps
|
|
|
116 curr_z = 1
|
|
|
117 bin_number_dict = {}
|
|
|
118 for z in sorted(set(bin_list)):
|
|
|
119 bin_number_dict[z] = curr_z
|
|
|
120 curr_z += 1
|
|
|
121
|
|
|
122
|
|
|
123 # Append the bin number to the original file
|
|
|
124 #binFile = binFile + str(final_number_of_bin)+".txt"
|
|
|
125 binFile = binFile + str(binNumber) + ".txt"
|
|
|
126 binOut = open(binFile, "w")
|
|
|
127
|
|
|
128 for exons in src:
|
|
|
129 exon = exons.split()
|
|
|
130 id = int(exon[0])
|
|
|
131 gene = exon[1]
|
|
|
132 exonNumber = exon[5]
|
|
|
133
|
|
|
134 target = targetList[int(id)-1]
|
|
|
135 if target.id == id:
|
|
|
136 chr = target.chr
|
|
|
137 oriStart = target.start
|
|
|
138 oriEnd = target.end
|
|
|
139
|
|
|
140 else:
|
|
|
141 print "ERROR..."
|
|
|
142
|
|
|
143 f_bin_number = str(bin_number_dict[bin_list[count]])
|
|
|
144 binOut.write("\t".join([exons.strip("\n"), f_bin_number,chr, oriStart, oriEnd]))
|
|
|
145 binOut.write("\n")
|
|
|
146 count += 1
|
|
|
147
|
|
|
148 binOut.close()
|
|
|
149 print "End of assign.bin.number.py with %s exons in %s bins" %(count, final_number_of_bin)
|
|
|
150
|
|
|
151
|
|
|
152 return final_number_of_bin
|
|
|
153
|
