def outputSpikeIns(filename_tags,
                   outfile,
                   output_filename_pattern,
                   filename_design=None,
                   foldchange_max=10.0,
                   expression_max=5.0,
                   max_counts_per_bin=100,
                   expression_bin_width=0.5,
                   foldchange_bin_width=0.5,
                   iterations=1):

    E.info("loading tag data from %s" % filename_tags)

    if filename_design is not None:
        # load all tag data
        counts_table, design_table = loadTagDataPandas(
            filename_tags, filename_design)

        # filter
        counts_table = filterTagDataPandas(counts_table, design_table)

    else:
        raise NotImplementedError("reading full table not implemented")

    nobservations, nsamples = counts_table.shape

    groups = list(set(design_table["group"]))
    if len(groups) != 2:
        raise NotImplementedError(
            "spike in only implemented for one pairwise comparison")

    # select group data
    group1 = counts_table[design_table[design_table.group == groups[0]].index]
    group2 = counts_table[design_table[design_table.group == groups[1]].index]

    outfile.write("interval_id\t%s\t%s\n" %
                  ("\t".join(group1.columns), "\t".join(group2.columns)))
    outf_info = IOTools.openFile(output_filename_pattern + "info.tsv.gz", "w")
    outf_info.write("interval_id\tl10_expression\tl2fold\n")

    # important: convert to matrixes, otherwise there will be a per-name lookup
    # when l10average or l2fold are computed.
    group1 = group1.as_matrix()
    group2 = group2.as_matrix()

    # compute bins
    expression_bins = numpy.arange(0, expression_max, expression_bin_width)
    fold_change_bins = numpy.arange(-foldchange_max,
                                    foldchange_max, foldchange_bin_width)

    E.info("l10expression bins=%s" % expression_bins)
    E.info("l2fold change bins=%s" % fold_change_bins)

    # output values
    output_counts = numpy.zeros(
        (len(expression_bins) + 1, len(fold_change_bins) + 1))

    for iteration in range(iterations):

        # randomize order
        group1 = numpy.random.permutation(group1)
        group2 = numpy.random.permutation(group2)

        # compute means and foldchanges
        group1_mean = group1.mean(1)
        group2_mean = group2.mean(1)

        # compute average expression by means
        l10average = numpy.log((group1_mean + group2_mean) / 2.0)

        # compute a fold change with pseudo count of 1
        l2fold = numpy.log2((group1_mean + 1.0) / (group2_mean + 1.0))

        # digitize arrays with bins
        l10average_idx = numpy.digitize(l10average, expression_bins)
        l2fold_idx = numpy.digitize(l2fold, fold_change_bins)

        interval_id = 0
        for idx, coord in enumerate(zip(l10average_idx, l2fold_idx)):
            # assert expression_bins[coord[0]-1] <= l10average[idx] < expression_bins[coord[0]]
            # assert fold_change_bins[coord[1]-1] <= l2fold[idx] <  expression_bins[coord[1]]

            if output_counts[coord] >= max_counts_per_bin:
                continue

            output_counts[coord] += 1
            outf_info.write("spike%i\t%f\t%f\n" % (interval_id,
                                                   l10average[idx],
                                                   l2fold[idx]))

            outfile.write("spike%i\t%s\t%s\n" %
                          (interval_id,
                           "\t".join(
                               map(str, list(group1[idx]))),
                           "\t".join(map(str, list(group2[idx])))))
            interval_id += 1

    outf_info.close()

    df = pandas.DataFrame(output_counts)
    df.index = list(expression_bins) + [expression_max]
    df.columns = list(fold_change_bins) + [foldchange_max]

    df.to_csv(IOTools.openFile(output_filename_pattern + "counts.tsv.gz", "w"),
              sep="\t")

    E.info("output %i spike in intervals" % interval_id)
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/outputSpikeIns/original