def runEdgeRPandas(counts,
                   design_table,
                   outfile,
                   outfile_prefix="edger.",
                   fdr=0.1,
                   prefix="",
                   dispersion=None,
                   ref_group=None):
    '''run EdgeR on countsTable.

    Results are stored in *outfile* and files prefixed by *outfile_prefix*.

    The dispersion is usually measuered from replicates. If there are no
    replicates, you need to set the *dispersion* explicitely.

    See page 13 of the EdgeR user guide::

       2. Simply pick a reasonable dispersion value, based on your
       experience with similar data, and use that. Although
       subjective, this is still more defensible than assuming Poisson
       variation. Typical values are dispersion=0.4 for human data,
       dispersion=0.1 for data on genetically identical model
       organisms or dispersion=0.01 for technical replicates.

    '''

    # load library
    R('''suppressMessages(library('edgeR'))''')

    (groups, pairs, conds, factors, has_replicates,
     has_pairs) = checkTagGroupsPandas(design_table, ref_group)

    if not has_replicates and dispersion is None:
        raise ValueError("no replicates and no dispersion")

    # output heatmap plot
    plotCorrelationHeatmapMatplot(counts,
                                  '%(outfile_prefix)sheatmap.png' % locals(),
                                  cor_method="spearman")

    E.info('running EdgeR: groups=%s, pairs=%s, replicates=%s, pairs=%s' %
           (groups, pairs, has_replicates, has_pairs))

    r_counts = pandas2ri.py2ri(counts)

    passDFtoRGlobalEnvironment = R('''function(df){
    countsTable <<- df}''')
    passDFtoRGlobalEnvironment(r_counts)

    if has_pairs:
        # output difference between groups
        # TS #####
        # this is performed on non-normalised data
        # should we use Counts.py to normalise first?
        # also, this isn't edgeR specific, should this be
        # moved to a seperate summary function?
        # also move the MDS plotting?
        # #####
        first = True
        pairs_df = pandas.DataFrame()
        nrows = len(counts.index)
        n = 0
        for g1, g2 in itertools.combinations(groups, 2):

            keep_a = [x == g1 for x in conds]
            counts_a = counts.iloc[:, keep_a]
            keep_b = [x == g2 for x in conds]
            counts_b = counts.iloc[:, keep_b]
            index = range(n, n+nrows)
            n += nrows
            a = counts_a.sum(axis=1)
            b = counts_b.sum(axis=1)
            diff = a-b
            diff.sort()
            temp_df = pandas.DataFrame({"cumsum": np.cumsum(diff).tolist(),
                                        "comb": "_vs_".join([g1, g2]),
                                        "id": range(0, nrows)},
                                       index=index)
            pairs_df = pairs_df.append(temp_df)
        plot_pairs = R('''function(df, outfile){
        suppressMessages(library('ggplot2'))
        p = ggplot(df, aes(y=cumsum, x=id)) +
        geom_line(aes(col=factor(comb))) +
        scale_color_discrete(name="Comparison") +
        xlab("index") + ylab("Cumulative sum")
        ggsave("%(outfile_prefix)sbalance_groups.png", plot = p)}
        ''' % locals())

        r_pairs_df = pandas2ri.py2ri(pairs_df)
        plot_pairs(r_pairs_df)

        # output difference between pairs within groups
        first = True
        legend = []
        n = 0
        pairs_in_groups_df = pandas.DataFrame()
        for pair in set(pairs):
            for g1, g2 in itertools.combinations(groups, 2):

                key = "pair-%s-%s-vs-%s" % (pair, g1, g2)
                legend.append(key)
                keep_a = [x == pair and y == g1 for x, y in zip(pairs, conds)]
                keep_b = [x == pair and y == g2 for x, y in zip(pairs, conds)]

                # check pair and group combination actually exists
                if sum(keep_a) > 0 and 0 > sum(keep_b):
                    counts_a = counts.iloc[:, keep_a]
                    counts_b = counts.iloc[:, keep_b]
                    index = range(n, n+nrows)
                    n += nrows
                    a = counts_a.sum(axis=1)
                    b = counts_b.sum(axis=1)
                    diff = a-b
                    diff.sort()
                    comparison = "pair-%s-%s-vs-%s" % (pair, g1, g2)
                    temp_df = pandas.DataFrame({"cumsum": np.cumsum(diff).tolist(),
                                                "comb": comparison,
                                                "id": range(0, nrows)},
                                               index=index)
                    pairs_in_groups_df = pairs_in_groups_df.append(temp_df)

        plot_pairs_in_groups = R('''function(df, outfile){
        suppressMessages(library('ggplot2'))
        p = ggplot(df, aes(y=cumsum, x=id)) +
        geom_line(aes(col=factor(comb))) +
        scale_color_discrete(name="Comparison") +
        xlab("index") + ylab("Cumulative sum")
        ggsave("%(outfile_prefix)sbalance_pairs.png", plot = p)}
        ''' % locals())

        r_pairs_in_groups_df = pandas2ri.py2ri(pairs_in_groups_df)
        plot_pairs_in_groups(r_pairs_in_groups_df)

    # create r objects
    r_counts = pandas2ri.py2ri(counts)
    r_groups = ro.StrVector(conds)
    r_pairs = ro.StrVector(pairs)
    r_has_pairs = ro.default_py2ri(has_pairs)
    r_has_replicates = ro.default_py2ri(has_replicates)

    if dispersion is not None:
        r_dispersion = ro.default_py2ri(dispersion)
    else:
        r_dispersion = ro.default_py2ri(False)

    if ref_group is not None:
        r_ref_group = ro.default_py2ri(ref_group)
    else:
        r_ref_group = ro.default_py2ri(groups[0])

    # build DGEList object
    buildDGEList = R('''
    suppressMessages(library('edgeR'))
    function(counts, groups, ref_group){
    countsTable = DGEList(counts, group=groups)
    countsTable$samples$group <- relevel(countsTable$samples$group,
    ref = ref_group)
    countsTable = calcNormFactors(countsTable)
    return(countsTable)
    }''')

    r_countsTable = buildDGEList(r_counts, r_groups, r_ref_group)

    # output MDS plot
    # TT - should this be performed on the normalised counts table?(see above)
    R.png('''%(outfile_prefix)smds.png''' % locals())
    try:
        MDSplot = R('''function(counts){
        plotMDS(counts)}''')
        MDSplot(r_counts)
    except rpy2.rinterface.RRuntimeError, msg:
        E.warn("can not plot mds: %s" % msg)
    R['dev.off']()

    # build design matrix
    buildDesign = R('''function(countsTable, has_pairs, pairs){
    if (has_pairs==TRUE) {
        design <- model.matrix( ~pairs + countsTable$samples$group ) }
    else {
        design <- model.matrix( ~countsTable$samples$group ) }
    return(design)
    }''')

    r_design = buildDesign(r_countsTable, r_has_pairs, r_pairs)

    fitModel = R('''function(countsTable, design, has_replicates, dispersion){
    if (has_replicates == TRUE) {
        # estimate common dispersion
        countsTable = estimateGLMCommonDisp( countsTable, design )
        # estimate tagwise dispersion
        countsTable = estimateGLMTagwiseDisp( countsTable, design )
        # fitting model to each tag
        fit = glmFit( countsTable, design ) }
    else {
        # fitting model to each tag
        fit = glmFit(countsTable, design, dispersion=dispersion) }
    return(fit)
    }''')

    r_fit = fitModel(r_countsTable, r_design, r_has_replicates, r_dispersion)

    E.info("Generating output")

    # perform LR test
    lrtTest = R('''function(fit, prefix){
    lrt = glmLRT(fit)
    # save image for access to the whole of the lrt object
    save.image(paste0(prefix,"lrt.RData"))
    return(lrt)
    }''')
    r_lrt = lrtTest(r_fit, outfile_prefix)

    # return statistics table - must be a better way to do this?
    extractTable = R('''function(lrt){ return(lrt$table)}''')
    r_lrt_table = extractTable(r_lrt)

    # output cpm table
    outputCPMTable = R('''function(countsTable, outfile_prefix, lrt){
    suppressMessages(library(reshape2))
    countsTable.cpm <- cpm(countsTable, normalized.lib.sizes=TRUE)
    melted <- melt(countsTable.cpm)
    names(melted) <- c("test_id", "sample", "ncpm")

    # melt columns are factors - convert to string for sorting
    melted$test_id = levels(melted$test_id)[as.numeric(melted$test_id)]
    melted$sample = levels(melted$sample)[as.numeric(melted$sample)]

    # sort cpm table by test_id and sample
    sorted <- melted[with(melted, order(test_id, sample)),]
    gz <- gzfile(paste0(outfile_prefix,"cpm.tsv.gz"), "w" )
    write.table(sorted, file=gz, sep = "\t", row.names=FALSE, quote=FALSE)
    close(gz)}''')

    outputCPMTable(r_countsTable, outfile_prefix, r_lrt)

    # output differences between pairs
    if len(groups) == 2:
        plotSmear = R('''function(countsTable, outfile_prefix){
        png(paste0(outfile_prefix,"smaplot.png"
        plotSmear(countsTable, pair=c('%s'))
        abline(h=c(-2, 2), col='dodgerblue')
        dev.off()}''' % "','".join(groups))

    lrt_table = pandas2ri.ri2py(r_lrt_table)

    n_rows = lrt_table.shape[0]
    df_dict = makeEmptyDataFrameDict()

    # import r stats module for BH adjustment
    stats = importr('stats')

    df_dict["control_name"].extend((groups[0],)*n_rows)
    df_dict["treatment_name"].extend((groups[1],)*n_rows)
    df_dict["test_id"].extend(lrt_table.index)
    df_dict["control_mean"].extend(lrt_table['logCPM'])
    df_dict["treatment_mean"].extend(lrt_table['logCPM'])
    df_dict["control_std"].extend((0,)*n_rows)
    df_dict["treatment_std"].extend((0,)*n_rows)
    df_dict["p_value"].extend(lrt_table['PValue'])
    df_dict["p_value_adj"].extend(
        list(stats.p_adjust(FloatVector(df_dict["p_value"]), method='BH')))
    df_dict["significant"].extend(
        [int(float(x) < fdr) for x in df_dict["p_value_adj"]])
    df_dict["l2fold"].extend(list(numpy.log2(lrt_table['logFC'])))

    # TS -note: the transformed log2 fold change is not transformed!
    df_dict["transformed_l2fold"].extend(list(numpy.log2(lrt_table['logFC'])))

    # TS -note: check what happens when no fold change is available
    # TS -may need an if/else in list comprehension. Raise E.warn too?
    df_dict["fold"].extend([math.pow(2, float(x)) for x in lrt_table['logFC']])

    # set all status values to "OK"
    # TS - again, may need an if/else, check...
    df_dict["status"].extend(("OK",)*n_rows)

    results = pandas.DataFrame(df_dict)
    results.set_index("test_id", inplace=True)
    results.to_csv(outfile, sep="\t", header=True, index=True)

    counts = E.Counter()
    counts.signficant = sum(results['significant'])
    counts.insignficant = (len(results['significant']) - counts.signficant)
    counts.all_over = sum([x > 0 for x in results['l2fold']])
    counts.all_under = sum([x < 0 for x in results['l2fold']])
    counts.signficant_over = sum([results['significant'][x] == 1 and
                                 results['l2fold'][x] > 1 for
                                 x in range(0, n_rows)])
    counts.signficant_under = sum([results['significant'][x] == 1 and
                                   results['l2fold'][x] < 1 for
                                   x in range(0, n_rows)])

    outf = IOTools.openFile("%(outfile_prefix)ssummary.tsv" % locals(), "w")
    outf.write("category\tcounts\n%s\n" % counts.asTable())
    outf.close()
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/runEdgeRPandas/swapped operands of `>`