docs: def downsample(work_dir, sample_name, fastq_left_fpath, fastq_right_fpath, downsample_to, num_pairs=None):
      """ get N random headers from a fastq file without reading the
      whole thing into memory
      modified from: http://www.biostars.org/p/6544/
      """
      sample_name = sample_name or splitext(''.join(lc if lc == rc else '' for lc, rc in zip(fastq_left_fpath, fastq_right_fpath)))[0]
  
      l_out_fpath = make_downsampled_fpath(work_dir, fastq_left_fpath)
      r_out_fpath = make_downsampled_fpath(work_dir, fastq_right_fpath)
      if can_reuse(l_out_fpath, [fastq_left_fpath, fastq_right_fpath]):
          return l_out_fpath, r_out_fpath
  
      info('Processing ' + sample_name)
      if num_pairs is None:
          info(sample_name + ': counting number of reads in fastq...')
          num_pairs = _count_records_in_fastq(fastq_left_fpath)
      if num_pairs > LIMIT:
          info(sample_name + ' the number of reads is higher than ' + str(LIMIT) +
               ', sampling from only first ' + str(LIMIT))
          num_pairs = LIMIT
      info(sample_name + ': ' + str(num_pairs) + ' reads')
      num_downsample_pairs = int(downsample_to * num_pairs) if isinstance(downsample_to, float) else downsample_to
      if num_pairs <= num_downsample_pairs:
          info(sample_name + ': and it is less than ' + str(num_downsample_pairs) + ', so no downsampling.')
          return fastq_left_fpath, fastq_right_fpath
      else:
          info(sample_name + ': downsampling to ' + str(num_downsample_pairs))
          rand_records = sorted(random.sample(range(num_pairs), num_downsample_pairs))
  
      info('Opening ' + fastq_left_fpath)
      fh1 = open_gzipsafe(fastq_left_fpath)
      info('Opening ' + fastq_right_fpath)
      fh2 = open_gzipsafe(fastq_right_fpath) if fastq_right_fpath else None
  
      out_files = (l_out_fpath, r_out_fpath) if r_out_fpath else (l_out_fpath,)
  
      written_records = 0
      with file_transaction(work_dir, out_files) as tx_out_files:
          if isinstance(tx_out_files, six.string_types):
              tx_out_f1 = tx_out_files
          else:
              tx_out_f1, tx_out_f2 = tx_out_files
          info('Opening ' + str(tx_out_f1) + ' to write')
          sub1 = open_gzipsafe(tx_out_f1, "w")
          info('Opening ' + str(tx_out_f2) + ' to write')
          sub2 = open_gzipsafe(tx_out_f2, "w") if r_out_fpath else None
          rec_no = -1
          for rr in rand_records:
              while rec_no < rr:
                  rec_no += 1
                  for i in range(4): fh1.readline()
                  if fh2:
                      for i in range(4): fh2.readline()
              for i in range(4):
                  sub1.write(fh1.readline())
                  if sub2:
                      sub2.write(fh2.readline())
              written_records += 1
              if written_records % 10000 == 0:
                  info(sample_name + ': written ' + str(written_records) + ', rec_no ' + str(rec_no + 1))
              if rec_no > num_pairs:
                  info(sample_name + ' reached the limit of ' + str(num_pairs), ' read lines, stopping.')
                  break
          info(sample_name + ': done, written ' + str(written_records) + ', rec_no ' + str(rec_no))
          fh1.close()
          sub1.close()
          if fastq_right_fpath:
              fh2.close()
              sub2.close()
  
      info(sample_name + ': done downsampling, saved to ' + l_out_fpath + ' and ' + r_out_fpath + ', total ' + str(written_records) + ' paired reads written')
      return l_out_fpath, r_out_fpath