annotate Contra/scripts/assign_bin_number_v2.py @ 11:b6ce842a633d

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