comparison hd.py @ 25:9e384b0741f1 draft

planemo upload for repository https://github.com/monikaheinzl/duplexanalysis_galaxy/tree/master/tools/hd commit b8a2f7b7615b2bcd3b602027af31f4e677da94f6-dirty

24:3bc67ac46740

# and finally a CSV file with the data of the plots.
# It is also possible to perform the HD analysis with shortened tags with given sizes as input.
# The tool can run on a certain number of processors, which can be defined by the user.

# USAGE: python hd.py --inputFile filename --inputName1 filename --sample_size int /
#        --only_DCS True --FamilySize3 True --subset_tag True --nproc int --minFS int --maxFS int --nr_above_bars True/False --output_pdf outputfile_name_pdf --output_tabular outputfile_name_tabular --output_chimeras_tabular outputfile_name_chimeras_tabular

import argparse
import itertools
import operator
import sys
from collections import Counter, defaultdict
from functools import partial
from multiprocessing.pool import Pool

import matplotlib.pyplot as plt
import numpy
from matplotlib.backends.backend_pdf import PdfPages

            else:
                plt.annotate("{:,}\n{:.3f}".format(x_label, float(x_label) / sum(counts), 1),
                             xy=(label, x_label + len(con_list1) * 0.01),
                             xycoords="data", color="#000066", fontsize=10)

    legend = "sample size= {:,} against {:,}".format(sum(counts), lenTags)
    plt.text(0.14, -0.01, legend, size=12, transform=plt.gcf().transFigure)
    if nr_unique_chimeras != 0 and len_sample != 0:
        if relative == True:
            legend = "nr. of unique chimeric tags= {:,} ({:.5f}) (rel.diff=1)".format(nr_unique_chimeras,
                                                                         int(nr_unique_chimeras) / float(len_sample))
        else:
            legend = "nr. of unique chimeric tags= {:,} ({:.5f})".format(nr_unique_chimeras, int(nr_unique_chimeras) / float(len_sample))
        plt.text(0.14, -0.05, legend, size=12, transform=plt.gcf().transFigure)

    pdf.savefig(fig, bbox_inches="tight")
    plt.close("all")
    plt.clf()


def plotHDwithinSeq_Sum2(sum1, sum1min, sum2, sum2min, min_value, lenTags, title_file1, pdf):
    fig = plt.figure(figsize=(6, 8))
    plt.subplots_adjust(bottom=0.1)

    ham_partial = [sum1, sum1min, sum2, sum2min, numpy.array(min_value)]  # new hd within tags

    if len(range(minimumX, maximumX)) == 0:
        range1 = minimumX
    else:
        range1 = range(minimumX, maximumX + 2)

    plt.hist(ham_partial, align="left", rwidth=0.8, stacked=False, label=[ "HD a", "HD b'", "HD b", "HD a'", "HD a+b"], bins=range1, color=["#58ACFA", "#0404B4", "#FE642E", "#B40431", "#585858"], edgecolor='black', linewidth=1)

    plt.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(1.55, 1))
    plt.suptitle('Hamming distances within tags', fontsize=14)
    # plt.title(title_file1, fontsize=12)
    plt.xlabel("HD", fontsize=14)

    plt.grid(b=True, which='major', color='#424242', linestyle=':')

    plt.axis((minimumX - 1, maximumX + 1, 0, maximumY * 1.2))
    plt.xticks(numpy.arange(0, maximumX + 1, 1.0))
    # plt.ylim(0, maximumY * 1.2)

    legend = "sample size= {:,} against {:,}".format(len(numpy.concatenate(ham_partial)), lenTags)
    plt.text(0.14, -0.01, legend, size=12, transform=plt.gcf().transFigure)
    pdf.savefig(fig, bbox_inches="tight")
    plt.close("all")
    plt.clf()

    elif mate_b is True:  # HD calculation for all b's
        half1_mate1 = array1_half2
        half2_mate1 = array1_half
        half1_mate2 = array2_half2
        half2_mate2 = array2_half

    for a, b, tag in zip(half1_mate1, half2_mate1, array1):
        # exclude identical tag from array2, to prevent comparison to itself
        sameTag = numpy.where(array2 == tag)[0]
        indexArray2 = numpy.arange(0, len(array2), 1)

        array2_withoutSame = array2[index_withoutSame]  # whole tag (=not splitted into 2 halfs)

        dist = numpy.array([sum(itertools.imap(operator.ne, a, c)) for c in
                            array2_half_withoutSame])  # calculate HD of "a" in the tag to all "a's" or "b" in the tag to all "b's"
        min_index = numpy.where(dist == dist.min())[0]  # get index of min HD
        min_value = dist[min_index]  # get minimum HDs
        min_tag_half2 = array2_half2_withoutSame[min_index]  # get all "b's" of the tag or all "a's" of the tag with minimum HD
        min_tag_array2 = array2_withoutSame[min_index]  # get whole tag with min HD

        dist_second_half = numpy.array([sum(itertools.imap(operator.ne, b, e)) for e in
                             min_tag_half2])  # calculate HD of "b" to all "b's" or "a" to all "a's"
        dist2 = [dist_second_half.max()]
        min_value = [dist.min()]
        max_index = numpy.where(dist_second_half == dist_second_half.max())[0]  # get index of max HD
        max_tag = min_tag_array2[max_index]

        for d, d2 in zip(min_value, dist2):
            if mate_b is True:  # half2, corrects the variable of the HD from both halfs if it is a or b
                ham2.append(d)
                ham2min.append(d2)
            else:  # half1, corrects the variable of the HD from both halfs if it is a or b
                ham1.append(d)
                ham1min.append(d2)

            min_valueList.append(d + d2)
            min_tagsList.append(tag)
            difference1 = abs(d - d2)
            diff11.append(difference1)
            rel_difference = round(float(difference1) / (d + d2), 1)
            relativeDiffList.append(rel_difference)

            # tags which have identical parts:
            if d == 0 or d2 == 0:
                min_tagsList_zeros.append(numpy.array(tag))
                difference1_zeros = abs(d - d2)  # hd of non-identical part
                diff11_zeros.append(difference1_zeros)
                max_tag_list.append(numpy.array(max_tag))

        i += 1

    # print(i)
    # diff11 = [st for st in diff11 if st != 999]

    parser.add_argument('--minFS', default=1, type=int,
                        help='Only tags, which have a family size greater or equal than specified, are included in the HD analysis')
    parser.add_argument('--maxFS', default=0, type=int,
                        help='Only tags, which have a family size smaller or equal than specified, are included in the HD analysis')
    parser.add_argument('--nr_above_bars', action="store_true",
                        help='If no, values above bars in the histrograms are removed')

    parser.add_argument('--output_tabular', default="data.tabular", type=str,
                        help='Name of the tabular file.')
    parser.add_argument('--output_pdf', default="data.pdf", type=str,
                        help='Name of the pdf file.')
    parser.add_argument('--output_chimeras_tabular', default="data_chimeras.tabular", type=str,
                        help='Name of the tabular file with all chimeric tags.')
    return parser


def Hamming_Distance_Analysis(argv):
    parser = make_argparser()
    args = parser.parse_args(argv[1:])
    file1 = args.inputFile
    name1 = args.inputName1
    index_size = args.sample_size
    title_savedFile_pdf = args.output_pdf
    title_savedFile_csv = args.output_tabular
    output_chimeras_tabular = args.output_chimeras_tabular

    sep = "\t"
    onlyDuplicates = args.only_DCS
    minFS = args.minFS
    maxFS = args.maxFS
    nr_above_bars = args.nr_above_bars
    subset = args.subset_tag
    nproc = args.nproc

    # input checks
    if index_size < 0:

    with open(title_savedFile_csv, "w") as output_file, PdfPages(title_savedFile_pdf) as pdf:
        print("dataset: ", name1)
        integers, data_array = readFileReferenceFree(file1)
        data_array = numpy.array(data_array)
        print("total nr of tags:", len(data_array))
        n = [i for i, x in enumerate(data_array[:, 1]) if "N" in x]
        if len(n) != 0:  # delete tags with N in the tag from data
            print("nr of tags with N's within tag:", len(n), float(len(n)) / len(data_array))
            index_whole_array = numpy.arange(0, len(data_array), 1)
            index_withoutN_inTag = numpy.delete(index_whole_array, n)
            data_array = data_array[index_withoutN_inTag, :]
            integers = integers[index_withoutN_inTag]
            print("total nr of tags without Ns:", len(data_array))

        data_array_whole_dataset = data_array

        int_f = numpy.array(data_array[:, 0]).astype(int)
        data_array = data_array[numpy.where(int_f >= minFS)]
        integers = integers[integers >= minFS]

        # select family size for tags
        if maxFS > 0:
            int_f2 = numpy.array(data_array[:, 0]).astype(int)
            data_array = data_array[numpy.where(int_f2 <= maxFS)]
            integers = integers[integers <= maxFS]

        print("min FS", min(integers))
        print("max FS", max(integers))

        tags = data_array[:, 2]
        seq = data_array[:, 1]

        if onlyDuplicates is True:
            # 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

            duplTagsBA = duplTags_double[1::2]  # ba of DCS

            duplTags_tag = tags[numpy.in1d(seq, d)][0::2]  # ab
            duplTags_seq = seq[numpy.in1d(seq, d)][0::2]  # ab - tags

            data_array = numpy.column_stack((duplTags, duplTags_seq))
            data_array = numpy.column_stack((data_array, duplTags_tag))
            integers = numpy.array(data_array[:, 0]).astype(int)
            print("DCS in whole dataset", len(data_array))

        # HD analysis for a subset of the tag
        if subset > 0:
            tag1 = numpy.array([i[0:(len(i)) / 2] for i in data_array[:, 1]])
            tag2 = numpy.array([i[len(i) / 2:len(i)] for i in data_array[:, 1]])

                    [i[flanking_region:len(i) - flanking_region_rounded_end] for i in tag2])

            data_array_tag = numpy.array([i + j for i, j in zip(tag1_shorten, tag2_shorten)])
            data_array = numpy.column_stack((data_array[:, 0], data_array_tag, data_array[:, 2]))

        print("length of tag:", len(data_array[0, 1]))
        # select sample: if no size given --> all vs. all comparison
        if index_size == 0:
            result = numpy.arange(0, len(data_array), 1)
        else:
            result = numpy.random.choice(len(integers), size=index_size,
                                         replace=False)  # array of random sequences of size=index.size
            # unique_tags, unique_indices = numpy.unique(data_array[:, 1], return_index=True) # get only unique tags
            # result = numpy.random.choice(unique_indices, size=index_size,
            #                        replace=False)  # array of random sequences of size=index.size

        # with open("index_result1_{}.pkl".format(app_f), "wb") as o:
        #     pickle.dump(result, o, pickle.HIGHEST_PROTOCOL)

        # comparison random tags to whole dataset
        result1 = data_array[result, 1]  # random tags
        result2 = data_array[:, 1]  # all tags
        print("sample size:", len(result1))

        # HD analysis of whole tag
        proc_pool = Pool(nproc)
        chunks_sample = numpy.array_split(result1, nproc)
        ham = proc_pool.map(partial(hamming, array2=result2), chunks_sample)

        ham = numpy.concatenate(ham).astype(int)
        # with open("HD_whole dataset_{}.txt".format(app_f), "w") as output_file1:
        # for h, tag in zip(ham, result1):
        #     output_file1.write("{}\t{}\n".format(tag, h))

        # HD analysis for chimeric reads
        # result2 = data_array_whole_dataset[:,1]

        proc_pool_b = Pool(nproc)
        diff_list_a = proc_pool_b.map(partial(hamming_difference, array2=result2, mate_b=False), chunks_sample)
        diff_list_b = proc_pool_b.map(partial(hamming_difference, array2=result2, mate_b=True), chunks_sample)
        proc_pool_b.close()
        proc_pool_b.join()
        diff = numpy.concatenate((numpy.concatenate([item[0] for item in diff_list_a]),
                                  numpy.concatenate([item_b[0] for item_b in diff_list_b]))).astype(int)
        HDhalf1 = numpy.concatenate((numpy.concatenate([item[1] for item in diff_list_a]),
                                     numpy.concatenate([item_b[1] for item_b in diff_list_b]))).astype(int)
        HDhalf2 = numpy.concatenate((numpy.concatenate([item[2] for item in diff_list_a]),
                                     numpy.concatenate([item_b[2] for item_b in diff_list_b]))).astype(int)
        minHDs = numpy.concatenate((numpy.concatenate([item[3] for item in diff_list_a]),
                                    numpy.concatenate([item_b[3] for item_b in diff_list_b]))).astype(int)
        minHD_tags = numpy.concatenate((numpy.concatenate([item[4] for item in diff_list_a]),
                                        numpy.concatenate([item_b[4] for item_b in diff_list_b])))
        rel_Diff = numpy.concatenate((numpy.concatenate([item[5] for item in diff_list_a]),
                                      numpy.concatenate([item_b[5] for item_b in diff_list_b])))
        diff_zeros = numpy.concatenate((numpy.concatenate([item[6] for item in diff_list_a]),
                                        numpy.concatenate([item_b[6] for item_b in diff_list_b]))).astype(int)
        minHD_tags_zeros = numpy.concatenate((numpy.concatenate([item[7] for item in diff_list_a]),
                                              numpy.concatenate([item_b[7] for item_b in diff_list_b])))
        HDhalf1min = numpy.concatenate((numpy.concatenate([item[8] for item in diff_list_a]),
                                        numpy.concatenate([item_b[8] for item_b in diff_list_b]))).astype(int)
        HDhalf2min = numpy.concatenate((numpy.concatenate([item[9] for item in diff_list_a]),
                                        numpy.concatenate([item_b[9] for item_b in diff_list_b]))).astype(int)

        chimera_tags = numpy.concatenate(([item[10] for item in diff_list_a],
                                       [item_b[10] for item_b in diff_list_b]))

        chimera_tags = [x for x in chimera_tags if x != []]
        chimera_tags_new = []

        for i in chimera_tags:
            if len(i) > 1:
                for t in i:
                    chimera_tags_new.append(t)
            else:
                chimera_tags_new.extend(i)

        chimeras_dic = defaultdict(list)
        for t1, t2 in zip(minHD_tags_zeros, chimera_tags_new):
            chimeras_dic[t1].append(t2)

        lst_unique_chimeras = []
        with open(output_chimeras_tabular, "w") as output_file1:
            unique_chimeras = numpy.unique(minHD_tags_zeros)
            sample_half_a = numpy.array([i[0:(len(i)) / 2] for i in unique_chimeras])  # mate1 part1
            sample_half_b = numpy.array([i[len(i) / 2:len(i)] for i in unique_chimeras])  # mate1 part 2

            output_file1.write("sample tag\tsimilar tag\n")
            for tag1, a, b in zip(unique_chimeras, sample_half_a, sample_half_b):
                max_tags = numpy.concatenate(chimeras_dic.get(tag1))

                if tag1 in chimeras_dic.values():
                    continue
                else:
                    lst_unique_chimeras.append(tag1)

                chimera_half_a = numpy.array([i[0:(len(i)) / 2] for i in max_tags])  # mate1 part1
                chimera_half_b = numpy.array([i[len(i) / 2:len(i)] for i in max_tags])  # mate1 part 2

                new_format = []
                for i in range(len(max_tags)):
                    if a == chimera_half_a[i]:
                        max_tag = "*{}* {}".format(chimera_half_a[i], chimera_half_b[i])
                        new_format.append(max_tag)

                    elif b == chimera_half_b[i]:
                        max_tag = "{} *{}*".format(chimera_half_a[i], chimera_half_b[i])
                        new_format.append(max_tag)

                sample_tag = "{} {}".format(a, b)
                output_file1.write("{}\t{}\n".format(sample_tag, ", ".join(new_format)))
            output_file1.write(
                "This file contains all tags that were identified as chimeras as the first column and the corresponding tags which returned a Hamming distance of zero in either the first or the second half of the sample tag as the second column.\n "
                "The tags were separated by an empty space into their halves and the * marks the identical half.")

        nr_chimeric_tags = len(lst_unique_chimeras)
        print("nr of unique chimeras:", nr_chimeric_tags)

        lenTags = len(data_array)
        len_sample = len(result1)

        quant = numpy.array(data_array[result, 0]).astype(int)  # family size for sample of tags

            for s, q in zip(seq, quant):
                seqDic[s].append(q)
        else:
            seqDic = dict(zip(seq, quant))


        lst_minHD_tags = []
        for i in minHD_tags:
            lst_minHD_tags.append(seqDic.get(i))

        if onlyDuplicates:
            lst_minHD_tags = numpy.concatenate(([item[0] for item in lst_minHD_tags],
                                        [item_b[1] for item_b in lst_minHD_tags])).astype(int)
        # else:
        #     lst_minHD_tags = numpy.concatenate(lst_minHD_tags)

        # histogram with absolute and relative difference between HDs of both parts of the tag
        listDifference1, maximumXDifference, minimumXDifference = hammingDistanceWithFS(lst_minHD_tags, diff)
        listRelDifference1, maximumXRelDifference, minimumXRelDifference = hammingDistanceWithFS(lst_minHD_tags,
                                                                                                 rel_Diff)

        # family size distribution separated after the difference between HDs of both parts of the tag
        # familySizeList1_diff, hammingDistances_diff, maximumXFS_diff, minimumXFS_diff = familySizeDistributionWithHD(
        #   lst_minHD_tags, diff, diff=True, rel=False)
        # familySizeList1_reldiff, hammingDistances_reldiff, maximumXFS_reldiff, minimumXFS_reldiff = familySizeDistributionWithHD(
        #    lst_minHD_tags, rel_Diff, diff=True, rel=True)

        # chimeric read analysis: tags which have HD=0 in one of the halfs
        if len(minHD_tags_zeros) != 0:
            lst_minHD_tags_zeros = []
            for i in minHD_tags_zeros:
                lst_minHD_tags_zeros.append(seqDic.get(i))  # get family size for tags of chimeric reads
            if onlyDuplicates:
                lst_minHD_tags_zeros = numpy.concatenate(([item[0] for item in lst_minHD_tags_zeros],
                                                [item_b[1] for item_b in lst_minHD_tags_zeros])).astype(int)
            
            # histogram with HD of non-identical half
            listDifference1_zeros, maximumXDifference_zeros, minimumXDifference_zeros = hammingDistanceWithFS(
            lst_minHD_tags_zeros, diff_zeros)
            
        # plot Hamming Distance with Family size distribution
        plotHDwithFSD(list1=list1, maximumX=maximumX, minimumX=minimumX, pdf=pdf,
                      subtitle="Hamming distance separated by family size", title_file1=name1, lenTags=lenTags,
                      xlabel="HD", nr_above_bars=nr_above_bars)

        # Plot FSD with separation after
        plotFSDwithHD2(familySizeList1, maximumXFS, minimumXFS,
                       originalCounts=quant, subtitle="Family size distribution separated by Hamming distance",
                       pdf=pdf, relative=False, title_file1=name1, diff=False)

        # Plot HD within tags
        plotHDwithinSeq_Sum2(HDhalf1, HDhalf1min, HDhalf2, HDhalf2min, minHDs, pdf=pdf, lenTags=lenTags,
                             title_file1=name1)

        # Plot difference between HD's separated after FSD
        plotHDwithFSD(listDifference1, maximumXDifference, minimumXDifference, pdf=pdf,
                      subtitle="Delta Hamming distance within tags",
                      title_file1=name1, lenTags=lenTags,
                      xlabel="absolute delta HD", relative=False, nr_above_bars=nr_above_bars)

        plotHDwithFSD(listRelDifference1, maximumXRelDifference, minimumXRelDifference, pdf=pdf,
                      subtitle="Chimera Analysis: relative delta Hamming distances",
                      title_file1=name1, lenTags=lenTags,
                      xlabel="relative delta HD", relative=True, nr_above_bars=nr_above_bars, nr_unique_chimeras=nr_chimeric_tags, len_sample=len_sample)

        # plots for chimeric reads
        if len(minHD_tags_zeros) != 0:
            # HD
            plotHDwithFSD(listDifference1_zeros, maximumXDifference_zeros, minimumXDifference_zeros, pdf=pdf,subtitle="Hamming distance of the non-identical half of chimeras",
                          title_file1=name1, lenTags=lenTags, xlabel="HD", relative=False,
                          nr_above_bars=nr_above_bars, nr_unique_chimeras=nr_chimeric_tags, len_sample=len_sample)

        # print all data to a CSV file
        # HD

        output_file.write(
            "relative frequency:{}{}\n\n".format(sep, float(max_fs[len(max_fs) - 1]) / sum(max_fs)))

        # HD within tags
        output_file.write(
            "The hamming distances were calculated by comparing each half of all tags against the tag(s) with the minimum Hamming distance per half.\n"
            "Since this calculation was repeated, but starting with the second half to find all possible chimeras in the data, the actual number of tags in the plots differs from the sample size entered by the user.\n"
            "In addition, both family sizes of one tag will be included in the plots if only tags of reads that can form a DCS were allowed.\n")

        output_file.write("length of one part of the tag = {}\n\n".format(len(data_array[0, 1]) / 2))

        createFileHDwithinTag(summary9, sumCol9, overallSum9, output_file,
                              "Hamming distance of each half in the tag", sep)

        createFileHD(summary13, sumCol13, overallSum13, output_file,
                     "Chimera analysis: relative delta Hamming distances", sep)

        if len(minHD_tags_zeros) != 0:
            output_file.write(
                "Chimeras:\nAll tags were filtered: only those tags where at least one half is identical with the half of the min. tag are kept.\nSo the Hamming distance of the non-identical half is shown.\n")
            output_file.write(
                "Be aware that the real number of chimeric tags (where rel. diff = 1) is not shown in the plot because of the above reasons.\n")
            output_file.write("real number of chimeric tags{}{}{}{}\n".format(sep, nr_chimeric_tags, sep, int(nr_chimeric_tags) / float(len_sample)))
            createFileHD(summary15, sumCol15, overallSum15, output_file,
                         "Hamming distances of non-zero half", sep)

        output_file.write("\n")


if __name__ == '__main__':
    sys.exit(Hamming_Distance_Analysis(sys.argv))

25:9e384b0741f1

# and finally a CSV file with the data of the plots.
# It is also possible to perform the HD analysis with shortened tags with given sizes as input.
# The tool can run on a certain number of processors, which can be defined by the user.

# USAGE: python hd.py --inputFile filename --inputName1 filename --sample_size int /
#        --only_DCS True --FamilySize3 True --subset_tag True --nproc int --minFS int --maxFS int --nr_above_bars True/False --output_tabular outptufile_name_tabular

import argparse
import itertools
import operator
import sys
from collections import Counter, defaultdict
from functools import partial
from multiprocessing.pool import Pool
import random
import os

import matplotlib.pyplot as plt
import numpy
from matplotlib.backends.backend_pdf import PdfPages

            else:
                plt.annotate("{:,}\n{:.3f}".format(x_label, float(x_label) / sum(counts), 1),
                             xy=(label, x_label + len(con_list1) * 0.01),
                             xycoords="data", color="#000066", fontsize=10)

    legend = "nr. of tags = {:,}\nsample size = {:,}\nnr. of data points = {:,}".format(lenTags, len_sample, sum(counts))
    plt.text(0.14, -0.05, legend, size=12, transform=plt.gcf().transFigure)

    # if nr_unique_chimeras != 0 and len_sample != 0:
    #     if relative == True:
    #         legend = "nr. of unique chimeric tags= {:,} ({:.5f}) (rel.diff=1)".format(nr_unique_chimeras,
    #                                                                      int(nr_unique_chimeras) / float(len_sample))
    #     else:
    #         legend = "nr. of unique chimeric tags= {:,} ({:.5f})".format(nr_unique_chimeras, int(nr_unique_chimeras) / float(len_sample))
    #     plt.text(0.14, -0.09, legend, size=12, transform=plt.gcf().transFigure)

    pdf.savefig(fig, bbox_inches="tight")
    plt.close("all")
    plt.clf()


def plotHDwithinSeq_Sum2(sum1, sum1min, sum2, sum2min, min_value, lenTags, title_file1, pdf, len_sample):
    fig = plt.figure(figsize=(6, 8))
    plt.subplots_adjust(bottom=0.1)

    ham_partial = [sum1, sum1min, sum2, sum2min, numpy.array(min_value)]  # new hd within tags

    if len(range(minimumX, maximumX)) == 0:
        range1 = minimumX
    else:
        range1 = range(minimumX, maximumX + 2)

    plt.hist(ham_partial, align="left", rwidth=0.8, stacked=False, label=["HD a", "HD b'", "HD b", "HD a'", "HD a+b"], bins=range1, color=["#58ACFA", "#0404B4", "#FE642E", "#B40431", "#585858"], edgecolor='black', linewidth=1)

    plt.legend(loc='upper right', fontsize=14, frameon=True, bbox_to_anchor=(1.55, 1))
    plt.suptitle('Hamming distances within tags', fontsize=14)
    # plt.title(title_file1, fontsize=12)
    plt.xlabel("HD", fontsize=14)

    plt.grid(b=True, which='major', color='#424242', linestyle=':')

    plt.axis((minimumX - 1, maximumX + 1, 0, maximumY * 1.2))
    plt.xticks(numpy.arange(0, maximumX + 1, 1.0))
    # plt.ylim(0, maximumY * 1.2)
    legend = "nr. of tags = {:,}\nsample size = {:,}\nnr. of data points = {:,}".format(lenTags, len_sample, len(numpy.concatenate(ham_partial)))

    # legend = "sample size= {:,} against {:,}".format(len(numpy.concatenate(ham_partial)), lenTags)
    plt.text(0.14, -0.05, legend, size=12, transform=plt.gcf().transFigure)
    pdf.savefig(fig, bbox_inches="tight")
    plt.close("all")
    plt.clf()

    elif mate_b is True:  # HD calculation for all b's
        half1_mate1 = array1_half2
        half2_mate1 = array1_half
        half1_mate2 = array2_half2
        half2_mate2 = array2_half
    # half1_mate1, index_halves = numpy.unique(half1_mate1, return_index=True)
    # print(len(half1_mate1))
    # half2_mate1 = half2_mate1[index_halves]
    # array1 = array1[index_halves]

    for a, b, tag in zip(half1_mate1, half2_mate1, array1):
        # exclude identical tag from array2, to prevent comparison to itself
        sameTag = numpy.where(array2 == tag)[0]
        indexArray2 = numpy.arange(0, len(array2), 1)

        array2_withoutSame = array2[index_withoutSame]  # whole tag (=not splitted into 2 halfs)

        dist = numpy.array([sum(itertools.imap(operator.ne, a, c)) for c in
                            array2_half_withoutSame])  # calculate HD of "a" in the tag to all "a's" or "b" in the tag to all "b's"
        min_index = numpy.where(dist == dist.min())[0]  # get index of min HD
        min_value = dist.min()
        # min_value = dist[min_index]  # get minimum HDs
        min_tag_half2 = array2_half2_withoutSame[min_index]  # get all "b's" of the tag or all "a's" of the tag with minimum HD
        min_tag_array2 = array2_withoutSame[min_index]  # get whole tag with min HD

        dist_second_half = numpy.array([sum(itertools.imap(operator.ne, b, e)) for e in min_tag_half2])  # calculate HD of "b" to all "b's" or "a" to all "a's"
        max_value = dist_second_half.max()
        max_index = numpy.where(dist_second_half == dist_second_half.max())[0]  # get index of max HD
        max_tag = min_tag_array2[max_index]

        # for d, d2 in zip(min_value, max_value):
        if mate_b is True:  # half2, corrects the variable of the HD from both halfs if it is a or b
            ham2.append(min_value)
            ham2min.append(max_value)
        else:  # half1, corrects the variable of the HD from both halfs if it is a or b
            ham1.append(min_value)
            ham1min.append(max_value)

        min_valueList.append(min_value + max_value)
        min_tagsList.append(tag)
        difference1 = abs(min_value - max_value)
        diff11.append(difference1)
        rel_difference = round(float(difference1) / (min_value + max_value), 1)
        relativeDiffList.append(rel_difference)

        # tags which have identical parts:
        if min_value == 0 or max_value == 0:
            min_tagsList_zeros.append(numpy.array(tag))
            difference1_zeros = abs(min_value - max_value)  # hd of non-identical part
            diff11_zeros.append(difference1_zeros)
            max_tag_list.append(max_tag)
        else:
            min_tagsList_zeros.append(None)
            diff11_zeros.append(None)
            max_tag_list.append(numpy.array(["None"]))

            # max_tag_list.append(numpy.array(max_tag))

        i += 1

    # print(i)
    # diff11 = [st for st in diff11 if st != 999]

    parser.add_argument('--minFS', default=1, type=int,
                        help='Only tags, which have a family size greater or equal than specified, are included in the HD analysis')
    parser.add_argument('--maxFS', default=0, type=int,
                        help='Only tags, which have a family size smaller or equal than specified, are included in the HD analysis')
    parser.add_argument('--nr_above_bars', action="store_true",
                        help='If no, values above bars in the histograms are removed')

    parser.add_argument('--output_tabular', default="data.tabular", type=str,
                        help='Name of the tabular file.')
    parser.add_argument('--output_pdf', default="data.pdf", type=str,
                        help='Name of the pdf file.')
    parser.add_argument('--output_chimeras_tabular', default="data.tabular", type=str,
                        help='Name of the tabular file with all chimeric tags.')

    return parser


def Hamming_Distance_Analysis(argv):

    parser = make_argparser()
    args = parser.parse_args(argv[1:])

    file1 = args.inputFile
    name1 = args.inputName1

    index_size = args.sample_size
    title_savedFile_pdf = args.output_pdf
    title_savedFile_csv = args.output_tabular
    output_chimeras_tabular = args.output_chimeras_tabular

    sep = "\t"
    onlyDuplicates = args.only_DCS
    minFS = args.minFS
    maxFS = args.maxFS
    nr_above_bars = args.nr_above_bars

    subset = args.subset_tag
    nproc = args.nproc

    # input checks
    if index_size < 0:

    with open(title_savedFile_csv, "w") as output_file, PdfPages(title_savedFile_pdf) as pdf:
        print("dataset: ", name1)
        integers, data_array = readFileReferenceFree(file1)
        data_array = numpy.array(data_array)
        print("total nr of tags with Ns:", len(data_array))
        n = [i for i, x in enumerate(data_array[:, 1]) if "N" in x]
        if len(n) != 0:  # delete tags with N in the tag from data
            print("nr of tags with N's within tag:", len(n), float(len(n)) / len(data_array))
            index_whole_array = numpy.arange(0, len(data_array), 1)
            index_withoutN_inTag = numpy.delete(index_whole_array, n)
            data_array = data_array[index_withoutN_inTag, :]
            integers = integers[index_withoutN_inTag]
            print("total nr of tags without Ns:", len(data_array))

        int_f = numpy.array(data_array[:, 0]).astype(int)
        data_array = data_array[numpy.where(int_f >= minFS)]
        integers = integers[integers >= minFS]

        # select family size for tags
        if maxFS > 0:
            int_f2 = numpy.array(data_array[:, 0]).astype(int)
            data_array = data_array[numpy.where(int_f2 <= maxFS)]
            integers = integers[integers <= maxFS]

        if onlyDuplicates is True:
            tags = data_array[:, 2]
            seq = data_array[:, 1]

            # 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

            duplTagsBA = duplTags_double[1::2]  # ba of DCS

            duplTags_tag = tags[numpy.in1d(seq, d)][0::2]  # ab
            duplTags_seq = seq[numpy.in1d(seq, d)][0::2]  # ab - tags

            if minFS > 1:
                duplTags_tag = duplTags_tag[(duplTags >= 3) & (duplTagsBA >= 3)]
                duplTags_seq = duplTags_seq[(duplTags >= 3) & (duplTagsBA >= 3)]
                duplTags = duplTags[(duplTags >= 3) & (duplTagsBA >= 3)]  # ab+ba with FS>=3

            data_array = numpy.column_stack((duplTags, duplTags_seq))
            data_array = numpy.column_stack((data_array, duplTags_tag))
            integers = numpy.array(data_array[:, 0]).astype(int)
            print("DCS in whole dataset", len(data_array))

        print("min FS", min(integers))
        print("max FS", max(integers))

        # HD analysis for a subset of the tag
        if subset > 0:
            tag1 = numpy.array([i[0:(len(i)) / 2] for i in data_array[:, 1]])
            tag2 = numpy.array([i[len(i) / 2:len(i)] for i in data_array[:, 1]])

                    [i[flanking_region:len(i) - flanking_region_rounded_end] for i in tag2])

            data_array_tag = numpy.array([i + j for i, j in zip(tag1_shorten, tag2_shorten)])
            data_array = numpy.column_stack((data_array[:, 0], data_array_tag, data_array[:, 2]))

        print("length of tag= ", len(data_array[0, 1]))
        # select sample: if no size given --> all vs. all comparison
        if index_size == 0:
            result = numpy.arange(0, len(data_array), 1)
        else:
            numpy.random.shuffle(data_array)
            unique_tags, unique_indices = numpy.unique(data_array[:, 1], return_index=True)  # get only unique tags
            result = numpy.random.choice(unique_indices, size=index_size, replace=False)  # array of random sequences of size=index.size

            # result = numpy.random.choice(len(integers), size=index_size,
            #                             replace=False)  # array of random sequences of size=index.size
            # result = numpy.where(numpy.array(random_tags) == numpy.array(data_array[:,1]))[0]

        # with open("index_result1_{}.pkl".format(app_f), "wb") as o:
        #     pickle.dump(result, o, pickle.HIGHEST_PROTOCOL)

        # comparison random tags to whole dataset
        result1 = data_array[result, 1]  # random tags
        result2 = data_array[:, 1]  # all tags
        print("sample size= ", len(result1))

        # HD analysis of whole tag
        proc_pool = Pool(nproc)
        chunks_sample = numpy.array_split(result1, nproc)
        ham = proc_pool.map(partial(hamming, array2=result2), chunks_sample)

        ham = numpy.concatenate(ham).astype(int)
        # with open("HD_whole dataset_{}.txt".format(app_f), "w") as output_file1:
        # for h, tag in zip(ham, result1):
        #     output_file1.write("{}\t{}\n".format(tag, h))

        # # HD analysis for chimeric reads
        # result2 = data_array_whole_dataset[:,1]

        proc_pool_b = Pool(nproc)
        diff_list_a = proc_pool_b.map(partial(hamming_difference, array2=result2, mate_b=False), chunks_sample)
        diff_list_b = proc_pool_b.map(partial(hamming_difference, array2=result2, mate_b=True), chunks_sample)
        proc_pool_b.close()
        proc_pool_b.join()
        HDhalf1 = numpy.concatenate((numpy.concatenate([item[1] for item in diff_list_a]),
                                     numpy.concatenate([item_b[1] for item_b in diff_list_b]))).astype(int)
        HDhalf2 = numpy.concatenate((numpy.concatenate([item[2] for item in diff_list_a]),
                                     numpy.concatenate([item_b[2] for item_b in diff_list_b]))).astype(int)
        minHDs = numpy.concatenate((numpy.concatenate([item[3] for item in diff_list_a]),
                                    numpy.concatenate([item_b[3] for item_b in diff_list_b]))).astype(int)
        HDhalf1min = numpy.concatenate((numpy.concatenate([item[8] for item in diff_list_a]),
                                        numpy.concatenate([item_b[8] for item_b in diff_list_b]))).astype(int)
        HDhalf2min = numpy.concatenate((numpy.concatenate([item[9] for item in diff_list_a]),
                                        numpy.concatenate([item_b[9] for item_b in diff_list_b]))).astype(int)

        rel_Diff1 = numpy.concatenate([item[5] for item in diff_list_a])
        rel_Diff2 = numpy.concatenate([item[5] for item in diff_list_b])
        diff1 = numpy.concatenate([item[0] for item in diff_list_a])
        diff2 = numpy.concatenate([item[0] for item in diff_list_b])

        diff_zeros1 = numpy.concatenate([item[6] for item in diff_list_a])
        diff_zeros2 = numpy.concatenate([item[6] for item in diff_list_b])
        minHD_tags = numpy.concatenate([item[4] for item in diff_list_a])
        minHD_tags_zeros1 = numpy.concatenate([item[7] for item in diff_list_a])
        minHD_tags_zeros2 = numpy.concatenate([item[7] for item in diff_list_b])
        chim_tags = [item[10] for item in diff_list_a]
        chim_tags2 = [item[10] for item in diff_list_b]
        chimera_tags1 = [ii if isinstance(i, list) else i for i in chim_tags for ii in i]
        chimera_tags2 = [ii if isinstance(i, list) else i for i in chim_tags2 for ii in i]

        rel_Diff = []
        diff_zeros = []
        minHD_tags_zeros = []
        diff = []
        chimera_tags = []
        for d1, d2, rel1, rel2, zeros1, zeros2, tag1, tag2, ctag1, ctag2 in \
                zip(diff1, diff2, rel_Diff1, rel_Diff2, diff_zeros1, diff_zeros2, minHD_tags_zeros1, minHD_tags_zeros2, chimera_tags1, chimera_tags2):
            rel_Diff.append(max(rel1, rel2))
            diff.append(max(d1, d2))

            if all(i is not None for i in [zeros1, zeros2]):
                diff_zeros.append(max(zeros1, zeros2))
                minHD_tags_zeros.append(str(tag1))
                tags = [ctag1, ctag2]
                chimera_tags.append(tags)
            elif zeros1 is not None and zeros2 is None:
                diff_zeros.append(zeros1)
                minHD_tags_zeros.append(str(tag1))
                chimera_tags.append(ctag1)
            elif zeros1 is None and zeros2 is not None:
                diff_zeros.append(zeros2)
                minHD_tags_zeros.append(str(tag2))
                chimera_tags.append(ctag2)

        chimera_tags_new = chimera_tags
        #data_chimeraAnalysis = numpy.column_stack((minHD_tags_zeros, chimera_tags_new))
        # chimeras_dic = defaultdict(list)
        #
        # for t1, t2 in zip(minHD_tags_zeros, chimera_tags_new):
        #     if len(t2) >1 and type(t2) is not numpy.ndarray:
        #         t2 = numpy.concatenate(t2)
        #     chimeras_dic[t1].append(t2)

        with open(output_chimeras_tabular, "w") as output_file1:
            output_file1.write("chimera tag\tsimilar tag with HD=0\n")
            for i in range(len(minHD_tags_zeros)):
                tag1 = minHD_tags_zeros[i]
                sample_half_a = tag1[0:(len(tag1)) / 2]
                sample_half_b = tag1[len(tag1) / 2:len(tag1)]

                max_tags = chimera_tags_new[i]
                if isinstance(max_tags, list) and len(max_tags) > 1:
                    max_tags = numpy.concatenate(max_tags)
                #if isinstance(max_tags, list): #and type(max_tags) is not numpy.ndarray:
                #    print(max_tags)
                #    max_tags = numpy.concatenate(max_tags)
                max_tags = numpy.unique(max_tags)

                chimera_half_a = numpy.array([i[0:(len(i)) / 2] for i in max_tags])  # mate1 part1
                chimera_half_b = numpy.array([i[len(i) / 2:len(i)] for i in max_tags])  # mate1 part 2

                new_format = []
                for j in range(len(max_tags)):
                    if sample_half_a == chimera_half_a[j]:
                        max_tag = "*{}* {}".format(chimera_half_a[j], chimera_half_b[j])
                        new_format.append(max_tag)

                    elif sample_half_b == chimera_half_b[j]:
                        max_tag = "{} *{}*".format(chimera_half_a[j], chimera_half_b[j])
                        new_format.append(max_tag)

                sample_tag = "{} {}".format(sample_half_a, sample_half_b)
                output_file1.write("{}\t{}\n".format(sample_tag, ", ".join(new_format)))
            output_file1.write(
                "This file contains all tags that were identified as chimeras as the first column and the corresponding tags which returned a Hamming distance of zero in either the first or the second half of the sample tag as the second column.\n "
                "The tags were separated by an empty space into their halves and the * marks the identical half.")

            # unique_chimeras = numpy.array(minHD_tags_zeros)
            #
            # sample_half_a = numpy.array([i[0:(len(i)) / 2] for i in unique_chimeras])  # mate1 part1
            # sample_half_b = numpy.array([i[len(i) / 2:len(i)] for i in unique_chimeras])  # mate1 part 2
            #
            # output_file1.write("sample tag\tsimilar tag\n")
            # for tag1, a, b in zip(unique_chimeras, sample_half_a, sample_half_b):
            #     max_tags = numpy.concatenate(chimeras_dic.get(tag1))
            #     max_tags = numpy.unique(max_tags)
            #
            #     chimera_half_a = numpy.array([i[0:(len(i)) / 2] for i in max_tags])  # mate1 part1
            #     chimera_half_b = numpy.array([i[len(i) / 2:len(i)] for i in max_tags])  # mate1 part 2
            #
            #     new_format = []
            #     for i in range(len(max_tags)):
            #         if a == chimera_half_a[i]:
            #             max_tag = "*{}* {}".format(chimera_half_a[i], chimera_half_b[i])
            #             new_format.append(max_tag)
            #
            #         elif b == chimera_half_b[i]:
            #             max_tag = "{} *{}*".format(chimera_half_a[i], chimera_half_b[i])
            #             new_format.append(max_tag)
            #
            #     sample_tag = "{} {}".format(a, b)
            #     output_file1.write("{}\t{}\n".format(sample_tag, ", ".join(new_format)))
            # output_file1.write(
            #     "This file contains all tags that were identified as chimeras as the first column and the corresponding tags which returned a Hamming distance of zero in either the first or the second half of the sample tag as the second column.\n "
            #     "The tags were separated by an empty space into their halves and the * marks the identical half.")

        nr_chimeric_tags = len(minHD_tags_zeros)
        print("nr of unique chimeras", nr_chimeric_tags)

        lenTags = len(data_array)
        len_sample = len(result1)

        quant = numpy.array(data_array[result, 0]).astype(int)  # family size for sample of tags

            for s, q in zip(seq, quant):
                seqDic[s].append(q)
        else:
            seqDic = dict(zip(seq, quant))

        lst_minHD_tags = []
        for i in minHD_tags:
            lst_minHD_tags.append(seqDic.get(i))

        if onlyDuplicates:
            lst_minHD_tags = numpy.concatenate(([item[0] for item in lst_minHD_tags], [item_b[1] for item_b in lst_minHD_tags])).astype(int)

        # histogram with absolute and relative difference between HDs of both parts of the tag
        listDifference1, maximumXDifference, minimumXDifference = hammingDistanceWithFS(lst_minHD_tags, diff)
        listRelDifference1, maximumXRelDifference, minimumXRelDifference = hammingDistanceWithFS(lst_minHD_tags,
                                                                                                 rel_Diff)
        # chimeric read analysis: tags which have HD=0 in one of the halfs
        if len(minHD_tags_zeros) != 0:
            lst_minHD_tags_zeros = []
            for i in minHD_tags_zeros:
                lst_minHD_tags_zeros.append(seqDic.get(i))  # get family size for tags of chimeric reads
            if onlyDuplicates:
                lst_minHD_tags_zeros = numpy.concatenate(([item[0] for item in lst_minHD_tags_zeros], [item_b[1] for item_b in lst_minHD_tags_zeros])).astype(int)

            # histogram with HD of non-identical half
            listDifference1_zeros, maximumXDifference_zeros, minimumXDifference_zeros = hammingDistanceWithFS(lst_minHD_tags_zeros, diff_zeros)

        # plot Hamming Distance with Family size distribution
        plotHDwithFSD(list1=list1, maximumX=maximumX, minimumX=minimumX, pdf=pdf,
                      subtitle="Hamming distance separated by family size", title_file1=name1, lenTags=lenTags,
                      xlabel="HD", nr_above_bars=nr_above_bars, len_sample=len_sample)

        # Plot FSD with separation after
        plotFSDwithHD2(familySizeList1, maximumXFS, minimumXFS,
                       originalCounts=quant, subtitle="Family size distribution separated by Hamming distance",
                       pdf=pdf, relative=False, title_file1=name1, diff=False)

        # Plot HD within tags
        plotHDwithinSeq_Sum2(HDhalf1, HDhalf1min, HDhalf2, HDhalf2min, minHDs, pdf=pdf, lenTags=lenTags,
                             title_file1=name1, len_sample=len_sample)

        # Plot difference between HD's separated after FSD
        plotHDwithFSD(listDifference1, maximumXDifference, minimumXDifference, pdf=pdf,
                      subtitle="Delta Hamming distance within tags",
                      title_file1=name1, lenTags=lenTags,
                      xlabel="absolute delta HD", relative=False, nr_above_bars=nr_above_bars, len_sample=len_sample)

        plotHDwithFSD(listRelDifference1, maximumXRelDifference, minimumXRelDifference, pdf=pdf,
                      subtitle="Chimera Analysis: relative delta Hamming distances",
                      title_file1=name1, lenTags=lenTags,
                      xlabel="relative delta HD", relative=True, nr_above_bars=nr_above_bars, nr_unique_chimeras=nr_chimeric_tags, len_sample=len_sample)

        # plots for chimeric reads
        if len(minHD_tags_zeros) != 0:
            # HD
            plotHDwithFSD(listDifference1_zeros, maximumXDifference_zeros, minimumXDifference_zeros, pdf=pdf, subtitle="Hamming distance of chimeras",
                          title_file1=name1, lenTags=lenTags, xlabel="HD", relative=False,
                          nr_above_bars=nr_above_bars, nr_unique_chimeras=nr_chimeric_tags, len_sample=len_sample)

        # print all data to a CSV file
        # HD

        output_file.write(
            "relative frequency:{}{}\n\n".format(sep, float(max_fs[len(max_fs) - 1]) / sum(max_fs)))

        # HD within tags
        output_file.write(
            "The Hamming distances were calculated by comparing the first halve against all halves and selected the minimum value (HD a).\n"
            "For the second half of the tag, we compared them against all tags which resulted in the minimum HD of the previous step and selected the maximum value (HD b').\n"
            "Finally, it was possible to calculate the absolute and relative differences between the HDs (absolute and relative delta HD).\n"
            "These calculations were repeated, but starting with the second half in the first step to find all possible chimeras in the data (HD b and HD  For simplicity we used the maximum value between the delta values in the end.\n"
            "When only tags that can form DCS were allowed in the analysis, family sizes for the forward and reverse (ab and ba) will be included in the plots.\n")

        output_file.write("length of one part of the tag = {}\n\n".format(len(data_array[0, 1]) / 2))

        createFileHDwithinTag(summary9, sumCol9, overallSum9, output_file,
                              "Hamming distance of each half in the tag", sep)

        createFileHD(summary13, sumCol13, overallSum13, output_file,
                     "Chimera analysis: relative delta Hamming distances", sep)

        if len(minHD_tags_zeros) != 0:
            output_file.write(
                "Chimeras:\nAll tags were filtered: only those tags where at least one half was identical (HD=0) and therefore, had a relative delta of 1 were kept. These tags are considered as chimeric.\nSo the Hamming distances of the chimeric tags are shown.\n")
            createFileHD(summary15, sumCol15, overallSum15, output_file,
                         "Hamming distances of chimeras", sep)

        output_file.write("\n")


if __name__ == '__main__':
    sys.exit(Hamming_Distance_Analysis(sys.argv))