def runDESeq2(outfile,
              outfile_prefix="deseq2",
              fdr=0.1,
              ref_group=None,
              model=None,
              contrasts=None,
              plot=1,
              ):
    """
    Run DESeq2 on counts table.

    If no model is passed, then defaults to the group column in design file

    Does not make use of group tag data bc function doesn't accomodate
    multi-factor designs

    To Do: Parse results into standard output format.
    KB: I have done this but I'm not sure if it is compatible with complex
    design tables
    Fix fact that plotMA is hardcoded.
    """

    # load libraries
    R('''suppressMessages(library('DESeq2'))''')

    # Create metadata... this will eventually be a pandas dataframe
    if isinstance(R('''factors'''), rpy2.robjects.vectors.DataFrame):
        E.info("DESeq2: Merging additional factors in design file to"
               "create metadata table")
        R('''mdata <- cbind(groups, factors)''')
        mdata = tuple(R('''names(mdata)'''))
    else:
        R('''mdata <- data.frame(group=groups)''')
        mdata = "group"
    E.info(mdata % "DESeq2 colData headers are: %s")

    # Check for model and that model terms are in metadata table
    if model:
        assert contrasts, "Must specifiy contrasts is model design provided"
        terms = set([x for x in re.split("\W", model) if x != ''])
        assert terms.issubset(mdata), \
            "DESeq2: design formula has terms not present in colData"
    else:
        if mdata != "group":
            E.warn("DESeq2 model specified, with no metadata in design file")
        terms = ["group", ]
        model = "~ group"
    E.info("DESeq2 design formula is: %s" % model)

    # Create DESeqDataSet, using countsTable, mdata, model
    R('''suppressMessages(dds <- DESeqDataSetFromMatrix(countData=countsTable,
                                       colData=mdata,
                                       design=%(model)s))''' % locals())
    # WARNING: This is not done automatically... I don't know why?
    R('''colnames(dds) <- colnames(countsTable)''')
    E.info("Combined colata, design formula and counts table to create"
           " DESeqDataSet instance")

    # Rlog transform
    R('''suppressMessages(rld <- rlog(dds))''')

    if plot == 1:
        # Plot PCA of rlog transformed count data for top 500
        for factor in terms:
            outf = outfile_prefix + factor + "_PCAplot500.tiff"
            E.info("Creating PCA plot for factor: %s" % outf)
            R('''x <- plotPCA(rld, intgroup="%(factor)s")''' % locals())
            # R('''saveRDS(x, '%(outf)s')''' % locals())
            R('''tiff("%(outf)s")''' % locals())
            R('''plot(x)''')
            R('''dev.off()''')

    # Extract rlog transformed count data...

    R('''rlogtab = as.data.frame(assay(rld))''')
    R('''rlogtab$test_id = rownames(rlogtab)''')
    R('''rlogtab = rlogtab[, c(ncol(rlogtab), 1:ncol(rlogtab)-1)]''')
    R('''rlogtab = as.data.frame(rlogtab)''')
    R.data('rlogtab')
    rlog_out = R('rlogtab')
    rlogoutf = outfile_prefix + "rlog.tsv"

    rlog_out.to_csv(rlogoutf, sep="\t", index=False)
    os.system("gzip %s" % rlogoutf)

    # Run DESeq2
    R('''suppressMessages(dds <- DESeq(dds))''')
    E.info("Completed DESeq2 differential expression analysis")

    # Extract contrasts...
    if contrasts:
        contrasts = (x.split(":") for x in contrasts.split(","))
    else:
        # created by loadTagData...
        groups = R('''levels(groups)''')
        contrasts = (("group",) + x for x in itertools.combinations(groups, 2))

    df_final = pandas.DataFrame()
    raw_final = pandas.DataFrame()
    all_results = []
    for combination in contrasts:
        variable, control, treatment = combination

        # Fetch results
        gfix = "%s_%s_vs_%s" % (variable, control, treatment)
        outfile_groups_prefix = outfile_prefix + gfix + "_MAplot.png"
        R('''res <- results(dds, contrast=c("%(variable)s",
                                            "%(treatment)s",
                                            "%(control)s"))''' % locals())
        E.info("Extracting contrast for levels %s (treatment) vs %s (control)"
               " for factor %s" % (treatment, control, variable))

        # plot MA plot
        if plot == 1:
            R('''png("%(outfile_groups_prefix)s")''' % locals())
            R('''plotMA(res, alpha=%f)''' % fdr)
            R('''dev.off()''')
            E.info("Plotted MA plot for levels %s (treatment) vs %s (control)"
                   " for factor %s" % (treatment, control, variable))
        R('''res_df <- as.data.frame(res)''')
        R('''res_df$test_id = rownames(res_df)''')
        R('''res_df = res_df[, c(ncol(res_df), 1:ncol(res_df)-1)]''')
        R.data('res_df')
        raw_out = R('res_df')

        #  manipulate output into standard format
        df_out = deseq2ParseResults(fdr, vsd=False)

        #  label the deseq2 raw output file and append it to the raw output tab
        raw_out["treatment"] = [treatment, ]*len(df_out.index)
        raw_out["control"] = [control, ]*len(df_out.index)
        raw_out["variable"] = [variable, ]*len(df_out.index)
        raw_final = raw_final.append(raw_out, ignore_index=True)

        #  write the standardised output table
        df_out.to_csv(IOTools.openFile(outfile_prefix + gfix + ".tsv.gz", "w"),
                      sep="\t",
                      index_label="test_id")

        E.info("Extracted results table for contrast  '%s' (treatment) vs '%s'"
               " (control) for factor '%s'" % (treatment, control, variable))

        # append to final dataframe
        df_out.reset_index(inplace=True)
        df_out.rename(columns={"index": "test_id"}, inplace=True)
        df_final = df_final.append(df_out, ignore_index=True)

    results = df_final.values.tolist()

    # write final dataframe into standard format
    writeExpressionResults(outfile, results)

    rawoutf = outfile_prefix + "raw.tsv"
    raw_final.to_csv(rawoutf, sep="\t", index=False)
    os.system("gzip %s" % rawoutf)
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/runDESeq2/swapped operands of `%`