python class for converting p-values to adjusted p-values (or q-values) for multiple comparison correction.

MCP_simulation.ipynb

MCPconverter.py

import pandas as pd
import numpy as np
import statsmodels as sms

class MCPConverter(object):
    """
    input: array of p-values.
    * convert p-value into adjusted p-value (or q-value)
    """
    def __init__(self, pvals, zscores=None):
        self.pvals = pvals
        self.zscores = zscores
        self.len = len(pvals)
        if zscores is not None:
            srted = np.array(sorted(zip(pvals.copy(), zscores.copy())))
            self.sorted_pvals = srted[:, 0]
            self.sorted_zscores = srted[:, 1]
        else:
            self.sorted_pvals = np.array(sorted(pvals.copy()))
        self.order = sorted(range(len(pvals)), key=lambda x: pvals[x])
    
    def adjust(self, method="holm"):
        """
        methods = ["bonferroni", "holm", "bh", "lfdr"]
         (local FDR method needs 'statsmodels' package)
        """
        if method is "bonferroni":
            return [np.min([1, i]) for i in self.sorted_pvals * self.len]
        elif method is "holm":
            return [np.min([1, i]) for i in (self.sorted_pvals * (self.len - np.arange(1, self.len+1) + 1))]
        elif method is "bh":
            p_times_m_i = self.sorted_pvals * self.len / np.arange(1, self.len+1)
            return [np.min([p, p_times_m_i[i+1]]) if i < self.len-1 else p for i, p in enumerate(p_times_m_i)]
        elif method is "lfdr":
            if self.zscores is None:
                raise ValueError("Z-scores were not provided.")
            return sms.stats.multitest.local_fdr(abs(self.sorted_zscores))
        else:
            raise ValueError("invalid method entered: '{}'".format(method))
            
    def adjust_many(self, methods=["bonferroni", "holm", "bh", "lfdr"]):
        if self.zscores is not None:
            df = pd.DataFrame(np.c_[self.sorted_pvals, self.sorted_zscores], columns=["p_values", "z_scores"])
            for method in methods:
                df[method] = self.adjust(method)
        else:
            df = pd.DataFrame(self.sorted_pvals, columns=["p_values"])
            for method in methods:
                if method is not "lfdr":
                    df[method] = self.adjust(method)
        return df