Mercurial > repos > mheinzl > fsd
diff fsd.py @ 45:6651e76baca1 draft
planemo upload for repository https://github.com/monikaheinzl/duplexanalysis_galaxy/tree/master/tools/fsd commit 033dd7b750f68e8aa68f327d7d72bd311ddbee4e-dirty
| author | mheinzl |
|---|---|
| date | Tue, 27 Aug 2019 07:36:53 -0400 |
| parents | a76af7fd9fca |
| children | 901827154779 |
line wrap: on
line diff
--- a/fsd.py Wed Aug 14 13:03:14 2019 -0400 +++ b/fsd.py Tue Aug 27 07:36:53 2019 -0400 @@ -41,6 +41,7 @@ parser.add_argument('--inputFile4', default=None, help='Tabular File with three columns: ab or ba, tag and family size.') parser.add_argument('--inputName4') parser.add_argument('--log_axis', action="store_false", help='Transform y axis in log scale.') + parser.add_argument('--rel_freq', action="store_false", help='If False, the relative frequencies are displayed.') parser.add_argument('--output_pdf', default="data.pdf", type=str, help='Name of the pdf file.') parser.add_argument('--output_tabular', default="data.tabular", type=str, help='Name of the tabular file.') return parser @@ -61,6 +62,7 @@ fourthFile = args.inputFile4 name4 = args.inputName4 log_axis = args.log_axis + rel_freq = args.rel_freq title_file = args.output_tabular title_file2 = args.output_pdf @@ -78,7 +80,7 @@ data_array_list = [] list_to_plot_original = [] colors = [] - bins = numpy.arange(1, 22) + bins = numpy.arange(1, 22) with open(title_file, "w") as output_file, PdfPages(title_file2) as pdf: fig = plt.figure() fig.subplots_adjust(left=0.12, right=0.97, bottom=0.23, top=0.94, hspace=0) @@ -98,6 +100,8 @@ # data1[bigFamilies] = 22 data1 = numpy.clip(integers, bins[0], bins[-1]) name1 = name1.split(".tabular")[0] + if len(name1) > 40: + name1 = name1[:40] list_to_plot.append(data1) label.append(name1) data_array_list.append(file1) @@ -106,21 +110,24 @@ fig.text(0.05, 0.11, legend, size=10, transform=plt.gcf().transFigure) fig2.text(0.05, 0.11, legend, size=10, transform=plt.gcf().transFigure) - legend1 = "singletons:\nnr. of tags\n{:,} ({:.3f})".format(numpy.bincount(data1)[1], float(numpy.bincount(data1)[1]) / len(data1)) + legend1 = "singletons:\nnr. of tags\n{:,} ({:.3f})".format(numpy.bincount(data1)[1], + float(numpy.bincount(data1)[1]) / len(data1)) fig.text(0.32, 0.11, legend1, size=10, transform=plt.gcf().transFigure) fig2.text(0.32, 0.11, legend1, size=10, transform=plt.gcf().transFigure) - legend3b = "PE reads\n{:,} ({:.3f})".format(numpy.bincount(data1)[1], float(numpy.bincount(data1)[1]) / sum(integers)) + legend3b = "PE reads\n{:,} ({:.3f})".format(numpy.bincount(data1)[1], + float(numpy.bincount(data1)[1]) / sum(integers)) fig.text(0.45, 0.11, legend3b, size=10, transform=plt.gcf().transFigure) fig2.text(0.45, 0.11, legend3b, size=10, transform=plt.gcf().transFigure) legend4 = "family size > 20:\nnr. of tags\n{:,} ({:.3f})".format(len(integers[integers > 20]), - float(sum(integers[integers > 20])) - / sum(integers)) + float(len(integers[integers > 20])) + / len(integers)) fig.text(0.58, 0.11, legend4, size=10, transform=plt.gcf().transFigure) fig2.text(0.58, 0.11, legend4, size=10, transform=plt.gcf().transFigure) - legend5 = "PE reads\n{:,} ({:.3f})".format(sum(integers[integers > 20]), float(sum(integers[integers > 20])) / sum(integers)) + legend5 = "PE reads\n{:,} ({:.3f})".format(sum(integers[integers > 20]), + float(sum(integers[integers > 20])) / sum(integers)) fig.text(0.70, 0.11, legend5, size=10, transform=plt.gcf().transFigure) fig2.text(0.70, 0.11, legend5, size=10, transform=plt.gcf().transFigure) @@ -145,6 +152,8 @@ data2 = numpy.clip(integers2, bins[0], bins[-1]) list_to_plot.append(data2) name2 = name2.split(".tabular")[0] + if len(name2) > 40: + name2 = name2[:40] label.append(name2) data_array_list.append(file2) @@ -159,11 +168,13 @@ fig.text(0.45, 0.09, legend3, size=10, transform=plt.gcf().transFigure) fig2.text(0.45, 0.09, legend3, size=10, transform=plt.gcf().transFigure) - legend4 = "{:,} ({:.3f})".format(len(integers2[integers2 > 20]), float(sum(integers2[integers2 > 20])) / sum(integers2)) + legend4 = "{:,} ({:.3f})".format(len(integers2[integers2 > 20]), + float(len(integers2[integers2 > 20])) / len(integers2)) fig.text(0.58, 0.09, legend4, size=10, transform=plt.gcf().transFigure) fig2.text(0.58, 0.09, legend4, size=10, transform=plt.gcf().transFigure) - legend5 = "{:,} ({:.3f})".format(sum(integers2[integers2 > 20]), float(sum(integers2[integers2 > 20])) / sum(integers2)) + legend5 = "{:,} ({:.3f})".format(sum(integers2[integers2 > 20]), + float(sum(integers2[integers2 > 20])) / sum(integers2)) fig.text(0.70, 0.09, legend5, size=10, transform=plt.gcf().transFigure) fig2.text(0.70, 0.09, legend5, size=10, transform=plt.gcf().transFigure) @@ -188,6 +199,8 @@ data3 = numpy.clip(integers3, bins[0], bins[-1]) list_to_plot.append(data3) name3 = name3.split(".tabular")[0] + if len(name3) > 40: + name3 = name3[:40] label.append(name3) data_array_list.append(file3) @@ -198,11 +211,13 @@ fig.text(0.32, 0.07, legend1, size=10, transform=plt.gcf().transFigure) fig2.text(0.32, 0.07, legend1, size=10, transform=plt.gcf().transFigure) - legend3b = "{:,} ({:.3f})".format(numpy.bincount(data3)[1], float(numpy.bincount(data3)[1]) / sum(integers3)) + legend3b = "{:,} ({:.3f})".format(numpy.bincount(data3)[1], + float(numpy.bincount(data3)[1]) / sum(integers3)) fig.text(0.45, 0.07, legend3b, size=10, transform=plt.gcf().transFigure) fig2.text(0.45, 0.07, legend3b, size=10, transform=plt.gcf().transFigure) - legend4 = "{:,} ({:.3f})".format(len(integers3[integers3 > 20]), float(sum(integers3[integers3 > 20])) / sum(integers3)) + legend4 = "{:,} ({:.3f})".format(len(integers3[integers3 > 20]), + float(len(integers3[integers3 > 20])) / len(integers3)) fig.text(0.58, 0.07, legend4, size=10, transform=plt.gcf().transFigure) fig2.text(0.58, 0.07, legend4, size=10, transform=plt.gcf().transFigure) @@ -231,6 +246,8 @@ data4 = numpy.clip(integers4, bins[0], bins[-1]) list_to_plot.append(data4) name4 = name4.split(".tabular")[0] + if len(name4) > 40: + name4 = name4[:40] label.append(name4) data_array_list.append(file4) @@ -241,11 +258,13 @@ fig.text(0.32, 0.05, legend1, size=10, transform=plt.gcf().transFigure) fig2.text(0.32, 0.05, legend1, size=10, transform=plt.gcf().transFigure) - legend3b = "{:,} ({:.3f})".format(numpy.bincount(data4)[1], float(numpy.bincount(data4)[1]) / sum(integers4)) + legend3b = "{:,} ({:.3f})".format(numpy.bincount(data4)[1], + float(numpy.bincount(data4)[1]) / sum(integers4)) fig.text(0.45, 0.05, legend3b, size=10, transform=plt.gcf().transFigure) fig2.text(0.45, 0.05, legend3b, size=10, transform=plt.gcf().transFigure) - legend4 = "{:,} ({:.3f})".format(len(integers4[integers4 > 20]), float(sum(integers4[integers4 > 20])) / sum(integers4)) + legend4 = "{:,} ({:.3f})".format(len(integers4[integers4 > 20]), + float(len(integers4[integers4 > 20])) / len(integers4)) fig.text(0.58, 0.05, legend4, size=10, transform=plt.gcf().transFigure) fig2.text(0.58, 0.05, legend4, size=10, transform=plt.gcf().transFigure) @@ -265,137 +284,129 @@ maximumX = numpy.amax(numpy.concatenate(list_to_plot)) minimumX = numpy.amin(numpy.concatenate(list_to_plot)) list_to_plot2 = list_to_plot - to_plot = ["Absolute frequencies", "Relative frequencies"] - plt.xticks([], []) - plt.yticks([], []) - fig.suptitle('Family Size Distribution (tags)', fontsize=14) - for l in range(len(to_plot)): - ax = fig.add_subplot(2, 1, l+1) - ticks = numpy.arange(1, 22, 1) - ticks1 = map(str, ticks) - ticks1[len(ticks1) - 1] = ">20" + if rel_freq: + ylab = "Relative Frequency" + else: + ylab = "Absolute Frequency" - if to_plot[l] == "Relative frequencies": - w = [numpy.zeros_like(data) + 1. / len(data) for data in list_to_plot2] - counts_rel = ax.hist(list_to_plot2, weights=w, - bins=numpy.arange(1, 23), stacked=False, edgecolor="black", - linewidth=1, label=label, align="left", alpha=0.7, rwidth=0.8) - ax.set_ylim(0, 1.07) - else: - counts = ax.hist(list_to_plot2, bins=numpy.arange(1, 23), stacked=False, edgecolor="black", linewidth=1, label=label, align="left", alpha=0.7, rwidth=0.8, color=colors) - ax.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(0.9, 1)) - - ax.set_xticks(numpy.array(ticks)) - ax.set_xticklabels(ticks1) - - ax.set_ylabel(to_plot[l], fontsize=14) - ax.set_xlabel("Family size", fontsize=14) - if log_axis: - ax.set_yscale('log') - ax.grid(b=True, which="major", color="#424242", linestyle=":") - ax.margins(0.01, None) - pdf.savefig(fig) - plt.close() + # PLOT FSD based on tags + fig.suptitle('Family Size Distribution (FSD) based on families', fontsize=14) + ax = fig.add_subplot(1, 1, 1) + ticks = numpy.arange(1, 22, 1) + ticks1 = map(str, ticks) + ticks1[len(ticks1) - 1] = ">20" + ax.set_xticks([], []) + if rel_freq: + w = [numpy.zeros_like(data) + 1. / len(data) for data in list_to_plot2] + counts = ax.hist(list_to_plot2, weights=w, + bins=numpy.arange(1, 23), stacked=False, edgecolor="black", color=colors, + linewidth=1, label=label, align="left", alpha=0.7, rwidth=0.8) + ax.set_ylim(0, 1.07) + else: + counts = ax.hist(list_to_plot2, bins=numpy.arange(1, 23), stacked=False, edgecolor="black", linewidth=1, + label=label, align="left", alpha=0.7, rwidth=0.8, color=colors) + ax.set_xticks(numpy.array(ticks)) + ax.set_xticklabels(ticks1) + ax.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(0.9, 1)) - fig2.suptitle('Family Size Distribution (PE reads)', fontsize=14) - for l in range(len(to_plot)): - ax = fig2.add_subplot(2, 1, l + 1) - ticks = numpy.arange(minimumX, maximumX + 1) - ticks = numpy.arange(1, 22) - ticks1 = map(str, ticks) - ticks1[len(ticks1) - 1] = ">20" - reads = [] - reads_rel = [] - - barWidth = 0 - (len(list_to_plot)+1)/2 * 1./(len(list_to_plot)+1) - - for i in range(len(list_to_plot2)): - x = list(numpy.arange(1, 22).astype(float)) - unique, c = numpy.unique(list_to_plot2[i], return_counts=True) - y = unique * c - if sum(list_to_plot_original[i] > 20) > 0: - y[len(y) - 1] = sum(list_to_plot_original[i][list_to_plot_original[i] > 20]) - y = [y[x[idx] == unique][0] if x[idx] in unique else 0 for idx in range(len(x))] - reads.append(y) - reads_rel.append(list(numpy.float_(y)) / sum(y)) - #x = [xi + barWidth for xi in x] + ax.set_ylabel(ylab, fontsize=14) + ax.set_xlabel("Family size", fontsize=14) + if log_axis: + ax.set_yscale('log') + ax.grid(b=True, which="major", color="#424242", linestyle=":") + ax.margins(0.01, None) + pdf.savefig(fig) + #plt.close() - if len(list_to_plot2) == 1: - x = [xi * 0.5 for xi in x] - w = 0.4 - else: - x = [xi + barWidth for xi in x] - w = 1./(len(list_to_plot) + 1) - if to_plot[l] == "Relative frequencies": - counts2_rel = ax.bar(x, list(numpy.float_(y)) / numpy.sum(y), align="edge", width=w, - edgecolor="black", label=label[i],linewidth=1, alpha=0.7, color=colors[i]) - ax.set_ylim(0, 1.07) - else: - #y = list(y.reshape((len(y)))) - counts2 = ax.bar(x, y, align="edge", width=w, edgecolor="black", label=label[i], linewidth=1, - alpha=0.7, color=colors[i]) - if i == len(list_to_plot2)-1: - barWidth += 1. / (len(list_to_plot) + 1) + 1. / (len(list_to_plot) + 1) - else: - barWidth += 1. / (len(list_to_plot) + 1) - - if to_plot[l] == "Absolute frequencies": - ax.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(0.9, 1)) - else: - ax.set_xlabel("Family size", fontsize=14) + # PLOT FSD based on PE reads + fig2.suptitle('Family Size Distribution (FSD) based on PE reads', fontsize=14) + ax2 = fig2.add_subplot(1, 1, 1) + ticks = numpy.arange(1, 22) + ticks1 = map(str, ticks) + ticks1[len(ticks1) - 1] = ">20" + reads = [] + reads_rel = [] + + barWidth = 0 - (len(list_to_plot) + 1) / 2 * 1. / (len(list_to_plot) + 1) + ax2.set_xticks([], []) + + for i in range(len(list_to_plot2)): + x = list(numpy.arange(1, 22).astype(float)) + unique, c = numpy.unique(list_to_plot2[i], return_counts=True) + y = unique * c + if sum(list_to_plot_original[i] > 20) > 0: + y[len(y) - 1] = sum(list_to_plot_original[i][list_to_plot_original[i] > 20]) + y = [y[x[idx] == unique][0] if x[idx] in unique else 0 for idx in range(len(x))] + reads.append(y) + reads_rel.append(list(numpy.float_(y)) / sum(y)) if len(list_to_plot2) == 1: - ax.set_xticks(numpy.array([xi + 0.2 for xi in x])) + x = [xi * 0.5 for xi in x] + w = 0.4 + else: + x = [xi + barWidth for xi in x] + w = 1. / (len(list_to_plot) + 1) + if rel_freq: + counts2_rel = ax2.bar(x, list(numpy.float_(y)) / numpy.sum(y), align="edge", width=w, + edgecolor="black", label=label[i], linewidth=1, alpha=0.7, color=colors[i]) + ax2.set_ylim(0, 1.07) + else: + counts2 = ax2.bar(x, y, align="edge", width=w, edgecolor="black", label=label[i], linewidth=1, + alpha=0.7, color=colors[i]) + + if i == len(list_to_plot2) - 1: + barWidth += 1. / (len(list_to_plot) + 1) + 1. / (len(list_to_plot) + 1) else: - ax.set_xticks(numpy.array(ticks)) - ax.set_xticklabels(ticks1) - ax.set_ylabel(to_plot[l], fontsize=14) - if log_axis: - ax.set_yscale('log') - ax.grid(b=True, which="major", color="#424242", linestyle=":") - ax.margins(0.01, None) + barWidth += 1. / (len(list_to_plot) + 1) + + ax2.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(0.9, 1)) + + if len(list_to_plot2) == 1: + ax2.set_xticks(numpy.array([xi + 0.2 for xi in x])) + else: + ax2.set_xticks(numpy.array(ticks)) + ax2.set_xticklabels(ticks1) + ax2.set_xlabel("Family size", fontsize=14) + ax2.set_ylabel(ylab, fontsize=14) + if log_axis: + ax2.set_yscale('log') + ax2.grid(b=True, which="major", color="#424242", linestyle=":") + ax2.margins(0.01, None) pdf.savefig(fig2) plt.close() # write data to CSV file tags - output_file.write("Values from family size distribution with all datasets (tags)\n") + counts = [numpy.bincount(d, minlength=22)[1:] for d in list_to_plot2] # original counts of family sizes + output_file.write("Values from family size distribution with all datasets based on families\n") output_file.write("\nFamily size") for i in label: output_file.write("{}{}".format(sep, i)) - # output_file.write("{}sum".format(sep)) output_file.write("\n") j = 0 - for fs in counts[1][0:len(counts[1]) - 1]: + for fs in bins: if fs == 21: fs = ">20" else: fs = "={}".format(fs) output_file.write("FS{}{}".format(fs, sep)) - if len(label) == 1: - output_file.write("{}{}".format(int(counts[0][j]), sep)) - else: - for n in range(len(label)): - output_file.write("{}{}".format(int(counts[0][n][j]), sep)) + for n in range(len(label)): + output_file.write("{}{}".format(int(counts[n][j]), sep)) output_file.write("\n") j += 1 output_file.write("sum{}".format(sep)) - if len(label) == 1: - output_file.write("{}{}".format(int(sum(counts[0])), sep)) - else: - for i in counts[0]: - output_file.write("{}{}".format(int(sum(i)), sep)) + for i in counts: + output_file.write("{}{}".format(int(sum(i)), sep)) # write data to CSV file PE reads - output_file.write("\n\nValues from family size distribution with all datasets (PE reads)\n") + output_file.write("\n\nValues from family size distribution with all datasets based on PE reads\n") output_file.write("\nFamily size") for i in label: output_file.write("{}{}".format(sep, i)) - # output_file.write("{}sum".format(sep)) output_file.write("\n") j = 0 - + for fs in bins: if fs == 21: fs = ">20" @@ -410,7 +421,7 @@ output_file.write("\n") j += 1 output_file.write("sum{}".format(sep)) - if len(label) == 1: + if len(label) == 1: output_file.write("{}{}".format(int(sum(numpy.concatenate(reads))), sep)) else: for i in reads: @@ -461,6 +472,7 @@ dataBA_o = ba_o[numpy.in1d(baSeq, d, invert=True)] list1 = [duplTags_double, dataAB, dataBA] # list for plotting + list1_o = [duplTags_double_o, dataAB_o, dataBA_o] # list for plotting # information for family size >= 3 dataAB_FS3 = dataAB[dataAB >= 3] @@ -483,20 +495,30 @@ fig = plt.figure() plt.subplots_adjust(left=0.12, right=0.97, bottom=0.3, top=0.94, hspace=0) - counts = plt.hist(list1, bins=numpy.arange(1, 23), stacked=True, label=["duplex", "ab", "ba"], - edgecolor="black", linewidth=1, align="left", color=["#FF0000", "#5FB404", "#FFBF00"], - rwidth=0.8) + + if rel_freq: + w = [numpy.zeros_like(d) + 1. / len(numpy.concatenate(list1)) for d in list1] + counts = plt.hist(list1, bins=numpy.arange(1, 23), stacked=True, label=["duplex", "ab", "ba"], weights=w, + edgecolor="black", linewidth=1, align="left", color=["#FF0000", "#5FB404", "#FFBF00"], + rwidth=0.8) + plt.ylim(0, 1.07) + else: + counts = plt.hist(list1, bins=numpy.arange(1, 23), stacked=True, label=["duplex", "ab", "ba"], + edgecolor="black", linewidth=1, align="left", color=["#FF0000", "#5FB404", "#FFBF00"], + rwidth=0.8) + # tick labels of x axis ticks = numpy.arange(1, 22, 1) ticks1 = map(str, ticks) ticks1[len(ticks1) - 1] = ">20" plt.xticks(numpy.array(ticks), ticks1) - singl = counts[0][2][0] # singletons - last = counts[0][2][len(counts[0][0]) - 1] # large families + # singl = counts[0][2][0] # singletons + singl = len(data_o[data_o == 1]) + last = len(data_o[data_o > 20]) # large families if log_axis: plt.yscale('log') plt.legend(loc='upper right', fontsize=14, bbox_to_anchor=(0.9, 1), frameon=True) - plt.title(name_file, fontsize=14) + plt.title("{}: FSD based on families".format(name_file), fontsize=14) plt.xlabel("Family size", fontsize=14) plt.ylabel("Absolute Frequency", fontsize=14) plt.margins(0.01, None) @@ -506,31 +528,47 @@ legend = "SSCS ab= \nSSCS ba= \nDCS (total)= \ntotal nr. of tags=" plt.text(0.1, 0.09, legend, size=10, transform=plt.gcf().transFigure) - 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)), (len(ab) + len(ba))) + 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)), + (len(ab) + len(ba))) plt.text(0.23, 0.09, legend, size=10, transform=plt.gcf().transFigure) - 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)), (sum(ab_o) + sum(ba_o))) + 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)), + (sum(ab_o) + sum(ba_o))) plt.text(0.38, 0.09, legend5, size=10, transform=plt.gcf().transFigure) - 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))) + 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))) plt.text(0.54, 0.09, legend, size=10, transform=plt.gcf().transFigure) - 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))) + 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))) plt.text(0.64, 0.09, legend, size=10, transform=plt.gcf().transFigure) legend1 = "\nsingletons:\nfamily size > 20:" plt.text(0.1, 0.03, legend1, size=10, transform=plt.gcf().transFigure) - legend4 = "{:,}\n{:,}".format(singl.astype(int), last.astype(int)) + legend4 = "{:,}\n{:,}".format(singl, last) plt.text(0.23, 0.03, legend4, size=10, transform=plt.gcf().transFigure) - - legend3 = "{:.3f}\n{:.3f}".format(singl / len(data), last / len(data)) + legend3 = "{:.3f}\n{:.3f}".format(float(singl) / len(data), float(last) / len(data)) plt.text(0.64, 0.03, legend3, size=10, transform=plt.gcf().transFigure) legend3 = "\n\n{:,}".format(sum(data_o[data_o > 20])) plt.text(0.38, 0.03, legend3, size=10, transform=plt.gcf().transFigure) - legend3 = "{:.3f}\n{:.3f}".format(float(singl)/sum(data_o), float(sum(data_o[data_o > 20])) / sum(data_o)) + legend3 = "{:.3f}\n{:.3f}".format(float(singl) / sum(data_o), float(sum(data_o[data_o > 20])) / sum(data_o)) plt.text(0.84, 0.03, legend3, size=10, transform=plt.gcf().transFigure) legend = "PE reads\nunique\n{:.3f}\n{:.3f}\n{:.3f}\n{:,}".format( @@ -550,8 +588,137 @@ pdf.savefig(fig) plt.close() + # PLOT FSD based on PE reads + fig3 = plt.figure() + plt.subplots_adjust(left=0.12, right=0.97, bottom=0.3, top=0.94, hspace=0) + + fig3.suptitle("{}: FSD based on PE reads".format(name_file), fontsize=14) + ax2 = fig3.add_subplot(1, 1, 1) + ticks = numpy.arange(1, 22) + ticks1 = map(str, ticks) + ticks1[len(ticks1) - 1] = ">20" + reads = [] + reads_rel = [] + + #barWidth = 0 - (len(list_to_plot) + 1) / 2 * 1. / (len(list_to_plot) + 1) + ax2.set_xticks([], []) + + list_y = [] + label = ["duplex", "ab", "ba"] + col = ["#FF0000", "#5FB404", "#FFBF00"] + for i in range(len(list1)): + x = list(numpy.arange(1, 22).astype(float)) + unique, c = numpy.unique(list1[i], return_counts=True) + y = unique * c + if sum(list1_o[i] > 20) > 0: + y[len(y) - 1] = sum(list1_o[i][list1_o[i] > 20]) + y = [y[x[idx] == unique][0] if x[idx] in unique else 0 for idx in range(len(x))] + reads.append(y) + reads_rel.append(list(numpy.float_(y)) / sum(numpy.concatenate(list1_o))) + + if rel_freq: + y = list(numpy.float_(y)) / sum(numpy.concatenate(list1_o)) + ax2.set_ylim(0, 1.07) + else: + y = y + + list_y.append(y) + if i == 0: + counts2 = ax2.bar(x, y, align="edge", width=0.8, + edgecolor="black", label=label[0], + linewidth=1, alpha=1, color=col[0]) + elif i == 1: + counts2 = ax2.bar(x, y, bottom=list_y[i-1], align="edge", width=0.8, + edgecolor="black", label=label[1], + linewidth=1, alpha=1, color=col[1]) + elif i == 2: + bars = numpy.add(list_y[0], list_y[1]).tolist() + + counts2 = ax2.bar(x, y, bottom=bars, align="edge", width=0.8, + edgecolor="black", label=label[2], + linewidth=1, alpha=1, color=col[2]) + + ax2.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(0.9, 1)) + + singl = len(data_o[data_o == 1]) + last = len(data_o[data_o > 20]) # large families + + ax2.set_xticks(numpy.array(ticks)) + ax2.set_xticklabels(ticks1) + ax2.set_xlabel("Family size", fontsize=14) + ax2.set_ylabel(ylab, fontsize=14) + if log_axis: + ax2.set_yscale('log') + ax2.grid(b=True, which="major", color="#424242", linestyle=":") + ax2.margins(0.01, None) + + # extra information beneath the plot + legend = "SSCS ab= \nSSCS ba= \nDCS (total)= \ntotal nr. of tags=" + plt.text(0.1, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + 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)), + (len(ab) + len(ba))) + plt.text(0.23, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + 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)), + (sum(ab_o) + sum(ba_o))) + plt.text(0.38, 0.09, legend5, size=10, transform=plt.gcf().transFigure) + + 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))) + plt.text(0.54, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + 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))) + plt.text(0.64, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + legend1 = "\nsingletons:\nfamily size > 20:" + plt.text(0.1, 0.03, legend1, size=10, transform=plt.gcf().transFigure) + + legend4 = "{:,}\n{:,}".format(singl, last) + plt.text(0.23, 0.03, legend4, size=10, transform=plt.gcf().transFigure) + legend3 = "{:.3f}\n{:.3f}".format(float(singl) / len(data), float(last) / len(data)) + plt.text(0.64, 0.03, legend3, size=10, transform=plt.gcf().transFigure) + + legend3 = "\n\n{:,}".format(sum(data_o[data_o > 20])) + plt.text(0.38, 0.03, legend3, size=10, transform=plt.gcf().transFigure) + + legend3 = "{:.3f}\n{:.3f}".format(float(singl) / sum(data_o), float(sum(data_o[data_o > 20])) / sum(data_o)) + plt.text(0.84, 0.03, legend3, size=10, transform=plt.gcf().transFigure) + + legend = "PE reads\nunique\n{:.3f}\n{:.3f}\n{:.3f}\n{:,}".format( + float(sum(dataAB_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), + float(sum(dataBA_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), + float(sum(duplTags_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), + (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o))) + plt.text(0.74, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + legend = "total\n{:.3f}\n{:.3f}\n{:.3f} ({:.3f})\n{:,}".format( + float(sum(dataAB_o)) / (sum(ab_o) + sum(ba_o)), + float(sum(dataBA_o)) / (sum(ab_o) + sum(ba_o)), + float(sum(duplTags_o)) / (sum(ab_o) + sum(ba_o)), + float(sum(duplTags_double_o)) / (sum(ab_o) + sum(ba_o)), (sum(ab_o) + sum(ba_o))) + plt.text(0.84, 0.09, legend, size=10, transform=plt.gcf().transFigure) + + pdf.savefig(fig3) + plt.close() + # write same information to a csv file count = numpy.bincount(data_o) # original counts of family sizes + output_file.write("\nDataset:{}{}\n".format(sep, name_file)) output_file.write("max. family size:{}{}\n".format(sep, max(data_o))) output_file.write("absolute frequency:{}{}\n".format(sep, count[len(count) - 1])) @@ -560,37 +727,58 @@ output_file.write("median family size:{}{}\n".format(sep, numpy.median(numpy.array(data_o)))) output_file.write("mean family size:{}{}\n\n".format(sep, numpy.mean(numpy.array(data_o)))) - output_file.write("{}singletons:{}{}{}family size > 20:{}{}{}{}length of dataset:\n".format(sep, sep, sep, sep, sep, sep, sep, sep)) - 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)) + output_file.write( + "{}singletons:{}{}{}family size > 20:{}{}{}{}length of dataset:\n".format(sep, sep, sep, sep, sep, sep, + sep, sep)) + 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)) output_file.write("{}{}{}{}{:.3f}{}{:.3f}{}{}{}{:.3f}{}{}{}{:.3f}{}{}{}{}\n\n".format( - name_file, sep, singl.astype(int), sep, singl / len(data), sep, float(singl)/sum(data_o), sep, - 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))) + name_file, sep, singl, sep, float(singl) / len(data), sep, float(singl) / sum(data_o), sep, + last, sep, float(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))) # information for FS >= 1 - output_file.write("The unique frequencies were calculated from the dataset where the tags occured only once (=ab without DCS, ba without DCS)\n" - "Whereas the total frequencies were calculated from the whole dataset (=including the DCS).\n\n") - 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)) + output_file.write( + "The unique frequencies were calculated from the dataset where the tags occured only once (=ab without DCS, ba without DCS)\n" + "Whereas the total frequencies were calculated from the whole dataset (=including the DCS).\n\n") + 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)) output_file.write("{}{}{}unique:{}total{}unique{}total:\n".format(sep, sep, sep, sep, sep, sep)) output_file.write("SSCS ab{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format( - sep, len(dataAB), sep, sum(dataAB_o), sep, float(len(dataAB)) / (len(dataAB) + len(dataBA) + len(duplTags)), + sep, len(dataAB), sep, sum(dataAB_o), sep, + float(len(dataAB)) / (len(dataAB) + len(dataBA) + len(duplTags)), sep, float(sum(dataAB_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, float(len(dataAB)) / (len(ab) + len(ba)), sep, float(sum(dataAB_o)) / (sum(ab_o) + sum(ba_o)))) output_file.write("SSCS ba{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format( - sep, len(dataBA), sep, sum(dataBA_o), sep, float(len(dataBA)) / (len(dataBA) + len(dataBA) + len(duplTags)), - sep, float(sum(dataBA_o)) / (sum(dataBA_o) + sum(dataBA_o) + sum(duplTags_o)), sep, float(len(dataBA)) / (len(ba) + len(ba)), + sep, len(dataBA), sep, sum(dataBA_o), sep, + float(len(dataBA)) / (len(dataBA) + len(dataBA) + len(duplTags)), + sep, float(sum(dataBA_o)) / (sum(dataBA_o) + sum(dataBA_o) + sum(duplTags_o)), sep, + float(len(dataBA)) / (len(ba) + len(ba)), sep, float(sum(dataBA_o)) / (sum(ba_o) + sum(ba_o)))) - output_file.write("DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format( - sep, len(duplTags), len(duplTags_double), sep, sum(duplTags_o), sum(duplTags_double_o), 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)), sep, - float(sum(duplTags_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, - float(sum(duplTags_o)) / (sum(ab_o) + sum(ba_o)), float(sum(duplTags_double_o)) / (sum(ab_o) + sum(ba_o)))) + output_file.write( + "DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format( + sep, len(duplTags), len(duplTags_double), sep, sum(duplTags_o), sum(duplTags_double_o), 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)), sep, + float(sum(duplTags_o)) / (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, + float(sum(duplTags_o)) / (sum(ab_o) + sum(ba_o)), + float(sum(duplTags_double_o)) / (sum(ab_o) + sum(ba_o)))) output_file.write("total nr. of tags{}{}{}{}{}{}{}{}{}{}{}{}\n".format( - sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, + sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, + (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, (len(dataAB) + len(dataBA) + len(duplTags)), sep, (len(ab) + len(ba)), sep, (sum(dataAB_o) + sum(dataBA_o) + sum(duplTags_o)), sep, (sum(ab_o) + sum(ba_o)))) # information for FS >= 3 - 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)) + 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)) output_file.write("{}{}{}unique:{}total{}unique{}total:\n".format(sep, sep, sep, sep, sep, sep)) output_file.write("SSCS ab{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format( sep, len(dataAB_FS3), sep, sum(dataAB_FS3_o), sep, @@ -604,29 +792,61 @@ sep, float(len(dataBA_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + duplTags_double_FS3), sep, float(sum(dataBA_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), sep, float(sum(dataBA_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o))) - output_file.write("DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format( - sep, len(duplTags_FS3), duplTags_double_FS3, sep, sum(duplTags_FS3_o), duplTags_double_FS3_o, sep, - float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, - float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), - float(duplTags_double_FS3) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), - sep, float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), sep, - float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o), - float(duplTags_double_FS3_o) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o))) + output_file.write( + "DCS (total){}{} ({}){}{} ({}){}{:.3f}{}{:.3f} ({:.3f}){}{:.3f}{}{:.3f} ({:.3f})\n".format( + sep, len(duplTags_FS3), duplTags_double_FS3, sep, sum(duplTags_FS3_o), duplTags_double_FS3_o, sep, + float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, + float(len(duplTags_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), + float(duplTags_double_FS3) / (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), + sep, float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), + sep, + float(sum(duplTags_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o), + float(duplTags_double_FS3_o) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o))) output_file.write("total nr. of tags{}{}{}{}{}{}{}{}{}{}{}{}\n".format( - sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), - sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), - 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))) + sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, + (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), + sep, (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, + (len(dataAB_FS3) + len(dataBA_FS3) + duplTags_double_FS3), + 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))) - output_file.write("\nValues from family size distribution\n") + counts = [numpy.bincount(d, minlength=22)[1:] for d in list1] # original counts of family sizes + output_file.write("\nValues from family size distribution based on families\n") output_file.write("{}duplex{}ab{}ba{}sum\n".format(sep, sep, sep, sep)) - for dx, ab, ba, fs in zip(counts[0][0], counts[0][1], counts[0][2], counts[1]): + + j = 0 + for fs in bins: if fs == 21: fs = ">20" else: fs = "={}".format(fs) - ab1 = ab - dx - ba1 = ba - ab - output_file.write("FS{}{}{}{}{}{}{}{}{}\n".format(fs, sep, int(dx), sep, int(ab1), sep, int(ba1), sep, int(ba))) + output_file.write("FS{}{}".format(fs, sep)) + for n in range(3): + output_file.write("{}{}".format(int(counts[n][j]), sep)) + output_file.write("{}\n".format(counts[0][j] + counts[1][j] + counts[2][j])) + j += 1 + output_file.write("sum{}".format(sep)) + for i in counts: + output_file.write("{}{}".format(int(sum(i)), sep)) + output_file.write("{}\n".format(sum(counts[0] + counts[1] + counts[2]))) + + output_file.write("\nValues from family size distribution based on PE reads\n") + output_file.write("{}duplex{}ab{}ba{}sum\n".format(sep, sep, sep, sep)) + j = 0 + for fs in bins: + if fs == 21: + fs = ">20" + else: + fs = "={}".format(fs) + output_file.write("FS{}{}".format(fs, sep)) + for n in range(3): + output_file.write("{}{}".format(int(reads[n][j]), sep)) + output_file.write("{}\n".format(reads[0][j] + reads[1][j] + reads[2][j])) + j += 1 + output_file.write("sum{}".format(sep)) + for i in reads: + output_file.write("{}{}".format(int(sum(i)), sep)) + output_file.write("{}\n".format(sum(reads[0] + reads[1] + reads[2]))) print("Files successfully created!")