def runEdgeR(outfile,
             outfile_prefix="edger.",
             fdr=0.1,
             prefix="",
             dispersion=None,
             ref_group=None,
             ref_regex=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, has_replicates, has_pairs = groupTagData(ref_group,
                                                            ref_regex)

    # output heatmap plot
    R.png('%(outfile_prefix)sheatmap.png' % locals())
    plotCorrelationHeatmap()
    R['dev.off']()

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

    if has_pairs:
        # output difference between groups
        R.png('''%(outfile_prefix)sbalance_groups.png''' % locals())
        first = True
        for g1, g2 in itertools.combinations(groups, 2):
            R('''a = rowSums( countsTable[groups == '%s'] ) ''' % g1)
            R('''b = rowSums( countsTable[groups == '%s'] ) ''' % g2)
            if first:
                R('''plot(cumsum(sort(a - b)), type = 'l') ''')
                first = False
            else:
                R('''lines(cumsum(sort(a - b))) ''')

        R['dev.off']()

        R('''suppressMessages(library('ggplot2'))''')
        R('''suppressMessages(library('reshape'))''')

        # output difference between pairs within groups
        first = True
        legend = []
        for pair in pairs:
            for g1, g2 in itertools.combinations(groups, 2):
                key = re.sub("-", "_", "pair_%s_%s_vs_%s" % (pair, g1, g2))
                legend.append(key)
                # print R('''colnames( countsTable) ''')
                # print R(''' pairs=='%s' ''' % pair)
                # print R(''' groups=='%s' ''' % g1)
                R('''a = rowSums( countsTable[pairs == '%s' & groups == '%s'])''' % (
                    pair, g1))
                R('''b = rowSums( countsTable[pairs == '%s' & groups == '%s'])''' % (
                    pair, g2))
                R('''c = cumsum( sort(a - b) )''')
                R('''c = c - min(c)''')
                if first:
                    data = R('''d = data.frame( %s = c)''' % key)
                    first = False
                else:
                    R('''d$%s = c''' % key)

        # remove row names (gene idenitifiers)
        R('''row.names(d) = NULL''')
        # add numbers of genes (x-axis)
        R('''d$genes=1:nrow(d)''')

        # merge data for ggplot
        R('''d = melt( d, 'genes', variable_name = 'comparison' )''')

        # plot
        R('''gp = ggplot(d)''')
        R('''pp = gp + \
            geom_line(aes(x=genes,y=value,group=comparison,color=comparison))''')

        try:
            R.ggsave('''%(outfile_prefix)sbalance_pairs.png''' % locals())
            R['dev.off']()
        except rpy2.rinterface.RRuntimeError, msg:
            E.warn("could not plot: %s" % msg)

    # build DGEList object
    # ignore message: "Calculating library sizes from column totals"
    R('''countsTable = suppressMessages(DGEList(countsTable, group=groups))''')

    # Relevel groups to make the results predictable - IMS
    if ref_group is not None:
        R('''countsTable$samples$group <- relevel(countsTable$samples$group,
        ref = "%s")''' % ref_group)
    else:
        # if no ref_group provided use first group in groups
        R('''countsTable$sample$group <- relevel(countsTable$samples$group,
        ref = "%s")''' % groups[0])

    # calculate normalisation factors
    E.info("calculating normalization factors")
    R('''countsTable = calcNormFactors( countsTable )''')
    E.info("output")

    # output MDS plot
    R.png('''%(outfile_prefix)smds.png''' % locals())
    try:
        R('''plotMDS( countsTable )''')
    except rpy2.rinterface.RRuntimeError, msg:
        E.warn("can not plot mds: %s" % msg)
    R['dev.off']()

    # build design matrix
    if has_pairs:
        R('''design = model.matrix(~pairs + countsTable$samples$group)''')
    else:
        R('''design = model.matrix(~countsTable$samples$group)''')

    # R('''rownames(design) = rownames( countsTable$samples )''')
    # R('''colnames(design)[length(colnames(design))] = "CD4" ''' )

    # fitting model to each tag
    if has_replicates:
        # estimate common dispersion
        R('''countsTable = estimateGLMCommonDisp(countsTable, design)''')
        # estimate tagwise dispersion
        R('''countsTable = estimateGLMTagwiseDisp(countsTable, design)''')
        # fitting model to each tag
        R('''fit = glmFit(countsTable, design)''')
    else:
        # fitting model to each tag
        if dispersion is None:
            raise ValueError("no replicates and no dispersion")
        E.warn("no replicates - using a fixed dispersion value")
        R('''fit = glmFit(countsTable, design, dispersion=%f)''' %
          dispersion)

    # perform LR test
    R('''lrt = glmLRT(fit)''')

    E.info("Generating output")

    # output cpm table
    R('''suppressMessages(library(reshape2))''')
    R('''countsTable.cpm <- cpm(countsTable, normalized.lib.sizes=TRUE)''')
    R('''melted <- melt(countsTable.cpm)''')
    R('''names(melted) <- c("test_id", "sample", "ncpm")''')
    # melt columns are factors - convert to string for sorting
    R('''melted$test_id = levels(melted$test_id)[as.numeric(melted$test_id)]''')
    R('''melted$sample = levels(melted$sample)[as.numeric(melted$sample)]''')
    # sort cpm table by test_id and sample
    R('''sorted = melted[with(melted, order(test_id, sample)),]''')
    R('''gz = gzfile("%(outfile_prefix)scpm.tsv.gz", "w" )''' % locals())
    R('''write.table(sorted, file=gz, sep = "\t",
                     row.names=FALSE, quote=FALSE)''')
    R('''close(gz)''')

    # compute adjusted P-Values
    R('''padj = p.adjust(lrt$table$PValue, 'BH')''')

    rtype = collections.namedtuple("rtype", "lfold logCPM LR pvalue")

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

    # I am assuming that logFC is the base 2 logarithm foldchange.
    # Parse results and parse to file
    results = []
    counts = E.Counter()

    df = R('''lrt$table''')
    df["padj"] = R('''padj''')

    counts.significant = sum(df.padj <= fdr)
    counts.insignificant = sum(df.padj > fdr)
    counts.significant_over = sum((df.padj <= fdr) & (df.logFC > 0))
    counts.significant_under = sum((df.padj <= fdr) & (df.logFC < 0))
    counts.input = len(df)
    counts.all_over = sum(df.logFC > 0)
    counts.all_under = sum(df.logFC < 0)
    counts.fail = sum(df.PValue.isnull())
    counts.ok = counts.input - counts.fail

    df["fold"] = df.logFC.pow(2.0)
    df["significant"] = df.padj <= fdr

    # TODO: use pandas throughout
    for interval, d in df.iterrows():
        # fold change is determined by the alphabetical order of the factors.
        # Is the following correct?
        results.append(GeneExpressionResult._make((
            interval,
            groups[1],
            d.logCPM,
            0,
            groups[0],
            d.logCPM,
            0,
            d.PValue,
            d.padj,
            d.logFC,
            d.fold,
            d.logFC,  # no transform of lfold
            str(int(d.significant)),
            "OK")))

    writeExpressionResults(outfile, results)

    outf = IOTools.openFile("%(outfile_prefix)ssummary.tsv" % locals(), "w")
    outf.write("category\tcounts\n%s\n" % counts.asTable())
    outf.close()

# needs to put into class
##
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/runEdgeR/original