docs: def qualNorm_filter_genes(data, qualitative, pval_threshold=0.05, smoothing=1e-5, eps=1e-5):
      """
      Does qualNorm but returns a filtered gene set, based on a p-value threshold.
      """
      genes, cells = data.shape
      clusters = qualitative.shape[1]
      output = np.zeros((genes, clusters))
      missing_indices = []
      genes_included = []
      qual_indices = []
      thresholds = qualitative.min(1) + (qualitative.max(1) - qualitative.min(1))/2.0
      pvals = np.zeros(genes)
      for i in range(genes):
          if qualitative[i,:].max() == -1 and qualitative[i,:].min() == -1:
              missing_indices.append(i)
              continue
          qual_indices.append(i)
          threshold = thresholds[i]
          data_i = data[i,:]
          if sparse.issparse(data):
              data_i = data_i.toarray().flatten()
          assignments, means = poisson_cluster(data_i.reshape((1, cells)), 2)
          means = means.flatten()
          high_i = 1
          low_i = 0
          if means[0]>means[1]:
              high_i = 0
              low_i = 1
          # do a p-value test
          p_val = poisson_test(data_i[assignments==low_i], data_i[assignments==high_i], smoothing=smoothing)
          pvals[i] = p_val
          if p_val <= pval_threshold:
              genes_included.append(i)
          else:
              continue
          high_mean = np.median(data_i[assignments==high_i])
          low_mean = np.median(data_i[assignments==low_i]) + eps
          for k in range(clusters):
              if qualitative[i,k]>threshold:
                  output[i,k] = high_mean
              else:
                  output[i,k] = low_mean
      output = output[genes_included,:]
      pvals = pvals[genes_included]
      return output, pvals, genes_included