def run(self,
            counts,
            design,
            model=None,
            contrasts=None,
            outfile_prefix=None,
            fdr=0.1):

        counts.table = counts.table[design.table.index]

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

        if design.factors is not None:
            r_factors_df = pandas2ri.py2ri(design.factors)
        else:
            r_factors_df = ro.default_py2ri(False)

        if contrasts is not None:
            DEtype = "GLM"

            # if model not included, use the column names from design.factors
            if not model:

                if design.factors is not None:
                    model = "~" + "+".join(design.factors.columns)
                    model_terms = design.factors.columns.values.tolist()

                else:
                    E.warn("need to supply a full model or else "
                           "additional columns in the design table "
                           "which will be taken as the full model")
            else:
                model_terms = [x for x in re.split("[\+~ ]+", model)[1:]
                               if x != "0"]

            r_model = ro.default_py2ri(model)

        else:
            DEtype = "pairwise"

        r_DEtype = ro.default_py2ri(DEtype)

        E.info('running DESeq2: groups=%s, pairs=%s, replicates=%s, pairs=%s, '
               'additional_factors:' %
               (design.groups, design.pairs, design.has_replicates,
                design.has_pairs))
        E.info(design.factors)

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

        # build design matrix
        if DEtype == "pairwise":
            buildDesign = R('''
              function(counts, groups){

                design = data.frame(row.names = colnames(counts),
                                    condition = groups)
                return(design)}''')

            r_design = buildDesign(r_counts, r_groups)

            buildCountDataSet = R('''
            function(counts, design){

            for(column in colnames(design)){
              design[[column]] = factor(design[[column]])
            }

            dds <- suppressMessages(DESeqDataSetFromMatrix(
                     countData= counts,
                     colData = design,
                     design = ~condition))

            return(dds)

            }''' % locals())

            r_dds = buildCountDataSet(r_counts, r_design)

            performDifferentialTesting = R('''
            function(dds){
            dds = suppressMessages(DESeq(dds))

            png("%(outfile_prefix)s_dispersion.png")
            plotDispEsts(dds)
            dev.off()

            res = suppressMessages(results(dds, addMLE=TRUE))
            png(paste0(c("%(outfile_prefix)s", "MA.png"), collapse="_"))
            plotMA(res, alpha=%(fdr)s)
            dev.off()

            res = as.data.frame(res)

            c = counts(dds, normalized = TRUE)

            contrast = "condition"
            res$contrast = contrast
            contrast_levels = levels(dds@colData[[contrast]])

            if(length(contrast_levels)==2){
              res$control = contrast_levels[1]
              res$treatment = contrast_levels[2]}
            else{
              res$control = contrast
              res$treatment = contrast}

            res['test_id'] = rownames(res)

            return(res)

            }''' % locals())

            results = pandas2ri.ri2py(performDifferentialTesting(r_dds))

        # DEtype == "GLM"
        else:
            r_design = pandas2ri.py2ri(design.table)

            buildCountDataSet = R('''
            function(counts, design, model){

            for(column in colnames(design)){
              design[[column]] = factor(design[[column]])
            }

            full_model <- formula("%(model)s")

            dds <- suppressMessages(DESeqDataSetFromMatrix(
                     countData= counts,
                     colData = design,
                     design = full_model))

            return(dds)
            }''' % locals())

            r_dds = buildCountDataSet(r_counts, r_design, r_model)
            results = pandas.DataFrame()

            n = 0

            for contrast in contrasts:
                assert contrast in design.table.columns, (
                    "contrast: %s not found in design factors columns" % contrast)
                model = [x for x in model_terms if x != contrast]
                model = "~" + "+".join(model)

                performDifferentialTesting = R('''
                function(dds){

                ddsLRT <- suppressMessages(
                DESeq(dds, test="LRT", reduced=formula(%(model)s),
                      betaPrior=TRUE))
                #DESeq(dds, betaPrior=TRUE))

                png("%(outfile_prefix)s_dispersion.png")
                plotDispEsts(ddsLRT)
                dev.off()

                contrast_levels = as.vector(levels(dds@colData$%(contrast)s))

                for(levels in combn(contrast_levels, 2, simplify=F)){


                    res = suppressMessages(results(ddsLRT, addMLE=TRUE,
                                  contrast=c("%(contrast)s",
                                  levels[1], levels[2])))

                    png(paste0(c("%(outfile_prefix)s", "%(contrast)s",
                               levels[1], levels[2], "MA.png"), collapse="_"))
                    plotMA(res, alpha=%(fdr)s)
                    dev.off()
                }

                res = as.data.frame(res)
                res$contrast = "%(contrast)s"

                if(length(contrast_levels)==2){
                  tmp_df = data.frame(contrast_levels)
                  res$control = tmp_df[1,1]
                  res$treatment = tmp_df[2,1]
                  }
                else{
                  res$control = "%(contrast)s"
                  res$treatment = "%(contrast)s"}

                res['test_id'] = rownames(res)
                return(res)
                }''' % locals())

                tmp_results = pandas2ri.ri2py(performDifferentialTesting(r_dds))
                # tmp_results['test_id'] = tmp_results.index

                # need to set index to sequence of ints to avoid duplications
                n2 = n+tmp_results.shape[0]
                tmp_results.index = range(n, n2)
                n = n2

                results = results.append(tmp_results)

        final_result = DEResult_DESeq2(testTable=results)

        return final_result
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/DEExperiment_DESeq2.run/original