code: def bh_fdr(pval):
      
      pval_array = np.array(pval)
      sorted_order = np.argsort(pval_array)
      original_order = np.argsort(sorted_order)
      pval_array = pval_array[sorted_order]
  
      # calculate the needed alpha
      n = float(len(pval))
      pval_adj = np.zeros(int(n))
      i = np.arange(1, int(n)+1, dtype=float)[::-1]  # largest to smallest
      pval_adj = np.minimum(1, cummin(n/i * pval_array[::-1]))[::-1]
      return pval_adj[original_order]
docs: A python implementation of the Benjamani-Hochberg FDR method.
  
      This code should always give precisely the same answer as using
      p.adjust(pval, method="BH") in R.
  
      Parameters
      ----------
      pval : list or array
          list/array of p-values
  
      Returns
      -------
      pval_adj : np.array
          adjusted p-values according the benjamani-hochberg method