code: def analyze_dists(dists, cutoff=1000, alpha=.1):
      
      peak0 = [d for d in dists if d < cutoff]
      peak1 = [d for d in dists if d >= cutoff]
      c0, c1 = len(peak0), len(peak1)
      logging.debug("Component counts: {0} {1}".format(c0, c1))
      if c0 == 0 or c1 == 0 or float(c1) / len(dists) < alpha:
          logging.debug("Single peak identified ({0} / {1} < {2})".\
                          format(c1, len(dists), alpha))
          return np.median(dists)
  
      peak0_median = np.median(peak0)
      peak1_median = np.median(peak1)
      logging.debug("Dual peaks identified: {0}bp ({1}), {2}bp ({3}) (selected)".\
                  format(int(peak0_median), c0, int(peak1_median), c1))
  
      return peak1_median
docs: The dists can show bimodal distribution if they come from a mate-pair
      library. Assume bimodal distribution and then separate the two peaks. Based
      on the percentage in each peak, we can decide if it is indeed one peak or
      two peaks, and report the median respectively.