Mercurial > repos > mheinzl > fsd

comparison fsd.py @ 17:2e517a54eedc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/monikaheinzl/duplexanalysis_galaxy/tree/master/tools/fsd commit b8a2f7b7615b2bcd3b602027af31f4e677da94f6-dirty
author mheinzl
date Tue, 02 Apr 2019 05:10:09 -0400
parents 6bd9ef49d013
children c825a29a7d9f
comparison
equal deleted inserted replaced
16:6bd9ef49d013 17:2e517a54eedc
43 parser.add_argument('--output_tabular', default="data.tabular", type=str, help='Name of the tabular file.') 43 parser.add_argument('--output_tabular', default="data.tabular", type=str, help='Name of the tabular file.')
44 return parser 44 return parser
45 45
46 46
47 def compare_read_families(argv): 47 def compare_read_families(argv):
48
48 parser = make_argparser() 49 parser = make_argparser()
49 args = parser.parse_args(argv[1:]) 50 args = parser.parse_args(argv[1:])
50
51 firstFile = args.inputFile1 51 firstFile = args.inputFile1
52 name1 = args.inputName1 52 name1 = args.inputName1
53
54 secondFile = args.inputFile2 53 secondFile = args.inputFile2
55 name2 = args.inputName2 54 name2 = args.inputName2
56 thirdFile = args.inputFile3 55 thirdFile = args.inputFile3
57 name3 = args.inputName3 56 name3 = args.inputName3
58 fourthFile = args.inputFile4 57 fourthFile = args.inputFile4
59 name4 = args.inputName4 58 name4 = args.inputName4
60
61 title_file = args.output_tabular 59 title_file = args.output_tabular
62 title_file2 = args.output_pdf 60 title_file2 = args.output_pdf
63 61
64 sep = "\t" 62 sep = "\t"
65 63
88 list_to_plot.append(data1) 86 list_to_plot.append(data1)
89 label.append(name1) 87 label.append(name1)
90 data_array_list.append(file1) 88 data_array_list.append(file1)
91 89
92 legend = "\n\n\n{}".format(name1) 90 legend = "\n\n\n{}".format(name1)
93 plt.text(0.1, 0.11, legend, size=12, transform=plt.gcf().transFigure) 91 plt.text(0.05, 0.11, legend, size=10, transform=plt.gcf().transFigure)
94 legend1 = "singletons:\nabsolute nr.\n{:,}".format(numpy.bincount(data1)[1]) 92 legend1 = "singletons:\nnr. of tags\n{:,}".format(numpy.bincount(data1)[1])
95 plt.text(0.4, 0.11, legend1, size=12, transform=plt.gcf().transFigure) 93 plt.text(0.32, 0.11, legend1, size=10, transform=plt.gcf().transFigure)
96 94
97 legend3 = "rel. freq\n{:.3f}".format(float(numpy.bincount(data1)[1]) / len(data1)) 95 legend3 = "freq. of tags\n{:.3f}".format(float(numpy.bincount(data1)[1]) / len(data1))
98 plt.text(0.5, 0.11, legend3, size=12, transform=plt.gcf().transFigure) 96 plt.text(0.41, 0.11, legend3, size=10, transform=plt.gcf().transFigure)
99 97
100 legend4 = "family size > 20:\nabsolute nr.\n{:,}".format(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1].astype(int)) 98 legend3b = "PE reads\n{:.3f}".format(float(numpy.bincount(data1)[1]) / sum(integers))
101 plt.text(0.6, 0.11, legend4, size=12, transform=plt.gcf().transFigure) 99 plt.text(0.5, 0.11, legend3b, size=10, transform=plt.gcf().transFigure)
102 100
103 legend5 = "rel. freq\n{:.3f}".format(float(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1]) / len(data1)) 101 legend4 = "family size > 20:\nnr. of tags\n{:,} ({:.3f})".format(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1].astype(int), float(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1]) / len(data1))
104 plt.text(0.7, 0.11, legend5, size=12, transform=plt.gcf().transFigure) 102 plt.text(0.58, 0.11, legend4, size=10, transform=plt.gcf().transFigure)
105 103
106 legend6 = "total length\n{:,}".format(len(data1)) 104 legend5 = "PE reads\n{:,} ({:.3f})".format(sum(integers[integers > 20]), float(sum(integers[integers > 20])) / sum(integers))
107 plt.text(0.8, 0.11, legend6, size=12, transform=plt.gcf().transFigure) 105 plt.text(0.70, 0.11, legend5, size=10, transform=plt.gcf().transFigure)
106
107 legend6 = "total nr. of\ntags\n{:,}".format(len(data1))
108 plt.text(0.82, 0.11, legend6, size=10, transform=plt.gcf().transFigure)
109
110 legend6b = "PE reads\n{:,}".format(sum(integers))
111 plt.text(0.89, 0.11, legend6b, size=10, transform=plt.gcf().transFigure)
108 112
109 if secondFile != str(None): 113 if secondFile != str(None):
110 file2 = readFileReferenceFree(secondFile) 114 file2 = readFileReferenceFree(secondFile)
115 integers2 = numpy.array(file2[:, 0]).astype(int) # keep original family sizes
116
111 data2 = numpy.asarray(file2[:, 0]).astype(int) 117 data2 = numpy.asarray(file2[:, 0]).astype(int)
112 bigFamilies2 = numpy.where(data2 > 20)[0] 118 bigFamilies2 = numpy.where(data2 > 20)[0]
113 data2[bigFamilies2] = 22 119 data2[bigFamilies2] = 22
114 120
115 list_to_plot.append(data2) 121 list_to_plot.append(data2)
116 name2 = name2.split(".tabular")[0] 122 name2 = name2.split(".tabular")[0]
117 label.append(name2) 123 label.append(name2)
118 data_array_list.append(file2) 124 data_array_list.append(file2)
119 125
120 plt.text(0.1, 0.09, name2, size=12, transform=plt.gcf().transFigure) 126 plt.text(0.05, 0.09, name2, size=10, transform=plt.gcf().transFigure)
121 127
122 legend1 = "{:,}".format(numpy.bincount(data2)[1]) 128 legend1 = "{:,}".format(numpy.bincount(data2)[1])
123 plt.text(0.4, 0.09, legend1, size=12, transform=plt.gcf().transFigure) 129 plt.text(0.32, 0.09, legend1, size=10, transform=plt.gcf().transFigure)
124 130
125 legend3 = "{:.3f}".format(float(numpy.bincount(data2)[1]) / len(data2)) 131 legend3 = "{:.3f}".format(float(numpy.bincount(data2)[1]) / len(data2))
126 plt.text(0.5, 0.09, legend3, size=12, transform=plt.gcf().transFigure) 132 plt.text(0.41, 0.09, legend3, size=10, transform=plt.gcf().transFigure)
127 133
128 legend4 = "{:,}".format(numpy.bincount(data2)[len(numpy.bincount(data2)) - 1].astype(int)) 134 legend3b = "{:.3f}".format(float(numpy.bincount(data2)[1]) / sum(integers2))
129 plt.text(0.6, 0.09, legend4, size=12, transform=plt.gcf().transFigure) 135 plt.text(0.5, 0.09, legend3b, size=10, transform=plt.gcf().transFigure)
130 136
131 legend5 = "{:.3f}".format(float(numpy.bincount(data2)[len(numpy.bincount(data2)) - 1]) / len(data2)) 137 legend4 = "{:,} ({:.3f})".format(
132 plt.text(0.7, 0.09, legend5, size=12, transform=plt.gcf().transFigure) 138 numpy.bincount(data2)[len(numpy.bincount(data2)) - 1].astype(int),
139 float(numpy.bincount(data2)[len(numpy.bincount(data2)) - 1]) / len(data2))
140 plt.text(0.58, 0.09, legend4, size=10, transform=plt.gcf().transFigure)
141
142 legend5 = "{:,} ({:.3f})".format(sum(integers2[integers2 > 20]), float(sum(integers2[integers2 > 20])) / sum(integers2))
143 plt.text(0.70, 0.09, legend5, size=10, transform=plt.gcf().transFigure)
133 144
134 legend6 = "{:,}".format(len(data2)) 145 legend6 = "{:,}".format(len(data2))
135 plt.text(0.8, 0.09, legend6, size=12, transform=plt.gcf().transFigure) 146 plt.text(0.82, 0.09, legend6, size=10, transform=plt.gcf().transFigure)
147
148 legend6b = "{:,}".format(sum(integers2))
149 plt.text(0.89, 0.09, legend6b, size=10, transform=plt.gcf().transFigure)
136 150
137 if thirdFile != str(None): 151 if thirdFile != str(None):
138 file3 = readFileReferenceFree(thirdFile) 152 file3 = readFileReferenceFree(thirdFile)
153 integers3 = numpy.array(file3[:, 0]).astype(int) # keep original family sizes
139 154
140 data3 = numpy.asarray(file3[:, 0]).astype(int) 155 data3 = numpy.asarray(file3[:, 0]).astype(int)
141 bigFamilies3 = numpy.where(data3 > 20)[0] 156 bigFamilies3 = numpy.where(data3 > 20)[0]
142 data3[bigFamilies3] = 22 157 data3[bigFamilies3] = 22
143 158
144 list_to_plot.append(data3) 159 list_to_plot.append(data3)
145 name3 = name3.split(".tabular")[0] 160 name3 = name3.split(".tabular")[0]
146 label.append(name3) 161 label.append(name3)
147 data_array_list.append(file3) 162 data_array_list.append(file3)
148 163
149 plt.text(0.1, 0.07, name3, size=12, transform=plt.gcf().transFigure) 164 plt.text(0.05, 0.07, name3, size=10, transform=plt.gcf().transFigure)
150 165
151 legend1 = "{:,}".format(numpy.bincount(data3)[1]) 166 legend1 = "{:,}".format(numpy.bincount(data3)[1])
152 plt.text(0.4, 0.07, legend1, size=12, transform=plt.gcf().transFigure) 167 plt.text(0.32, 0.07, legend1, size=10, transform=plt.gcf().transFigure)
153 168
154 legend3 = "{:.3f}".format(float(numpy.bincount(data3)[1]) / len(data3)) 169 legend3 = "{:.3f}".format(float(numpy.bincount(data3)[1]) / len(data3))
155 plt.text(0.5, 0.07, legend3, size=12, transform=plt.gcf().transFigure) 170 plt.text(0.41, 0.07, legend3, size=10, transform=plt.gcf().transFigure)
156 171
157 legend4 = "{:,}".format(numpy.bincount(data3)[len(numpy.bincount(data3)) - 1].astype(int)) 172 legend3b = "{:.3f}".format(float(numpy.bincount(data3)[1]) / sum(integers3))
158 plt.text(0.6, 0.07, legend4, size=12, transform=plt.gcf().transFigure) 173 plt.text(0.5, 0.07, legend3b, size=10, transform=plt.gcf().transFigure)
159 174
160 legend5 = "{:.3f}".format(float(numpy.bincount(data3)[len(numpy.bincount(data3)) - 1]) / len(data3)) 175 legend4 = "{:,} ({:.3f})".format(
161 plt.text(0.7, 0.07, legend5, size=12, transform=plt.gcf().transFigure) 176 numpy.bincount(data3)[len(numpy.bincount(data3)) - 1].astype(int),
177 float(numpy.bincount(data3)[len(numpy.bincount(data3)) - 1]) / len(data3))
178 plt.text(0.58, 0.07, legend4, size=10, transform=plt.gcf().transFigure)
179
180 legend5 = "{:,} ({:.3f})".format(sum(integers3[integers3 > 20]),
181 float(sum(integers3[integers3 > 20])) / sum(integers3))
182 plt.text(0.70, 0.07, legend5, size=10, transform=plt.gcf().transFigure)
162 183
163 legend6 = "{:,}".format(len(data3)) 184 legend6 = "{:,}".format(len(data3))
164 plt.text(0.8, 0.07, legend6, size=12, transform=plt.gcf().transFigure) 185 plt.text(0.82, 0.07, legend6, size=10, transform=plt.gcf().transFigure)
186
187 legend6b = "{:,}".format(sum(integers3))
188 plt.text(0.89, 0.07, legend6b, size=10, transform=plt.gcf().transFigure)
165 189
166 if fourthFile != str(None): 190 if fourthFile != str(None):
167 file4 = readFileReferenceFree(fourthFile) 191 file4 = readFileReferenceFree(fourthFile)
192 integers4 = numpy.array(file4[:, 0]).astype(int) # keep original family sizes
168 193
169 data4 = numpy.asarray(file4[:, 0]).astype(int) 194 data4 = numpy.asarray(file4[:, 0]).astype(int)
170 195
171 bigFamilies4 = numpy.where(data4 > 20)[0] 196 bigFamilies4 = numpy.where(data4 > 20)[0]
172 data4[bigFamilies4] = 22 197 data4[bigFamilies4] = 22
174 list_to_plot.append(data4) 199 list_to_plot.append(data4)
175 name4 = name4.split(".tabular")[0] 200 name4 = name4.split(".tabular")[0]
176 label.append(name4) 201 label.append(name4)
177 data_array_list.append(file4) 202 data_array_list.append(file4)
178 203
179 plt.text(0.1, 0.05, name4, size=12, transform=plt.gcf().transFigure) 204 plt.text(0.05, 0.05, name4, size=10, transform=plt.gcf().transFigure)
180 205
181 legend1 = "{:,}".format(numpy.bincount(data4)[1]) 206 legend1 = "{:,}".format(numpy.bincount(data4)[1])
182 plt.text(0.4, 0.05, legend1, size=12, transform=plt.gcf().transFigure) 207 plt.text(0.32, 0.05, legend1, size=10, transform=plt.gcf().transFigure)
183 208
184 legend4 = "{:.3f}".format(float(numpy.bincount(data4)[1]) / len(data4)) 209 legend3 = "{:.3f}".format(float(numpy.bincount(data4)[1]) / len(data4))
185 plt.text(0.5, 0.05, legend4, size=12, transform=plt.gcf().transFigure) 210 plt.text(0.41, 0.05, legend3, size=10, transform=plt.gcf().transFigure)
186 211
187 legend4 = "{:,}".format(numpy.bincount(data4)[len(numpy.bincount(data4)) - 1].astype(int)) 212 legend3b = "{:.3f}".format(float(numpy.bincount(data4)[1]) / sum(integers4))
188 plt.text(0.6, 0.05, legend4, size=12, transform=plt.gcf().transFigure) 213 plt.text(0.5, 0.05, legend3b, size=10, transform=plt.gcf().transFigure)
189 214
190 legend5 = "{:.3f}".format(float(numpy.bincount(data4)[len(numpy.bincount(data4)) - 1]) / len(data4)) 215 legend4 = "{:,} ({:.3f})".format(
191 plt.text(0.7, 0.05, legend5, size=12, transform=plt.gcf().transFigure) 216 numpy.bincount(data4)[len(numpy.bincount(data4)) - 1].astype(int),
217 float(numpy.bincount(data4)[len(numpy.bincount(data4)) - 1]) / len(data4))
218 plt.text(0.58, 0.05, legend4, size=10, transform=plt.gcf().transFigure)
219
220 legend5 = "{:,} ({:.3f})".format(sum(integers4[integers4 > 20]),
221 float(sum(integers4[integers4 > 20])) / sum(integers4))
222 plt.text(0.70, 0.05, legend5, size=10, transform=plt.gcf().transFigure)
192 223
193 legend6 = "{:,}".format(len(data4)) 224 legend6 = "{:,}".format(len(data4))
194 plt.text(0.8, 0.05, legend6, size=12, transform=plt.gcf().transFigure) 225 plt.text(0.82, 0.05, legend6, size=10, transform=plt.gcf().transFigure)
226
227 legend6b = "{:,}".format(sum(integers4))
228 plt.text(0.89, 0.05, legend6b, size=10, transform=plt.gcf().transFigure)
195 229
196 maximumX = numpy.amax(numpy.concatenate(list_to_plot)) 230 maximumX = numpy.amax(numpy.concatenate(list_to_plot))
197 minimumX = numpy.amin(numpy.concatenate(list_to_plot)) 231 minimumX = numpy.amin(numpy.concatenate(list_to_plot))
198 232
199 counts = plt.hist(list_to_plot, bins=range(minimumX, maximumX + 1), stacked=False, edgecolor="black", 233 counts = plt.hist(list_to_plot, bins=range(minimumX, maximumX + 1), stacked=False, edgecolor="black",
200 linewidth=1, label=label, align="left", alpha=0.7, rwidth=0.8) 234 linewidth=1, label=label, align="left", rwidth=0.8, alpha=0.7)
201 235
202 ticks = numpy.arange(minimumX - 1, maximumX, 1) 236 ticks = numpy.arange(minimumX - 1, maximumX, 1)
203 ticks1 = map(str, ticks) 237 ticks1 = map(str, ticks)
204 ticks1[len(ticks1) - 1] = ">20" 238 ticks1[len(ticks1) - 1] = ">20"
205 plt.xticks(numpy.array(ticks), ticks1) 239 plt.xticks(numpy.array(ticks), ticks1)
240 else: 274 else:
241 for i in counts[0]: 275 for i in counts[0]:
242 output_file.write("{}{}".format(int(sum(i)), sep)) 276 output_file.write("{}{}".format(int(sum(i)), sep))
243 277
244 # Family size distribution after DCS and SSCS 278 # Family size distribution after DCS and SSCS
245 for dataset, data, name_file in zip(list_to_plot, data_array_list, label): 279 for dataset, data_o, name_file in zip(list_to_plot, data_array_list, label):
246 maximumX = numpy.amax(dataset) 280 maximumX = numpy.amax(dataset)
247 minimumX = numpy.amin(dataset) 281 minimumX = numpy.amin(dataset)
248 282
249 tags = numpy.array(data[:, 2]) 283 tags = numpy.array(data_o[:, 2])
250 seq = numpy.array(data[:, 1]) 284 seq = numpy.array(data_o[:, 1])
251 data = numpy.array(dataset) 285 data = numpy.array(dataset)
252 286 data_o = numpy.array(data_o[:, 0]).astype(int)
253 # find all unique tags and get the indices for ALL tags, but only once 287 # find all unique tags and get the indices for ALL tags, but only once
254 u, index_unique, c = numpy.unique(numpy.array(seq), return_counts=True, return_index=True) 288 u, index_unique, c = numpy.unique(numpy.array(seq), return_counts=True, return_index=True)
255 d = u[c > 1] 289 d = u[c > 1]
256 290
257 # get family sizes, tag for duplicates 291 # get family sizes, tag for duplicates
258 duplTags_double = data[numpy.in1d(seq, d)] 292 duplTags_double = data[numpy.in1d(seq, d)]
293 duplTags_double_o = data_o[numpy.in1d(seq, d)]
294
259 duplTags = duplTags_double[0::2] # ab of DCS 295 duplTags = duplTags_double[0::2] # ab of DCS
296 duplTags_o = duplTags_double_o[0::2] # ab of DCS
297
260 duplTagsBA = duplTags_double[1::2] # ba of DCS 298 duplTagsBA = duplTags_double[1::2] # ba of DCS
261 299 duplTagsBA_o = duplTags_double_o[1::2] # ba of DCS
262 # duplTags_double_tag = tags[numpy.in1d(seq, d)]
263 # duplTags_double_seq = seq[numpy.in1d(seq, d)]
264 300
265 # get family sizes for SSCS with no partner 301 # get family sizes for SSCS with no partner
266 ab = numpy.where(tags == "ab")[0] 302 ab = numpy.where(tags == "ab")[0]
267 abSeq = seq[ab] 303 abSeq = seq[ab]
304 ab_o = data_o[ab]
268 ab = data[ab] 305 ab = data[ab]
306
269 ba = numpy.where(tags == "ba")[0] 307 ba = numpy.where(tags == "ba")[0]
270 baSeq = seq[ba] 308 baSeq = seq[ba]
309 ba_o = data_o[ba]
271 ba = data[ba] 310 ba = data[ba]
272 311
273 dataAB = ab[numpy.in1d(abSeq, d, invert=True)] 312 dataAB = ab[numpy.in1d(abSeq, d, invert=True)]
313 dataAB_o = ab_o[numpy.in1d(abSeq, d, invert=True)]
314
274 dataBA = ba[numpy.in1d(baSeq, d, invert=True)] 315 dataBA = ba[numpy.in1d(baSeq, d, invert=True)]
316 dataBA_o = ba_o[numpy.in1d(baSeq, d, invert=True)]
275 317
276 list1 = [duplTags_double, dataAB, dataBA] # list for plotting 318 list1 = [duplTags_double, dataAB, dataBA] # list for plotting
277 319
278 # information for family size >= 3 320 # information for family size >= 3
279 dataAB_FS3 = dataAB[dataAB >= 3] 321 dataAB_FS3 = dataAB[dataAB >= 3]
322 dataAB_FS3_o = dataAB_o[dataAB_o >= 3]
280 dataBA_FS3 = dataBA[dataBA >= 3] 323 dataBA_FS3 = dataBA[dataBA >= 3]
324 dataBA_FS3_o = dataBA_o[dataBA_o >= 3]
281 ab_FS3 = ab[ab >= 3] 325 ab_FS3 = ab[ab >= 3]
282 ba_FS3 = ba[ba >= 3] 326 ba_FS3 = ba[ba >= 3]
327 ab_FS3_o = ab_o[ab_o >= 3]
328 ba_FS3_o = ba_o[ba_o >= 3]
283 329
284 duplTags_FS3 = duplTags[(duplTags >= 3) & (duplTagsBA >= 3)] # ab+ba with FS>=3 330 duplTags_FS3 = duplTags[(duplTags >= 3) & (duplTagsBA >= 3)] # ab+ba with FS>=3
285 duplTags_FS3_BA = duplTagsBA[(duplTags >= 3) & (duplTagsBA >= 3)] # ba+ab with FS>=3 331 duplTags_FS3_BA = duplTagsBA[(duplTags >= 3) & (duplTagsBA >= 3)] # ba+ab with FS>=3
286 duplTags_double_FS3 = len(duplTags_FS3) + len(duplTags_FS3_BA) # both ab and ba strands with FS>=3 332 duplTags_double_FS3 = len(duplTags_FS3) + len(duplTags_FS3_BA) # both ab and ba strands with FS>=3
287 333
334 # original FS
335 duplTags_FS3_o = duplTags_o[(duplTags_o >= 3) & (duplTagsBA_o >= 3)] # ab+ba with FS>=3
336 duplTags_FS3_BA_o = duplTagsBA_o[(duplTags_o >= 3) & (duplTagsBA_o >= 3)] # ba+ab with FS>=3
337 duplTags_double_FS3_o = sum(duplTags_FS3_o) + sum(duplTags_FS3_BA_o) # both ab and ba strands with FS>=3
338
288 fig = plt.figure() 339 fig = plt.figure()
289
290 plt.subplots_adjust(bottom=0.3) 340 plt.subplots_adjust(bottom=0.3)
291 counts = plt.hist(list1, bins=range(minimumX, maximumX + 1), stacked=True, label=["duplex", "ab", "ba"], edgecolor="black", linewidth=1, align="left", color=["#FF0000", "#5FB404", "#FFBF00"]) 341 counts = plt.hist(list1, bins=range(minimumX, maximumX + 1), stacked=True, label=["duplex", "ab", "ba"],
342 edgecolor="black", linewidth=1, align="left", color=["#FF0000", "#5FB404", "#FFBF00"],
343 rwidth=0.8)
292 # tick labels of x axis 344 # tick labels of x axis
293 ticks = numpy.arange(minimumX - 1, maximumX, 1) 345 ticks = numpy.arange(minimumX - 1, maximumX, 1)
294 ticks1 = map(str, ticks) 346 ticks1 = map(str, ticks)
295 ticks1[len(ticks1) - 1] = ">20" 347 ticks1[len(ticks1) - 1] = ">20"
296 plt.xticks(numpy.array(ticks), ticks1) 348 plt.xticks(numpy.array(ticks), ticks1)
297 singl = counts[0][2][0] # singletons 349 singl = counts[0][2][0] # singletons
298 last = counts[0][2][len(counts[0][0]) - 1] # large families 350 last = counts[0][2][len(counts[0][0]) - 1] # large families
299 351
300 plt.legend(loc='upper right', fontsize=14, bbox_to_anchor=(0.9, 1), frameon=True) 352 plt.legend(loc='upper right', fontsize=14, bbox_to_anchor=(0.9, 1), frameon=True)
301 # plt.title(name1, fontsize=14) 353 plt.title(name_file, fontsize=14)
302 plt.xlabel("Family size", fontsize=14) 354 plt.xlabel("Family size", fontsize=14)
303 plt.ylabel("Absolute Frequency", fontsize=14) 355 plt.ylabel("Absolute Frequency", fontsize=14)
304 plt.margins(0.01, None) 356 plt.margins(0.01, None)
305 plt.grid(b=True, which="major", color="#424242", linestyle=":") 357 plt.grid(b=True, which="major", color="#424242", linestyle=":")
306 358
307 # extra information beneath the plot 359 # extra information beneath the plot
308 legend = "SSCS ab= \nSSCS ba= \nDCS (total)= \nlength of dataset=" 360 legend = "SSCS ab= \nSSCS ba= \nDCS (total)= \ntotal nr. of tags="
309 plt.text(0.1, 0.09, legend, size=12, transform=plt.gcf().transFigure) 361 plt.text(0.1, 0.09, legend, size=10, transform=plt.gcf().transFigure)
310 362
311 legend = "absolute numbers\n\n{:,}\n{:,}\n{:,} ({:,})\n{:,}".format(len(dataAB), len(dataBA), len(duplTags), len(duplTags_double), (len(dataAB) + len(dataBA) + len(duplTags))) 363 legend = "nr. of tags\n\n{:,}\n{:,}\n{:,} ({:,})\n{:,}".format(len(dataAB), len(dataBA), len(duplTags), len(duplTags_double), (len(dataAB) + len(dataBA) + len(duplTags)))
312 plt.text(0.35, 0.09, legend, size=12, transform=plt.gcf().transFigure) 364 plt.text(0.23, 0.09, legend, size=10, transform=plt.gcf().transFigure)
313 365
314 legend = "relative frequencies\nunique\n{:.3f}\n{:.3f}\n{:.3f}\n{:,}".format(float(len(dataAB)) / (len(dataAB) + len(dataBA) + len(duplTags)), float(len(dataBA)) / (len(dataAB) + len(dataBA) + len(duplTags)), float(len(duplTags)) / (len(dataAB) + len(dataBA) + len(duplTags)), (len(dataAB) + len(dataBA) + len(duplTags))) 366 legend5 = "PE reads\n\n{:,}\n{:,}\n{:,} ({:,})\n{:,}".format(sum(dataAB_o), sum(dataBA_o), sum(duplTags_o), sum(duplTags_double_o), (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)))
315 plt.text(0.54, 0.09, legend, size=12, transform=plt.gcf().transFigure) 367 plt.text(0.38, 0.09, legend5, size=10, transform=plt.gcf().transFigure)
368
369 legend = "rel. freq. of tags\nunique\n{:.3f}\n{:.3f}\n{:.3f}\n{:,}".format(float(len(dataAB)) / (len(dataAB) + len(dataBA) + len(duplTags)), float(len(dataBA)) / (len(dataAB) + len(dataBA) + len(duplTags)), float(len(duplTags)) / (len(dataAB) + len(dataBA) + len(duplTags)), (len(dataAB) + len(dataBA) + len(duplTags)))
370 plt.text(0.54, 0.09, legend, size=10, transform=plt.gcf().transFigure)
316 371
317 legend = "total\n{:.3f}\n{:.3f}\n{:.3f} ({:.3f})\n{:,}".format(float(len(dataAB)) / (len(ab) + len(ba)), float(len(dataBA)) / (len(ab) + len(ba)), float(len(duplTags)) / (len(ab) + len(ba)), float(len(duplTags_double)) / (len(ab) + len(ba)), (len(ab) + len(ba))) 372 legend = "total\n{:.3f}\n{:.3f}\n{:.3f} ({:.3f})\n{:,}".format(float(len(dataAB)) / (len(ab) + len(ba)), float(len(dataBA)) / (len(ab) + len(ba)), float(len(duplTags)) / (len(ab) + len(ba)), float(len(duplTags_double)) / (len(ab) + len(ba)), (len(ab) + len(ba)))
318 plt.text(0.64, 0.09, legend, size=12, transform=plt.gcf().transFigure) 373 plt.text(0.64, 0.09, legend, size=10, transform=plt.gcf().transFigure)
319 374
320 legend1 = "\nsingletons:\nfamily size > 20:" 375 legend1 = "\nsingletons:\nfamily size > 20:"
321 plt.text(0.1, 0.03, legend1, size=12, transform=plt.gcf().transFigure) 376 plt.text(0.1, 0.03, legend1, size=10, transform=plt.gcf().transFigure)
322 377
323 legend4 = "{:,}\n{:,}".format(singl.astype(int), last.astype(int)) 378 legend4 = "{:,}\n{:,}".format(singl.astype(int), last.astype(int))
324 plt.text(0.35, 0.03, legend4, size=12, transform=plt.gcf().transFigure) 379 plt.text(0.23, 0.03, legend4, size=10, transform=plt.gcf().transFigure)
325 380
326 legend3 = "{:.3f}\n{:.3f}".format(singl / len(data), last / len(data)) 381 legend3 = "{:.3f}\n{:.3f}".format(singl / len(data), last / len(data))
327 plt.text(0.54, 0.03, legend3, size=12, transform=plt.gcf().transFigure) 382 plt.text(0.64, 0.03, legend3, size=10, transform=plt.gcf().transFigure)
383
384 legend3 = "\n\n{:,}".format(sum(data_o[data_o > 20]))
385 plt.text(0.38, 0.03, legend3, size=10, transform=plt.gcf().transFigure)
386
387 legend3 = "{:.3f}\n{:.3f}".format(float(singl)/sum(data_o), float(sum(data_o[data_o > 20])) / sum(data_o))
388 plt.text(0.84, 0.03, legend3, size=10, transform=plt.gcf().transFigure)
389
390 legend = "PE reads\nunique\n{:.3f}\n{:.3f}\n{:.3f}\n{:,}".format(
391 float(sum(dataAB_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)),
392 float(sum(dataBA_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)),
393 float(sum(duplTags_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)),
394 (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)))
395 plt.text(0.74, 0.09, legend, size=10, transform=plt.gcf().transFigure)
396
397 legend = "total\n{:.3f}\n{:.3f}\n{:.3f} ({:.3f})\n{:,}".format(
398 float(sum(dataAB_o)) / (sum(ab_o) + sum(ba_o)),
399 float(sum(dataBA_o)) / (sum(ab_o) + sum(ba_o)),
400 float(sum(duplTags_o)) / (sum(ab_o) + sum(ba_o)),
401 float(sum(duplTags_double_o)) / (sum(ab_o) + sum(ba_o)), (sum(ab_o) + sum(ba_o)))
402 plt.text(0.84, 0.09, legend, size=10, transform=plt.gcf().transFigure)
328 403
329 pdf.savefig(fig) 404 pdf.savefig(fig)
330 plt.close() 405 plt.close()
331 406
332 # write same information to a csv file 407 # write same information to a csv file
334 output_file.write("\nDataset:{}{}\n".format(sep, name_file)) 409 output_file.write("\nDataset:{}{}\n".format(sep, name_file))
335 output_file.write("max. family size:{}{}\n".format(sep, max(integers))) 410 output_file.write("max. family size:{}{}\n".format(sep, max(integers)))
336 output_file.write("absolute frequency:{}{}\n".format(sep, count[len(count) - 1])) 411 output_file.write("absolute frequency:{}{}\n".format(sep, count[len(count) - 1]))
337 output_file.write("relative frequency:{}{:.3f}\n\n".format(sep, float(count[len(count) - 1]) / sum(count))) 412 output_file.write("relative frequency:{}{:.3f}\n\n".format(sep, float(count[len(count) - 1]) / sum(count)))
338 413
339 output_file.write("{}singletons:{}{}family size > 20:\n".format(sep, sep, sep)) 414 output_file.write("{}singletons:{}{}{}family size > 20:\n".format(sep, sep, sep, sep))
340 output_file.write("{}absolute nr.{}rel. freq{}absolute nr.{}rel. freq{}total length\n".format(sep, sep, sep, sep, sep)) 415 output_file.write("{}nr. of tags{}rel. freq of tags{}rel.freq of PE reads{}nr. of tags{}rel. freq of tags{}nr. of PE reads{}rel. freq of PE reads{}total nr. of tags{}total nr. of PE reads\n".format(sep, sep, sep, sep, sep, sep, sep, sep, sep))
341 output_file.write("{}{}{}{}{:.3f}{}{}{}{:.3f}{}{}\n\n".format(name_file, sep, singl.astype(int), sep, singl / len(data), sep, last.astype(int), sep, last / len(data), sep, len(data))) 416 output_file.write("{}{}{}{}{:.3f}{}{:.3f}{}{}{}{:.3f}{}{}{}{:.3f}{}{}{}{}\n\n".format(
417 name_file, sep, singl.astype(int), sep, singl / len(data), sep, float(singl)/sum(data_o), sep,
418 last.astype(int), sep, last / len(data), sep, sum(data_o[data_o > 20]), sep, float(sum(data_o[data_o > 20])) / sum(data_o), sep, len(data), sep, sum(data_o)))
342 419
343 # information for FS >= 1 420 # information for FS >= 1
344 output_file.write("The unique frequencies were calculated from the dataset where the tags occured only once (=ab without DCS, ba without DCS)\nWhereas the total frequencies were calculated from the whole dataset (=including the DCS).\n\n") 421 output_file.write("The unique frequencies were calculated from the dataset where the tags occured only once (=ab without DCS, ba without DCS)\n"
345 output_file.write("FS >= 1{}{}unique:{}total:\n".format(sep, sep, sep)) 422 "Whereas the total frequencies were calculated from the whole dataset (=including the DCS).\n\n")
346 output_file.write("nr./rel. freq of ab={}{}{}{:.3f}{}{:.3f}\n".format(sep, len(dataAB), sep, float(len(dataAB)) / (len(dataAB) + len(dataBA) + len( duplTags)), sep, float(len(dataAB)) / (len(ab) + len(ba)))) 423 output_file.write("FS >= 1{}nr. of tags{}nr. of PE reads{}rel. freq of tags{}{}rel. freq of PE reads:\n".format(sep, sep, sep, sep, sep))
347 output_file.write("nr./rel. freq of ba={}{}{}{:.3f}{}{:.3f}\n".format(sep, len(dataBA), sep, float(len(dataBA)) / (len(dataBA) + len(dataBA) + len(duplTags)), sep, float(len(dataBA)) / (len(ba) + len(ba)))) 424 output_file.write("{}{}{}unique:{}total{}unique{}total:\n".format(sep, sep, sep, sep, sep, sep))
348 output_file.write("nr./rel. freq of DCS (total)={}{} ({}){}{:.3f}{}{:.3f} ({:.3f})\n".format(sep, len(duplTags), len(duplTags_double), sep, float(len(duplTags)) / (len(dataAB) + len(dataBA) + len(duplTags)), sep, float(len(duplTags)) / ( len(ab) + len(ba)), float(len(duplTags_double)) / (len(ab) + len(ba)))) 425 output_file.write("SSCS ab{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format(
349 output_file.write("length of dataset={}{}{}{}{}{}\n".format(sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, (len(ab) + len(ba)))) 426 sep, len(dataAB), sep, sum(dataAB_o), sep, float(len(dataAB)) / (len(dataAB) + len(dataBA) + len(duplTags)),
427 sep, float(sum(dataAB_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep,
428 float(len(dataAB)) / (len(ab) + len(ba)), sep, float(sum(dataAB_o)) / (sum(ab_o) + sum(ba_o))))
429 output_file.write("SSCS ba{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format(
430 sep, len(dataBA), sep, sum(dataBA_o), sep, float(len(dataBA)) / (len(dataBA) + len(dataBA) + len(duplTags)),
431 sep, float(sum(dataBA_o)) / (sum(dataBA_o) + sum(dataBA_o) + sum(duplTags_o)), sep, float(len(dataBA)) / (len(ba) + len(ba)),
432 sep, float(sum(dataBA_o)) / (sum(ba_o) + sum(ba_o))))
433 output_file.write("DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format(
434 sep, len(duplTags), len(duplTags_double), sep, sum(duplTags_o), sum(duplTags_double_o), sep,
435 float(len(duplTags)) / (len(dataAB) + len(dataBA) + len(duplTags)), sep,
436 float(len(duplTags)) / (len(ab) + len(ba)), float(len(duplTags_double)) / (len(ab) + len(ba)), sep,
437 float(sum(duplTags_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep,
438 float(sum(duplTags_o)) / (sum(ab_o) + sum(ba_o)), float(sum(duplTags_double_o)) / (sum(ab_o) + sum(ba_o))))
439 output_file.write("total nr. of tags{}{}{}{}{}{}{}{}{}{}{}{}\n".format(
440 sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep,
441 (len(dataAB) + len(dataBA) + len(duplTags)), sep, (len(ab) + len(ba)), sep,
442 (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, (sum(ab_o) + sum(ba_o))))
350 # information for FS >= 3 443 # information for FS >= 3
351 output_file.write("FS >= 3{}{}unique:{}total:\n".format(sep, sep, sep)) 444 output_file.write("\nFS >= 3{}nr. of tags{}nr. of PE reads{}rel. freq of tags{}{}rel. freq of PE reads:\n".format(sep, sep, sep, sep, sep))
352 output_file.write("nr./rel. freq of ab={}{}{}{:.3f}{}{:.3f}\n".format(sep, len(dataAB_FS3), sep, float(len(dataAB_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, float(len(dataAB_FS3)) / (len(ab_FS3) + len(ba_FS3)))) 445 output_file.write("{}{}{}unique:{}total{}unique{}total:\n".format(sep, sep, sep, sep, sep, sep))
353 output_file.write("nr./rel. freq of ba={}{}{}{:.3f}{}{:.3f}\n".format(sep, len(dataBA_FS3), sep, float(len(dataBA_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, float(len(dataBA_FS3)) / (len(ba_FS3) + len(ba_FS3)))) 446 output_file.write("SSCS ab{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format(
354 output_file.write("nr./rel. freq of DCS (total)={}{} ({}){}{:.3f}{}{:.3f} ({:.3f})\n".format(sep, len(duplTags_FS3), duplTags_double_FS3, sep, float(len( duplTags_FS3)) / (len(dataBA_FS3) + len(duplTags_FS3)), sep, float(len(duplTags_FS3)) / (len(ab_FS3) + len(ba_FS3)), float(duplTags_double_FS3) / (len(ab_FS3) + len(ba_FS3)))) 447 sep, len(dataAB_FS3), sep, sum(dataAB_FS3_o), sep,
355 output_file.write("length of dataset={}{}{}{}{}{}\n".format(sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (len(ab_FS3) + len(ba_FS3)))) 448 float(len(dataAB_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep,
449 float(len(dataAB_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + duplTags_double_FS3),
450 sep, float(sum(dataAB_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)),
451 sep, float(sum(dataAB_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o)))
452 output_file.write("SSCS ba{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format(
453 sep, len(dataBA_FS3), sep, sum(dataBA_FS3_o), sep,
454 float(len(dataBA_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + len(duplTags_FS3)),
455 sep, float(len(dataBA_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + duplTags_double_FS3),
456 sep, float(sum(dataBA_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)),
457 sep, float(sum(dataBA_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o)))
458 output_file.write("DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format(
459 sep, len(duplTags_FS3), duplTags_double_FS3, sep, sum(duplTags_FS3_o), duplTags_double_FS3_o, sep,
460 float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep,
461 float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3),
462 float(duplTags_double_FS3) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3),
463 sep, float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), sep,
464 float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o),
465 float(duplTags_double_FS3_o) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o)))
466 output_file.write("total nr. of tags{}{}{}{}{}{}{}{}{}{}{}{}\n".format(
467 sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)),
468 sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3),
469 sep, (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), sep, (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o)))
356 470
357 output_file.write("\nValues from family size distribution\n") 471 output_file.write("\nValues from family size distribution\n")
358 output_file.write("{}duplex{}ab{}ba{}sum\n".format(sep, sep, sep, sep)) 472 output_file.write("{}duplex{}ab{}ba{}sum\n".format(sep, sep, sep, sep))
359 for dx, ab, ba, fs in zip(counts[0][0], counts[0][1], counts[0][2], counts[1]): 473 for dx, ab, ba, fs in zip(counts[0][0], counts[0][1], counts[0][2], counts[1]):
360 if fs == 21: 474 if fs == 21: