def runDESeq(outfile,
             outfile_prefix="deseq.",
             fdr=0.1,
             prefix="",
             fit_type="parametric",
             dispersion_method="pooled",
             sharing_mode="maximum",
             ref_group=None,
             ref_regex=None,
             ):
    '''run DESeq on countsTable.

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

    The current analysis follows the analysis as outlined in version
    1.14.0

    DESeq ignores any pair information in the design matrix.

    The output is treatment and control. Fold change values are
    computed as treatment divided by control.

    '''

    # load library
    R('''suppressMessages(library('DESeq'))''')
    R('''suppressMessages(library('gplots'))''')
    R('''suppressMessages(library('RColorBrewer'))''')

    groups, pairs, has_replicates, has_pairs = groupTagData(ref_group,
                                                            ref_regex)

    # Run DESeq
    # Create Count data object
    E.info("running DESeq: replicates=%s" % (has_replicates))
    R('''cds <-newCountDataSet( countsTable, groups) ''')

    # Estimate size factors
    R('''cds <- estimateSizeFactors( cds )''')

    no_size_factors = R('''is.na(sum(sizeFactors(cds)))''')[0]
    if no_size_factors:
        E.warn("no size factors - can not estimate - no output")
        return

    # estimate variance
    if has_replicates:
        E.info("replicates - estimating variance from replicates")
    else:
        E.info("no replicates - estimating variance with method='blind'")
        dispersion_method = "blind"

    E.info("dispersion_method=%s, fit_type=%s, sharing_mode=%s" %
           (dispersion_method, fit_type, sharing_mode))
    R('''cds <- estimateDispersions( cds,
    method='%(dispersion_method)s',
    fitType='%(fit_type)s',
    sharingMode='%(sharing_mode)s')''' % locals())

    # bring into python namespace
    cds = R('''cds''')

    # plot fit - if method == "pooled":
    if dispersion_method == "pooled":
        R.png('%sdispersion_estimates_pooled.png' %
              outfile_prefix)
        R.plotDispEsts(cds)
        R['dev.off']()
    elif not has_replicates:
        # without replicates the following error appears
        # in the rpy2.py2ri conversion:
        #   'dims' cannot be of length 0
        pass
    else:
        dispersions = R('''ls(cds@fitInfo)''')
        for dispersion in dispersions:
            R.png('%sdispersion_estimates_%s.png' %
                  (outfile_prefix, dispersion))
        R.plotDispEsts(cds, name=dispersion)
        R['dev.off']()

    # plot size factors
    deseqPlotSizeFactors('%(outfile_prefix)ssize_factors.png' % locals())

    # output size factors
    deseqOutputSizeFactors("%(outfile_prefix)ssize_factors.tsv" % locals())

    # plot scatter plots of pairs
    deseqPlotPairs('%(outfile_prefix)spairs.png' % locals())

    if dispersion_method not in ("blind",):
        # also do a blind dispersion estimate for
        # a variance stabilizing transform
        R('''cds_blind <- estimateDispersions( cds,
        method='blind',
        fitType='%(fit_type)s',
        sharingMode='%(sharing_mode)s')''' % locals())
    else:
        R('''cds_blind = cds''')

    # perform variance stabilization for log2 fold changes
    vsd = R('''vsd = varianceStabilizingTransformation(cds_blind)''')
    # output normalized counts (in order)
    # gzfile does not work with rpy 2.4.2 in python namespace
    # using R.gzfile, so do it in R-space

    R('''t = counts(cds, normalized=TRUE);
    write.table(t[order(rownames(t)),],
    file=gzfile('%(outfile_prefix)scounts.tsv.gz'),
    row.names=TRUE,
    col.names=NA,
    quote=FALSE,
    sep='\t') ''' % locals())

    # output variance stabilized counts (in order)
    R('''t = exprs(vsd);
    write.table(t[order(rownames(t)),],
    file=gzfile('%(outfile_prefix)svsd.tsv.gz'),
    row.names=TRUE,
    col.names=NA,
    quote=FALSE,
    sep='\t')
    ''' % locals())

    # plot correlation heatmap of variance stabilized data
    deseqPlotCorrelationHeatmap(
        '%scorrelation_heatmap.png' % outfile_prefix,
        vsd)

    # plot PCA
    deseqPlotPCA('%spca.png' % outfile_prefix,
                 vsd)

    # plot gene heatmap for all genes - order by average expression
    # subtract one to get numpy indices
    select = R.order(R.rowMeans(R.counts(cds)), decreasing=True)
    # the following uses R-based indexing
    deseqPlotGeneHeatmap(
        '%sgene_heatmap.png' % outfile_prefix,
        R['as.matrix'](R.exprs(vsd).rx(select)))

    # plot heatmap of top 200 expressed genes
    deseqPlotGeneHeatmap(
        '%sgene_heatmap_top200.png' % outfile_prefix,
        R['as.matrix'](R.exprs(vsd).rx(select[:200])))

    # Call diffential expression for all pairings of groups included in the
    # design

    all_results = []
    for combination in itertools.combinations(groups, 2):
        control, treatment = combination
        gfix = "%s_vs_%s_" % (control, treatment)

        outfile_groups_prefix = outfile_prefix + gfix
        E.info(("calling differential expression for "
                "control=%s vs treatment=%s") %
               (control, treatment))
        res = R('''res = nbinomTest(cds, '%s', '%s')''' % (control, treatment))

        # plot significance
        R.png('''%(outfile_groups_prefix)ssignificance.png''' % locals())
        R('''plot(
        res$baseMean,
        res$log2FoldChange,
        log="x",
        pch=20, cex=.1,
        col = ifelse( res$padj < %(fdr)s, "red", "black"))''' % locals())
        R['dev.off']()

        # plot pvalues against rowsums
        deseqPlotPvaluesAgainstRowsums(
            '%(outfile_groups_prefix)spvalue_rowsums.png' % locals())

        E.info("Generating output (%s vs %s)" % (control, treatment))

        # get variance stabilized fold changes - note the reversal of
        # treatment/control
        R('''vsd_l2f =
        (rowMeans(exprs(vsd)[,conditions(cds) == '%s', drop=FALSE])
        - rowMeans( exprs(vsd)[,conditions(cds) == '%s', drop=FALSE]))''' %
          (treatment, control))

        # plot vsd correlation, see Figure 14 in the DESeq manual
        # if you also want to colour by expression level
        R.png('''%(outfile_groups_prefix)sfold_transformation.png''' %
              locals())
        R('''plot(
        res$log2FoldChange, vsd_l2f,
        pch=20, cex=.1,
        col = ifelse( res$padj < %(fdr)s, "red", "black" ) )''' % locals())
        R['dev.off']()

        # plot heatmap of differentially expressed genes
        # plot gene heatmap for all genes - order by average expression
        select = R('''select = res['padj'] < %f''' % fdr)

        if R('''sum(select)''')[0] > 0:
            E.info('%s vs %s: plotting %i genes in heatmap' %
                   (treatment, control, len(select)))
            data = R.exprs(vsd).rx(select)

            if not isinstance(data, rpy2.robjects.vectors.FloatVector):
                order = R.order(R.rowMeans(data), decreasing=True)
                deseqPlotGeneHeatmap(
                    '%sgene_heatmap.png' % outfile_groups_prefix,
                    R['as.matrix'](data[order]),
                    Colv=False,
                    Rowv=True)
            else:
                E.warn('can not plot differentially expressed genes')
        else:
            E.warn('no differentially expressed genes at fdr %f' % fdr)

        # Plot pvalue histogram
        R.png(locals() % '''%(outfile_groups_prefix)spvalue_histogram.png''')
        R('''pvalues = res$pval''')
        R('''hist(pvalues, breaks=50, col='skyblue' )''')
        R['dev.off']()

        # Plot diagnostic plots for FDR
        if has_replicates:
            R('''orderInPlot = order(pvalues)''')
            R('''showInPlot = (pvalues[orderInPlot] < 0.08)''')
            # Jethro - previously plotting x =
            # pvalues[orderInPlot][showInPlot]
            # pvalues[orderInPlot][showInPlot] contains all NA values
            # from pvalues which(showInPlot) doesn't... removing NA
            # values
            R('''true.pvalues <- pvalues[orderInPlot][showInPlot]''')
            R('''true.pvalues <- true.pvalues[is.finite(true.pvalues)]''')
            if R('''sum(showInPlot)''')[0] > 0:
                R.png('''%(outfile_groups_prefix)sfdr.png''' % locals())
                # failure when no replicates:
                # rpy2.rinterface.RRuntimeError:
                # Error in plot.window(...) : need finite 'xlim' values
                R('''plot( seq( along=which(showInPlot)),
                           true.pvalues,
                           pch='.',
                           xlab=expression(rank(p[i])),
                           ylab=expression(p[i]))''')
                R('''abline(a = 0, b = %(fdr)f / length(pvalues), col = "red")
                ''' % locals())
                R['dev.off']()
            else:
                E.warn('no p-values < 0.08')

        # Add log2 fold with variance stabilized l2fold value
        R('''res$transformed_log2FoldChange = vsd_l2f''')

        # Parse results and parse to file
        results, counts = deseqParseResults(control,
                                            treatment,
                                            fdr=fdr)

        all_results += results

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

    writeExpressionResults(outfile, all_results)
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/runDESeq/swapped operands of `%`