Mercurial > repos > fcaramia > contra

comparison Contra/scripts/assign_bin_number_v2.py @ 0:7564f3b1e675

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author fcaramia
date Thu, 13 Sep 2012 02:31:43 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:7564f3b1e675
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