def __call__(self, probesets,
                 treatments,
                 controls,
                 pattern=None,
                 fdr=0.10,
                 ngenes=None,
                 npermutations=1000,
                 ndelta=10,
                 method="ttest",
                 use_excel_sam=False,
                 treatment_label="treatment",
                 control_label="control"):

        if ngenes and fdr:
            raise ValueError("either supply ngenes or fdr, but not both.")

        R.library("siggenes")

        m = numpy.matrix(treatments + controls)
        m = numpy.transpose(m)

        E.debug("build expression matrix: %i x %i" % m.shape)

        labels = numpy.array([1] * len(treatments) + [0] * len(controls))
        # 1000 permutations for P-Values of down to 0.0001. Setting this
        # to a high value improved reproducibility of results.

        kwargs = {}
        # kwargs set to replicate excel SAM
        if use_excel_sam:
            kwargs.update(
                {"control":
                 R('''samControl( lambda = 0.5, n.delta = %(ndelta)s) ''' %
                   locals()),
                 "med": True,
                 "var.equal": True})
        else:
            kwargs.update({"control":
                           R('''samControl( n.delta = %(ndelta)s ) ''' %
                             locals())},)

        # the option B needs to be not set if wilc.stat is chosen

        if method == "ttest":
            kwargs["method"] = R('''d.stat''')
            kwargs["B"] = npermutations
        elif method == "wilc":
            kwargs["method"] = R('''wilc.stat''')
        elif method == "cat":
            kwargs["method"] = R('''cat.stat''')
        else:
            raise ValueError(method % "unknown statistic `%s`")

        E.info("running sam with the following options: %s" % str(kwargs))

        a = R.sam(numpy.array(m),
                  labels,
                  gene_names=numpy.array(probesets),
                  **kwargs)

        # E.debug("%s" % str(a))

        R.assign("a", a)

        fdr_data = collections.namedtuple("sam_fdr", (
            "delta", "p0", "false", "significant", "fdr", "cutlow",
            "cutup", "j2", "j1"))
        cutoff_data = collections.namedtuple(
            "sam_cutoff", ("delta", "significant", "fdr"))
        gene_data = collections.namedtuple(
            "sam_fdr", ("row", "dvalue", "stddev", "rawp", "qvalue", "rfold"))

        def _totable(robj):
            '''convert robj to a row-wise table.'''
            s = numpy.matrix(robj)
            t = [numpy.array(x).reshape(-1,) for x in s]
            return t

        # extract the fdr values
        # returns R matrix
        t = _totable(a.do_slot('mat.fdr'))
        assert len(t[0]) == len(fdr_data._fields)
        # for x in t: E.debug( "x=%s" % str(x))
        fdr_values = [fdr_data(*x) for x in t]

        # find d cutoff
        if fdr is not None and fdr > 0:
            s = numpy.matrix(R.findDelta(a, fdr))
            try:
                cutoffs = [cutoff_data(*numpy.array(x).reshape(-1,))
                           for x in s]
                E.debug("sam cutoffs for fdr %f: %s" % (fdr, str(cutoffs)))
                cutoff = cutoffs[-1]
            except TypeError:
                E.debug("could not get cutoff")
                cutoff = None
        elif ngenes:
            s = numpy.matrix(R.findDelta(a, ngenes))
            try:
                cutoffs = [cutoff_data(*numpy.array(x).reshape(-1,))
                           for x in s]
                E.debug("sam cutoffs for fdr %f: %s" % (fdr, str(cutoffs)))
                cutoff = cutoffs[-1]
            except TypeError:
                E.debug("could not get cutoff")
                cutoff = None
        else:
            raise ValueError("either supply ngenes or fdr")

        # collect (unadjusted) p-values and qvalues for all probesets
        pvalues = dict(zip(probesets, R('''a@p.value''')))
        qvalues = dict(zip(probesets, R('''a@q.value''')))

        if pattern:
            outfile = pattern % "sam.pdf"
            R.pdf(outfile)
            if cutoff:
                R.plot(a, cutoff.delta)
            else:
                R.plot(a)
            R['dev.off']()

        siggenes = {}
        significant_genes = set()
        if cutoff is not None:
            E.debug("using cutoff %s" % str(cutoff))

            summary = R('''summary( a, %f )''' % cutoff.delta)

            # summary = R.summary( a, cutoff.delta )
            R.assign("summary", summary)

            significant_genes = set(
                [probesets[int(x) - 1] for x
                 in R('''summary@row.sig.genes''')])
            # E.debug( "significant genes=%s" % str(significant_genes))

            r_result = zip(*_totable(summary.do_slot('mat.sig')))

            if len(r_result) > 0:

                assert len(r_result[0]) == 6, \
                    "expected six columns from siggenes module, got: %s" % \
                    len(r_result[0])

                for x in r_result:
                    if x[4] > fdr:
                        E.warn("%s has qvalue (%f) larger than cutoff, but "
                               "is significant significant." % (str(x), x[4]))

                # except TypeError:
                # only a single value
                #     x = [r_result[y] for y in ("Row", "d.value", "stdev", "rawp", "q.value", "R.fold") ]
                #     if x[4] > fdr:
                #         E.warn( "%s has qvalue (%f) larger than cutoff, but is called significant." % (str(x), x[4]))

                siggenes[probesets[int(x[0]) - 1]] = gene_data(*x)

        else:
            E.debug("no cutoff found - no significant genes.")

        genes = []
        for probeset, treatment, control in zip(
                probesets, zip(*treatments), zip(*controls)):

            mean1, mean2 = numpy.mean(treatment), numpy.mean(control)

            if probeset in siggenes:
                s = siggenes[probeset]
                pvalue = s.rawp
                qvalue = s.qvalue
            else:
                pvalue = pvalues[probeset]
                qvalue = qvalues[probeset]

            significant = (0, 1)[probeset in significant_genes]

            genes.append(GeneExpressionResult._make((probeset,
                                                     treatment_label,
                                                     mean1,
                                                     numpy.std(treatment),
                                                     control_label,
                                                     mean2,
                                                     numpy.std(control),
                                                     pvalue,
                                                     qvalue,
                                                     mean1 - mean2,
                                                     math.pow(
                                                         2, mean1 - mean2),
                                                     math.pow(
                                                         2, mean1 - mean2),
                                                     significant,
                                                     "OK")))

        return genes, cutoff, fdr_values


#########################################################################
#########################################################################
#########################################################################
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Expression.py/SAM.__call__/swapped operands of `%`