annotate Contra/scripts/average_count.py @ 0:7564f3b1e675

Uploaded
author fcaramia
date Thu, 13 Sep 2012 02:31:43 -0400
parents
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 : 28 Sep 2011 11:00AM
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
23
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
24 import sys
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
25 import math
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
26
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
27 def getAverage(list1):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
28 if len(list1) > 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
29 return float(sum(list1))/len(list1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
30
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
31 return 0.0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
32
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
33 def getStdDev(list1, avg):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
34 var = 0.0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
35 for x in list1:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
36 var += (avg - x) ** 2
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
37
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
38 if (len(list1)-1) > 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
39 var /= (len(list1)-1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
40
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
41 return math.sqrt(var)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
42
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
43 def getMinMax(list1):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
44 length = len(list1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
45 if length != 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
46 min = list1[0]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
47 max = list1[length-1]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
48 else:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
49 min = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
50 max = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
51
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
52 return min, max
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
53
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
54 def getMedian(list1):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
55 length = len(list1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
56 if length == 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
57 median = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
58 elif length % 2 == 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
59 median = (list1[length/2]+list1[(length/2) - 1])/2
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
60 else:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
61 median = list1[length/2]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
62 return median
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
63
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
64 def createDataDict(count, list1, r, offset, id_check, exon_check):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
65 tDict = {}
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
66 tDictOri = {}
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
67
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
68 while count < len(list1):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
69 t = list1[count].split()
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
70 tId = t[5]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
71 tExon = t[6]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
72
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
73 if (tId != id_check) or (tExon != exon_check):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
74 return count, tDict, tDictOri
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
75
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
76 tStart = int(t[2])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
77 tEnd = int(t[3])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
78 tCov = float(t[4]) / r + offset #GeoMean Normalisation
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
79 tCovOri = float(t[4]) + offset #without scaling
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
80
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
81 #filling dict
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
82 while tStart < tEnd:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
83 tDict[tStart] = tCov
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
84 tDictOri[tStart] = tCovOri #without scaling
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
85 tStart += 1
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
86
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
87 count += 1
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
88
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
89 return count, tDict, tDictOri
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
90
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
91 def getFactor (val1, val2):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
92 r = math.sqrt(val1 * val2)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
93 r1 = val1/r
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
94 r2 = val2/r
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
95 return r1, r2
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
96
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
97 def averageCount(tFile, nFile, averageOut, tReadCount, nReadCount, rd_threshold, minNBases):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
98 tList = file.readlines(open(tFile))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
99 nList = file.readlines(open(nFile))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
100 # constant & counter
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
101 OFF = 1
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
102 tCount = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
103 nCount = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
104
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
105 # create and open files
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
106 output = open(averageOut, "w")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
107
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
108 # Offset and Ratio for Geometric Mean Normalisation
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
109 r1, r2 = getFactor(tReadCount, nReadCount)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
110 if rd_threshold > 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
111 #OFF = 0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
112 OFF = 0.5
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
113
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
114 #big loop
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
115 while (nCount < len(nList)):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
116 # initialisation, get the chr, geneID, geneName
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
117 init = tList[tCount].split()
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
118 initial = init[5]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
119 _exon = init[6]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
120 chr = init[1]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
121 gene = init[0]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
122 _start = int(init[2])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
123
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
124 # check if t-gene and n-gene refer to the same gene
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
125 check_init = nList[nCount].split()
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
126 if check_init[5] != initial or check_init[6] != _exon:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
127 print "Initial: ", initial
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
128 print "Check_Init.id: ", check_init[5]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
129 print "_Exon: ", _exon
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
130 print "Check_Init.exon: ", check_init[6]
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
131 print "Error. Comparing different Gene"
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
132 sys.exit(1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
133
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
134 # create data dictionary for tumour and normal data (per each regions/ exon)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
135 tCount, tDict, tDictOri = createDataDict(tCount, tList, r1, OFF, initial, _exon)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
136 nCount, nDict, nDictOri = createDataDict(nCount, nList, r2, OFF, initial, _exon)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
137 # check number of bases in the both gene dict
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
138 if len(nDict) != len(tDict):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
139 print "N:", len(nDict)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
140 print "T:", len(tDict)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
141 print "Error. Different length of dict"
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
142 sys.exit(1)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
143
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
144 # compare coverage
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
145 count = _start
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
146 _max = max(nDict.keys())
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
147 ratioList = []
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
148 tumourList = []
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
149 normalList = []
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
150 tumourOriList = []
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
151 normalOriList = []
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
152 while count <= _max:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
153 # get ratio
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
154 if (nDict[count] < rd_threshold) and (tDict[count] < rd_threshold):
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
155 ratio = 0.0
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
156 else:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
157 if tDict[count] == 0:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
158 tDict[count] = 0.5
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
159
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
160 ratio = math.log((float(tDict[count]) / nDict[count]),2)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
161 tumourList.append(tDict[count])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
162 tumourOriList.append(tDictOri[count])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
163 normalList.append(nDict[count])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
164 normalOriList.append(nDictOri[count])
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
165 ratioList.append(ratio)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
166
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
167 count += 1
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
168
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
169 ratioLen = len(ratioList)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
170
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
171 # get average
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
172 avg = getAverage(ratioList)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
173 sd = getStdDev(ratioList, avg)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
174 tumourAvg= str(round(getAverage(tumourList),3))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
175 normalAvg= str(round(getAverage(normalList),3))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
176 tumourOriAvg = str(round(getAverage(tumourOriList),3))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
177 normalOriAvg = str(round(getAverage(normalOriList),3))
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
178
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
179 # get median
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
180 ratioList.sort()
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
181 min_logratio, max_logratio = getMinMax(ratioList)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
182 median = getMedian(ratioList)
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
183
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
184 # write output
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
185 if ratioLen >= minNBases:
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
186 output.write(initial + "\t" + gene + "\t" + str(ratioLen) + "\t")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
187 output.write(str(round(avg,3))+ "\t"+ str(count)+ "\t" + _exon + "\t")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
188 output.write(str(round(sd ,3))+ "\t"+ tumourAvg + "\t" + normalAvg +"\t")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
189 output.write(tumourOriAvg + "\t" + normalOriAvg + "\t")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
190 output.write(str(round(median,3)) + "\t" + str(round(min_logratio,3)) + "\t")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
191 output.write(str(round(max_logratio,3)) + "\n")
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
192
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
193 output.close()
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
194
7564f3b1e675 Uploaded
fcaramia
parents:
diff changeset
195 #print "End of averageCount.py with the last target = '%s'" %(initial)