Mercurial > repos > mheinzl > fsd
view fsd.py @ 33:55c7e49eea88 draft
planemo upload for repository https://github.com/monikaheinzl/duplexanalysis_galaxy/tree/master/tools/fsd commit b8a2f7b7615b2bcd3b602027af31f4e677da94f6-dirty
| author | mheinzl |
|---|---|
| date | Mon, 03 Jun 2019 07:02:48 -0400 |
| parents | 5e1450eda945 |
| children | faa56fbf6b32 |
line wrap: on
line source
#!/usr/bin/env python # Family size distribution of SSCSs # # Author: Monika Heinzl, Johannes-Kepler University Linz (Austria) # Contact: monika.heinzl@edumail.at # # Takes at least one TABULAR file with tags before the alignment to the SSCS, but up to 4 files can be provided, as input. # The program produces a plot which shows the distribution of family sizes of the all SSCSs from the input files and # a tabular file with the data of the plot, as well as a TXT file with all tags of the DCS and their family sizes. # If only one file is provided, then a family size distribution, which is separated after SSCSs without a partner and DCSs, is produced. # Whereas a family size distribution with multiple data in one plot is produced, when more than one file (up to 4) is given. # USAGE: python FSD_Galaxy_1.4_commandLine_FINAL.py --inputFile1 filename --inputName1 filename --inputFile2 filename2 --inputName2 filename2 --inputFile3 filename3 --inputName3 filename3 --inputFile4 filename4 --inputName4 filename4 --log_axis --output_tabular outptufile_name_tabular --output_pdf outptufile_name_pdf import argparse import sys import os import matplotlib.pyplot as plt import numpy from matplotlib.backends.backend_pdf import PdfPages plt.switch_backend('agg') def readFileReferenceFree(file): with open(file, 'r') as dest_f: data_array = numpy.genfromtxt(dest_f, skip_header=0, delimiter='\t', comments='#', dtype='string') return(data_array) def make_argparser(): parser = argparse.ArgumentParser(description='Family Size Distribution of duplex sequencing data') parser.add_argument('--inputFile1', help='Tabular File with three columns: ab or ba, tag and family size.') parser.add_argument('--inputName1') parser.add_argument('--inputFile2', default=None, help='Tabular File with three columns: ab or ba, tag and family size.') parser.add_argument('--inputName2') parser.add_argument('--inputFile3', default=None, help='Tabular File with three columns: ab or ba, tag and family size.') parser.add_argument('--inputName3') 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('--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 def compare_read_families(argv): parser = make_argparser() args = parser.parse_args(argv[1:]) firstFile = args.inputFile1 name1 = args.inputName1 secondFile = args.inputFile2 name2 = args.inputName2 thirdFile = args.inputFile3 name3 = args.inputName3 fourthFile = args.inputFile4 name4 = args.inputName4 log_axis = args.log_axis title_file = args.output_tabular title_file2 = args.output_pdf sep = "\t" plt.rc('figure', figsize=(11.69, 8.27)) # A4 format plt.rcParams['patch.edgecolor'] = "black" plt.rcParams['axes.facecolor'] = "E0E0E0" # grey background color plt.rcParams['xtick.labelsize'] = 14 plt.rcParams['ytick.labelsize'] = 14 list_to_plot = [] label = [] data_array_list = [] list_to_plot_original = [] colors = [] 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) fig2 = plt.figure() fig2.subplots_adjust(left=0.12, right=0.97, bottom=0.23, top=0.94, hspace=0) # plt.subplots_adjust(bottom=0.25) if firstFile != str(None): file1 = readFileReferenceFree(firstFile) integers = numpy.array(file1[:, 0]).astype(int) # keep original family sizes list_to_plot_original.append(integers) colors.append("#0000FF") # for plot: replace all big family sizes by 22 # data1 = numpy.array(file1[:, 0]).astype(int) # bigFamilies = numpy.where(data1 > 20)[0] # data1[bigFamilies] = 22 data1 = numpy.clip(integers, bins[0], bins[-1]) name1 = name1.split(".tabular")[0] list_to_plot.append(data1) label.append(name1) data_array_list.append(file1) legend = "\n\n\n{}".format(name1) 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)) 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)) 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(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1].astype(int), float(numpy.bincount(data1)[len(numpy.bincount(data1)) - 1]) / len(data1)) 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)) 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) legend6 = "total nr. of\ntags\n{:,}".format(len(data1)) fig.text(0.82, 0.11, legend6, size=10, transform=plt.gcf().transFigure) fig2.text(0.82, 0.11, legend6, size=10, transform=plt.gcf().transFigure) legend6b = "PE reads\n{:,}".format(sum(integers)) fig.text(0.89, 0.11, legend6b, size=10, transform=plt.gcf().transFigure) fig2.text(0.89, 0.11, legend6b, size=10, transform=plt.gcf().transFigure) if secondFile != str(None): file2 = readFileReferenceFree(secondFile) integers2 = numpy.array(file2[:, 0]).astype(int) # keep original family sizes list_to_plot_original.append(integers2) colors.append("#298A08") # data2 = numpy.asarray(file2[:, 0]).astype(int) # bigFamilies2 = numpy.where(data2 > 20)[0] # data2[bigFamilies2] = 22 data2 = numpy.clip(integers2, bins[0], bins[-1]) list_to_plot.append(data2) name2 = name2.split(".tabular")[0] label.append(name2) data_array_list.append(file2) fig.text(0.05, 0.09, name2, size=10, transform=plt.gcf().transFigure) fig2.text(0.05, 0.09, name2, size=10, transform=plt.gcf().transFigure) legend1 = "{:,} ({:.3f})".format(numpy.bincount(data2)[1], float(numpy.bincount(data2)[1]) / len(data2)) fig.text(0.32, 0.09, legend1, size=10, transform=plt.gcf().transFigure) fig2.text(0.32, 0.09, legend1, size=10, transform=plt.gcf().transFigure) legend3 = "{:,} ({:.3f})".format(numpy.bincount(data2)[1], float(numpy.bincount(data2)[1]) / sum(integers2)) 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( numpy.bincount(data2)[len(numpy.bincount(data2)) - 1].astype(int), float(numpy.bincount(data2)[len(numpy.bincount(data2)) - 1]) / len(data2)) 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)) 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) legend6 = "{:,}".format(len(data2)) fig.text(0.82, 0.09, legend6, size=10, transform=plt.gcf().transFigure) fig2.text(0.82, 0.09, legend6, size=10, transform=plt.gcf().transFigure) legend6b = "{:,}".format(sum(integers2)) fig.text(0.89, 0.09, legend6b, size=10, transform=plt.gcf().transFigure) fig2.text(0.89, 0.09, legend6b, size=10, transform=plt.gcf().transFigure) if thirdFile != str(None): file3 = readFileReferenceFree(thirdFile) integers3 = numpy.array(file3[:, 0]).astype(int) # keep original family sizes list_to_plot_original.append(integers3) colors.append("#DF0101") # data3 = numpy.asarray(file3[:, 0]).astype(int) # bigFamilies3 = numpy.where(data3 > 20)[0] # data3[bigFamilies3] = 22 data3 = numpy.clip(integers3, bins[0], bins[-1]) list_to_plot.append(data3) name3 = name3.split(".tabular")[0] label.append(name3) data_array_list.append(file3) fig.text(0.05, 0.07, name3, size=10, transform=plt.gcf().transFigure) fig2.text(0.05, 0.07, name3, size=10, transform=plt.gcf().transFigure) legend1 = "{:,} ({:.3f})".format(numpy.bincount(data3)[1], float(numpy.bincount(data3)[1]) / len(data3)) 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)) 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( numpy.bincount(data3)[len(numpy.bincount(data3)) - 1].astype(int), float(numpy.bincount(data3)[len(numpy.bincount(data3)) - 1]) / len(data3)) 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) legend5 = "{:,} ({:.3f})".format(sum(integers3[integers3 > 20]), float(sum(integers3[integers3 > 20])) / sum(integers3)) fig.text(0.70, 0.07, legend5, size=10, transform=plt.gcf().transFigure) fig2.text(0.70, 0.07, legend5, size=10, transform=plt.gcf().transFigure) legend6 = "{:,}".format(len(data3)) fig.text(0.82, 0.07, legend6, size=10, transform=plt.gcf().transFigure) fig2.text(0.82, 0.07, legend6, size=10, transform=plt.gcf().transFigure) legend6b = "{:,}".format(sum(integers3)) fig.text(0.89, 0.07, legend6b, size=10, transform=plt.gcf().transFigure) fig2.text(0.89, 0.07, legend6b, size=10, transform=plt.gcf().transFigure) if fourthFile != str(None): file4 = readFileReferenceFree(fourthFile) integers4 = numpy.array(file4[:, 0]).astype(int) # keep original family sizes list_to_plot_original.append(integers4) colors.append("#04cec7") # data4 = numpy.asarray(file4[:, 0]).astype(int) # bigFamilies4 = numpy.where(data4 > 20)[0] # data4[bigFamilies4] = 22 data4 = numpy.clip(integers4, bins[0], bins[-1]) list_to_plot.append(data4) name4 = name4.split(".tabular")[0] label.append(name4) data_array_list.append(file4) fig.text(0.05, 0.05, name4, size=10, transform=plt.gcf().transFigure) fig2.text(0.05, 0.05, name4, size=10, transform=plt.gcf().transFigure) legend1 = "{:,} ({:.3f})".format(numpy.bincount(data4)[1], float(numpy.bincount(data4)[1]) / len(data4)) 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)) 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( numpy.bincount(data4)[len(numpy.bincount(data4)) - 1].astype(int), float(numpy.bincount(data4)[len(numpy.bincount(data4)) - 1]) / len(data4)) 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) legend5 = "{:,} ({:.3f})".format(sum(integers4[integers4 > 20]), float(sum(integers4[integers4 > 20])) / sum(integers4)) fig.text(0.70, 0.05, legend5, size=10, transform=plt.gcf().transFigure) fig2.text(0.70, 0.05, legend5, size=10, transform=plt.gcf().transFigure) legend6 = "{:,}".format(len(data4)) fig.text(0.82, 0.05, legend6, size=10, transform=plt.gcf().transFigure) fig2.text(0.82, 0.05, legend6, size=10, transform=plt.gcf().transFigure) legend6b = "{:,}".format(sum(integers4)) fig.text(0.89, 0.05, legend6b, size=10, transform=plt.gcf().transFigure) fig2.text(0.89, 0.05, legend6b, size=10, transform=plt.gcf().transFigure) 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) if maximumX > 20: ticks1[len(ticks1) - 1] = ">20" if to_plot[l] == "Relative frequencies": counts_rel = ax.hist(list_to_plot2, bins=numpy.arange(minimumX, maximumX + 2), stacked=False, edgecolor="black", linewidth=1, label=label, align="left", alpha=0.8, rwidth=0.8, normed=True) else: counts = ax.hist(list_to_plot2, bins=numpy.arange(minimumX, maximumX + 2), stacked=False, edgecolor="black", linewidth=1, label=label, align="left", alpha=0.8, rwidth=0.8) 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() 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) ticks1 = map(str, ticks) if maximumX > 20: 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)): unique, c = numpy.unique(list_to_plot2[i], return_counts=True) new_c = [] new_unique = [] for t in ticks: if t not in unique: new_c.append(0) # add zero count of not occuring new_unique.append(t) else: c_idx = numpy.where(t == unique)[0] new_c.append(c[c_idx]) new_unique.append(unique[c_idx]) y = numpy.array(new_unique) * numpy.array(new_c) if len([list_to_plot_original > 20]) > 0: y[len(y) - 1] = sum(list_to_plot_original[i][list_to_plot_original[i] > 20]) reads.append(y) reads_rel.append(list(numpy.float_(y)) / sum(y)) x = list(numpy.arange(numpy.amin(unique), numpy.amax(unique) + 1).astype(float)) 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) #print(label[i]) #print(colors[i]) #print(w) # linewidth=1, #print(numpy.sum(y)) #print(float(numpy.sum(y))) #print(x) #print(new_y.shape) #new_y = numpy.array(new_y).reshape((len(new_y))) #print(new_y.shape) #print(new_y) #new_y_heights = [x[0] for x in new_y] #print(new_y_heights) #print(new_y.shape) if to_plot[l] == "Relative frequencies": print("relative") new_y = list(numpy.concatenate((numpy.array([yi / float(numpy.sum(y)) for yi in y])))) counts2_rel = ax.bar(x, new_y, align="edge", width=w, edgecolor="black", label=label[i], alpha=0.8, color=colors[i]) else: y = list(y) counts2 = ax.bar(x, y, align="edge", width=w, edgecolor="black", label=label[i], alpha=0.8, color=colors[i]) if i == len(list_to_plot2): 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) if len(list_to_plot2) == 1: ax.set_xticks(numpy.array([xi + 0.2 for xi in x])) 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) pdf.savefig(fig2) plt.close() # write data to CSV file tags output_file.write("Values from family size distribution with all datasets (tags)\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]: 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)) 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)) # 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("\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" else: fs = "={}".format(fs) output_file.write("FS{}{}".format(fs, sep)) if len(label) == 1: output_file.write("{}{}".format(int(reads[0][j]), sep)) else: for n in range(len(label)): output_file.write("{}{}".format(int(reads[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(numpy.concatenate(reads))), sep)) else: for i in reads: output_file.write("{}{}".format(int(sum(i)), sep)) output_file.write("\n") # Family size distribution after DCS and SSCS for dataset, data_o, name_file in zip(list_to_plot, data_array_list, label): maximumX = numpy.amax(dataset) minimumX = numpy.amin(dataset) tags = numpy.array(data_o[:, 2]) seq = numpy.array(data_o[:, 1]) data = numpy.array(dataset) data_o = numpy.array(data_o[:, 0]).astype(int) # find all unique tags and get the indices for ALL tags, but only once u, index_unique, c = numpy.unique(numpy.array(seq), return_counts=True, return_index=True) d = u[c > 1] # get family sizes, tag for duplicates duplTags_double = data[numpy.in1d(seq, d)] duplTags_double_o = data_o[numpy.in1d(seq, d)] duplTags = duplTags_double[0::2] # ab of DCS duplTags_o = duplTags_double_o[0::2] # ab of DCS duplTagsBA = duplTags_double[1::2] # ba of DCS duplTagsBA_o = duplTags_double_o[1::2] # ba of DCS # duplTags_double_tag = tags[numpy.in1d(seq, d)] # duplTags_double_seq = seq[numpy.in1d(seq, d)] # get family sizes for SSCS with no partner ab = numpy.where(tags == "ab")[0] abSeq = seq[ab] ab_o = data_o[ab] ab = data[ab] ba = numpy.where(tags == "ba")[0] baSeq = seq[ba] ba_o = data_o[ba] ba = data[ba] dataAB = ab[numpy.in1d(abSeq, d, invert=True)] dataAB_o = ab_o[numpy.in1d(abSeq, d, invert=True)] dataBA = ba[numpy.in1d(baSeq, d, invert=True)] dataBA_o = ba_o[numpy.in1d(baSeq, d, invert=True)] list1 = [duplTags_double, dataAB, dataBA] # list for plotting # information for family size >= 3 dataAB_FS3 = dataAB[dataAB >= 3] dataAB_FS3_o = dataAB_o[dataAB_o >= 3] dataBA_FS3 = dataBA[dataBA >= 3] dataBA_FS3_o = dataBA_o[dataBA_o >= 3] # ab_FS3 = ab[ab >= 3] # ba_FS3 = ba[ba >= 3] # ab_FS3_o = ab_o[ab_o >= 3] # ba_FS3_o = ba_o[ba_o >= 3] duplTags_FS3 = duplTags[(duplTags >= 3) & (duplTagsBA >= 3)] # ab+ba with FS>=3 duplTags_FS3_BA = duplTagsBA[(duplTags >= 3) & (duplTagsBA >= 3)] # ba+ab with FS>=3 duplTags_double_FS3 = len(duplTags_FS3) + len(duplTags_FS3_BA) # both ab and ba strands with FS>=3 # original FS duplTags_FS3_o = duplTags_o[(duplTags_o >= 3) & (duplTagsBA_o >= 3)] # ab+ba with FS>=3 duplTags_FS3_BA_o = duplTagsBA_o[(duplTags_o >= 3) & (duplTagsBA_o >= 3)] # ba+ab with FS>=3 duplTags_double_FS3_o = sum(duplTags_FS3_o) + sum(duplTags_FS3_BA_o) # both ab and ba strands with FS>=3 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(minimumX, maximumX + 2), 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) if maximumX > 20: 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 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.xlabel("Family size", fontsize=14) plt.ylabel("Absolute Frequency", fontsize=14) plt.margins(0.01, None) plt.grid(b=True, which="major", color="#424242", linestyle=":") # 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))) 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))) 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.astype(int), last.astype(int)) plt.text(0.23, 0.03, legend4, size=10, transform=plt.gcf().transFigure) legend3 = "{:.3f}\n{:.3f}".format(singl / len(data), 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(fig) plt.close() # write same information to a csv file count = numpy.bincount(integers) # original counts of family sizes output_file.write("\nDataset:{}{}\n".format(sep, name_file)) output_file.write("max. family size:{}{}\n".format(sep, max(integers))) output_file.write("absolute frequency:{}{}\n".format(sep, count[len(count) - 1])) output_file.write("relative frequency:{}{:.3f}\n\n".format(sep, float(count[len(count) - 1]) / sum(count))) 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))) # 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("{}{}{}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, 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, 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("total nr. of tags{}{}{}{}{}{}{}{}{}{}{}{}\n".format( 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("{}{}{}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, float(len(dataAB_FS3)) / (len(dataAB_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), sep, float(len(dataAB_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + duplTags_double_FS3), sep, float(sum(dataAB_FS3_o)) / (sum(dataAB_FS3_o) + sum(dataBA_FS3_o) + sum(duplTags_FS3_o)), sep, float(sum(dataAB_FS3_o)) / (sum(dataBA_FS3_o) + sum(dataBA_FS3_o) + duplTags_double_FS3_o))) output_file.write("SSCS ba{}{}{}{}{}{:.3f}{}{:.3f}{}{:.3f}{}{:.3f}\n".format( sep, len(dataBA_FS3), sep, sum(dataBA_FS3_o), sep, float(len(dataBA_FS3)) / (len(dataBA_FS3) + len(dataBA_FS3) + len(duplTags_FS3)), 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("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))) output_file.write("\nValues from family size distribution\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]): 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))) print("Files successfully created!") if __name__ == '__main__': sys.exit(compare_read_families(sys.argv))