docs: def find_contamination(pair, output_folder, databases_folder, forward_id='_R1', threads=1, keep_files=False,
                         quality_cutoff=20, base_cutoff=2, base_fraction_cutoff=0.05, cgmlst_db=None, Xmx=None, tmpdir=None,
                         data_type='Illumina', use_rmlst=False):
      """
      This needs some documentation fairly badly, so here we go.
      :param pair: This has become a misnomer. If the input reads are actually paired, needs to be a list
      with the full filepath to forward reads at index 0 and full path to reverse reads at index 1.
      If reads are unpaired, should be a list of length 1 with the only entry being the full filepath to read set.
      :param output_folder: Folder where outputs (confindr log and report, and other stuff) will be stored.
      This will be created if it does not exist. (I think - should write a test that double checks this).
      :param databases_folder: Full path to folder where ConFindr's databases live. These files can be
      downloaded from figshare in .tar.gz format (https://ndownloader.figshare.com/files/11864267), and
      will be automatically downloaded if the script is run from the command line.
      :param forward_id: Identifier that marks reads as being in the forward direction for paired reads.
      Defaults to _R1
      :param threads: Number of threads to run analyses with. All parts of this pipeline scale pretty well,
      so more is better.
      :param keep_files: Boolean that says whether or not to keep temporary files.
      :param quality_cutoff: Integer of the phred score required to have a base count towards a multiallelic site.
      :param base_cutoff: Integer of number of bases needed to have a base be part of a multiallelic site.
      :param base_fraction_cutoff: Float of fraction of bases needed to have a base be part of a multiallelic site.
      If specified will be used in parallel with base_cutoff
      :param cgmlst_db: if None, we're using rMLST, if True, using some sort of custom cgMLST database. This requires some
      custom parameters.
      :param Xmx: if None, BBTools will use auto memory detection. If string, BBTools will use what's specified as their
      memory request.
      :param tmpdir: if None, any genus-specifc databases that need to be created will be written to ConFindr DB location.
      Otherwise, genus-specific databases will be written here.
      """
      if os.path.isfile(os.path.join(databases_folder, 'download_date.txt')):
          with open(os.path.join(databases_folder, 'download_date.txt')) as f:
              database_download_date = f.readline().rstrip()
      else:
          database_download_date = 'NA'
      log = os.path.join(output_folder, 'confindr_log.txt')
      if len(pair) == 2:
          sample_name = os.path.split(pair[0])[-1].split(forward_id)[0]
          paired = True
          logging.debug('Sample is paired. Sample name is {}'.format(sample_name))
      else:
          sample_name = os.path.split(pair[0])[-1].split('.')[0]
          paired = False
          logging.debug('Sample is unpaired. Sample name is {}'.format(sample_name))
      sample_tmp_dir = os.path.join(output_folder, sample_name)
      if not os.path.isdir(sample_tmp_dir):
          os.makedirs(sample_tmp_dir)
  
      logging.info('Checking for cross-species contamination...')
      if paired:
          genus = find_cross_contamination(databases_folder, pair, tmpdir=sample_tmp_dir, log=log, threads=threads)
      else:
          genus = find_cross_contamination_unpaired(databases_folder, reads=pair[0], tmpdir=sample_tmp_dir, log=log, threads=threads)
      if len(genus.split(':')) > 1:
          write_output(output_report=os.path.join(output_folder, 'confindr_report.csv'),
                       sample_name=sample_name,
                       multi_positions=0,
                       genus=genus,
                       percent_contam='NA',
                       contam_stddev='NA',
                       total_gene_length=0,
                       database_download_date=database_download_date)
          logging.info('Found cross-contamination! Skipping rest of analysis...\n')
          if keep_files is False:
              shutil.rmtree(sample_tmp_dir)
          return
      # Setup genus-specific databases, if necessary.
      if cgmlst_db is not None:
          # Sanity check that the DB specified is actually a file, otherwise, quit with appropriate error message.
          if not os.path.isfile(cgmlst_db):
              logging.error('ERROR: Specified cgMLST file ({}) does not exist. Please check the path and try again.'.format(cgmlst_db))
              quit(code=1)
          sample_database = cgmlst_db
      else:
          db_folder = databases_folder if tmpdir is None else tmpdir
          if not os.path.isdir(db_folder):
              os.makedirs(db_folder)
          if genus != 'NA':
              # Logic here is as follows: users can either have both rMLST databases, which cover all of bacteria,
              # cgmlst-derived databases, which cover only Escherichia, Salmonella, and Listeria (may add more at some
              # point), or they can have both. They can also set priority to either always use rMLST, or to use my
              # core-genome derived stuff and fall back on rMLST if they're trying to look at a genus I haven't created
              # a scheme for.
  
              # In the event rmlst databases have priority, always use them.
              if use_rmlst is True:
                  sample_database = os.path.join(db_folder, '{}_db.fasta'.format(genus))
                  if not os.path.isfile(sample_database):
                      if os.path.isfile(os.path.join(db_folder, 'gene_allele.txt')) and os.path.isfile(os.path.join(db_folder, 'rMLST_combined.fasta')):
                          logging.info('Setting up genus-specific database for genus {}...'.format(genus))
                          allele_list = find_genusspecific_allele_list(os.path.join(db_folder, 'gene_allele.txt'), genus)
                          setup_allelespecific_database(fasta_file=sample_database,
                                                        database_folder=db_folder,
                                                        allele_list=allele_list)
              else:
                  # Check if a cgderived database is available. If not, try to use rMLST database.
                  sample_database = os.path.join(db_folder, '{}_db_cgderived.fasta'.format(genus))
                  if not os.path.isfile(sample_database):
                      sample_database = os.path.join(db_folder, '{}_db.fasta'.format(genus))
                      # Create genus specific database if it doesn't already exist and we have the necessary rMLST files.
                      if os.path.isfile(os.path.join(db_folder, 'rMLST_combined.fasta')) and os.path.isfile(os.path.join(db_folder, 'gene_allele.txt')) and not os.path.isfile(sample_database):
                          logging.info('Setting up genus-specific database for genus {}...'.format(genus))
                          allele_list = find_genusspecific_allele_list(os.path.join(db_folder, 'gene_allele.txt'), genus)
                          setup_allelespecific_database(fasta_file=sample_database,
                                                        database_folder=db_folder,
                                                        allele_list=allele_list)
  
          else:
              sample_database = os.path.join(db_folder, 'rMLST_combined.fasta')
  
      # If a user has gotten to this point and they don't have any database available to do analysis because
      # they don't have rMLST downloaded and we don't have a cg-derived database available, boot them with a helpful
      # message.
      if not os.path.isfile(sample_database):
          write_output(output_report=os.path.join(output_folder, 'confindr_report.csv'),
                       sample_name=sample_name,
                       multi_positions=0,
                       genus=genus,
                       percent_contam='NA',
                       contam_stddev='NA',
                       total_gene_length=0,
                       database_download_date=database_download_date)
          logging.info('Did not find databases for genus {genus}. You can download the rMLST database to get access to all '
                       'genera (see https://olc-bioinformatics.github.io/ConFindr/install/). Alternatively, if you have a '
                       'high-quality core-genome derived database for your genome of interest, we would be happy to '
                       'add it - open an issue at https://github.com/OLC-Bioinformatics/ConFindr/issues with the '
                       'title "Add genus-specific database: {genus}"\n'.format(genus=genus))
          if keep_files is False:
              shutil.rmtree(sample_tmp_dir)
          return
  
      # Extract rMLST reads and quality trim.
      logging.info('Extracting conserved core genes...')
      if paired:
          if Xmx is None:
              out, err, cmd = bbtools.bbduk_bait(reference=sample_database,
                                                 forward_in=pair[0],
                                                 reverse_in=pair[1],
                                                 forward_out=os.path.join(sample_tmp_dir, 'rmlst_R1.fastq.gz'),
                                                 reverse_out=os.path.join(sample_tmp_dir, 'rmlst_R2.fastq.gz'),
                                                 threads=threads,
                                                 returncmd=True)
          else:
              out, err, cmd = bbtools.bbduk_bait(reference=sample_database,
                                                 forward_in=pair[0],
                                                 reverse_in=pair[1],
                                                 forward_out=os.path.join(sample_tmp_dir, 'rmlst_R1.fastq.gz'),
                                                 reverse_out=os.path.join(sample_tmp_dir, 'rmlst_R2.fastq.gz'),
                                                 threads=threads,
                                                 Xmx=Xmx,
                                                 returncmd=True)
      else:
          if data_type == 'Nanopore':
              forward_out = os.path.join(sample_tmp_dir, 'trimmed.fastq.gz')
          else:
              forward_out = os.path.join(sample_tmp_dir, 'rmlst.fastq.gz')
          if Xmx is None:
              out, err, cmd = bbtools.bbduk_bait(reference=sample_database, forward_in=pair[0],
                                                 forward_out=forward_out,
                                                 returncmd=True, threads=threads)
          else:
              out, err, cmd = bbtools.bbduk_bait(reference=sample_database, forward_in=pair[0],
                                                 forward_out=forward_out, Xmx=Xmx,
                                                 returncmd=True, threads=threads)
  
      write_to_logfile(log, out, err, cmd)
      logging.info('Quality trimming...')
      if data_type == 'Illumina':
          if paired:
              if Xmx is None:
                  out, err, cmd = bbtools.bbduk_trim(forward_in=os.path.join(sample_tmp_dir, 'rmlst_R1.fastq.gz'),
                                                     reverse_in=os.path.join(sample_tmp_dir, 'rmlst_R2.fastq.gz'),
                                                     forward_out=os.path.join(sample_tmp_dir, 'trimmed_R1.fastq.gz'),
                                                     reverse_out=os.path.join(sample_tmp_dir, 'trimmed_R2.fastq.gz'),
                                                     threads=str(threads), returncmd=True)
              else:
                  out, err, cmd = bbtools.bbduk_trim(forward_in=os.path.join(sample_tmp_dir, 'rmlst_R1.fastq.gz'),
                                                     reverse_in=os.path.join(sample_tmp_dir, 'rmlst_R2.fastq.gz'),
                                                     forward_out=os.path.join(sample_tmp_dir, 'trimmed_R1.fastq.gz'),
                                                     reverse_out=os.path.join(sample_tmp_dir, 'trimmed_R2.fastq.gz'), Xmx=Xmx,
                                                     threads=str(threads), returncmd=True)
  
          else:
              if Xmx is None:
                  out, err, cmd = bbtools.bbduk_trim(forward_in=os.path.join(sample_tmp_dir, 'rmlst.fastq.gz'),
                                                     forward_out=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                                     returncmd=True, threads=threads)
              else:
                  out, err, cmd = bbtools.bbduk_trim(forward_in=os.path.join(sample_tmp_dir, 'rmlst.fastq.gz'),
                                                     forward_out=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                                     returncmd=True, threads=threads, Xmx=Xmx)
          write_to_logfile(log, out, err, cmd)
  
      logging.info('Detecting contamination...')
      # Now do mapping in two steps - first, map reads back to database with ambiguous reads matching all - this
      # will be used to get a count of number of reads aligned to each gene/allele so we can create a custom rmlst file
      # with only the most likely allele for each gene.
      if not os.path.isfile(sample_database + '.fai'):  # Don't bother re-indexing, this only needs to happen once.
          pysam.faidx(sample_database)
      kma_database = sample_database.replace('.fasta', '') + '_kma'
      kma_report = os.path.join(sample_tmp_dir, 'kma_rmlst')
      if not os.path.isfile(kma_database + '.name'):  # The .name is one of the files KMA creates when making a database.
          cmd = 'kma index -i {} -o {}'.format(sample_database, kma_database)  # NOTE: Need KMA >=1.2.0 for this to work.
          out, err = run_cmd(cmd)
          write_to_logfile(log, out, err, cmd)
  
      # Run KMA.
      if paired:
          cmd = 'kma -ipe {forward_in} {reverse_in} -t_db {kma_database} -o {kma_report} ' \
                '-t {threads}'.format(forward_in=os.path.join(sample_tmp_dir, 'trimmed_R1.fastq.gz'),
                                      reverse_in=os.path.join(sample_tmp_dir, 'trimmed_R2.fastq.gz'),
                                      kma_database=kma_database,
                                      kma_report=kma_report,
                                      threads=threads)
          out, err = run_cmd(cmd)
          write_to_logfile(log, out, err, cmd)
      else:
          if data_type == 'Illumina':
              cmd = 'kma -i {input_reads} -t_db {kma_database} -o {kma_report} ' \
                    '-t {threads}'.format(input_reads=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                          kma_database=kma_database,
                                          kma_report=kma_report,
                                          threads=threads)
          else:
              # Recommended Nanopore settings from KMA repo: https://bitbucket.org/genomicepidemiology/kma
              cmd = 'kma -i {input_reads} -t_db {kma_database} -o {kma_report} -mem_mode -mp 20 -mrs 0.0 -bcNano ' \
                    '-t {threads}'.format(input_reads=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                          kma_database=kma_database,
                                          kma_report=kma_report,
                                          threads=threads)
          out, err = run_cmd(cmd)
          write_to_logfile(log, out, err, cmd)
  
      rmlst_report = os.path.join(output_folder, sample_name + '_rmlst.csv')
      gene_alleles = find_rmlst_type(kma_report=kma_report + '.res',
                                     rmlst_report=rmlst_report)
  
      with open(os.path.join(sample_tmp_dir, 'rmlst.fasta'), 'w') as f:
          for contig in SeqIO.parse(sample_database, 'fasta'):
              if contig.id in gene_alleles:
                  f.write('>{}\n'.format(contig.id))
                  f.write(str(contig.seq) + '\n')
  
      rmlst_gene_length = find_total_sequence_length(os.path.join(sample_tmp_dir, 'rmlst.fasta'))
      logging.debug('Total gene length is {}'.format(rmlst_gene_length))
      # Second step of mapping - Do a mapping of our baited reads against a fasta file that has only one allele per
      # rMLST gene.
      pysam.faidx(os.path.join(sample_tmp_dir, 'rmlst.fasta'))
      if paired:
          cmd = 'bbmap.sh ref={ref} in={forward_in} in2={reverse_in} out={outbam} threads={threads} mdtag ' \
                'nodisk'.format(ref=os.path.join(sample_tmp_dir, 'rmlst.fasta'),
                                forward_in=os.path.join(sample_tmp_dir, 'trimmed_R1.fastq.gz'),
                                reverse_in=os.path.join(sample_tmp_dir, 'trimmed_R2.fastq.gz'),
                                outbam=os.path.join(sample_tmp_dir, 'out_2.bam'),
                                threads=threads)
          if cgmlst_db is not None:
              # Lots of core genes seem to have relatives within a genome that are at ~70 percent identity - this means
              # that reads that shouldn't really map do, and cause false positives. Adding in this subfilter means that
              # reads can only have one mismatch, so they actually have to be from the right gene for this to work.
              cmd += ' subfilter=1'
          if Xmx:
              cmd += ' -Xmx{}'.format(Xmx)
          out, err = run_cmd(cmd)
          write_to_logfile(log, out, err, cmd)
      else:
          if data_type == 'Illumina':
              cmd = 'bbmap.sh ref={ref} in={forward_in} out={outbam} threads={threads} mdtag ' \
                    'nodisk'.format(ref=os.path.join(sample_tmp_dir, 'rmlst.fasta'),
                                    forward_in=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                    outbam=os.path.join(sample_tmp_dir, 'out_2.bam'),
                                    threads=threads)
              if cgmlst_db is not None:
                  # Lots of core genes seem to have relatives within a genome that are at ~70 percent identity - this means
                  # that reads that shouldn't really map do, and cause false positives. Adding in this subfilter means that
                  # reads can only have one mismatch, so they actually have to be from the right gene for this to work.
                  cmd += ' subfilter=1'
              if Xmx:
                  cmd += ' -Xmx{}'.format(Xmx)
              out, err = run_cmd(cmd)
              write_to_logfile(log, out, err, cmd)
          else:
              cmd = 'minimap2 --MD -t {threads} -ax map-ont {ref} {reads} ' \
                    '> {outsam}'.format(ref=os.path.join(sample_tmp_dir, 'rmlst.fasta'),
                                        reads=os.path.join(sample_tmp_dir, 'trimmed.fastq.gz'),
                                        outsam=os.path.join(sample_tmp_dir, 'out_2.sam'),
                                        threads=threads)
              out, err = run_cmd(cmd)
              write_to_logfile(log, out, err, cmd)
              outbam = os.path.join(sample_tmp_dir, 'out_2.bam')
              # Apparently have to perform equivalent of a touch on this file for this to work.
              fh = open(outbam, 'w')
              fh.close()
              pysam.view('-b', '-o', outbam, os.path.join(sample_tmp_dir, 'out_2.sam'), save_stdout=outbam)
      pysam.sort('-o', os.path.join(sample_tmp_dir, 'contamination.bam'), os.path.join(sample_tmp_dir, 'out_2.bam'))
      pysam.index(os.path.join(sample_tmp_dir, 'contamination.bam'))
      # Now find number of multi-positions for each rMLST gene/allele combination
      multi_positions = 0
  
      # Run the BAM parsing in parallel! Some refactoring of the code would likely be a good idea so this
      # isn't quite so ugly, but it works.
      p = multiprocessing.Pool(processes=threads)
      bamfile_list = [os.path.join(sample_tmp_dir, 'contamination.bam')] * len(gene_alleles)
      # bamfile_list = [os.path.join(sample_tmp_dir, 'rmlst.bam')] * len(gene_alleles)
      reference_fasta_list = [os.path.join(sample_tmp_dir, 'rmlst.fasta')] * len(gene_alleles)
      quality_cutoff_list = [quality_cutoff] * len(gene_alleles)
      base_cutoff_list = [base_cutoff] * len(gene_alleles)
      base_fraction_list = [base_fraction_cutoff] * len(gene_alleles)
      multibase_dict_list = list()
      report_write_list = list()
      for multibase_dict, report_write in p.starmap(read_contig, zip(gene_alleles, bamfile_list, reference_fasta_list, quality_cutoff_list, base_cutoff_list, base_fraction_list), chunksize=1):
          multibase_dict_list.append(multibase_dict)
          report_write_list.append(report_write)
      p.close()
      p.join()
  
      # Write out report info.
      report_file = os.path.join(output_folder, sample_name + '_contamination.csv')
      with open(report_file, 'w') as r:
          r.write('{reference},{position},{bases},{coverage}\n'.format(reference='Gene',
                                                                       position='Position',
                                                                       bases='Bases',
                                                                       coverage='Coverage'))
          for item in report_write_list:
              for contamination_info in item:
                  r.write(contamination_info)
  
      # Total up the number of multibase positions.
      for multibase_position_dict in multibase_dict_list:
          multi_positions += sum([len(snp_positions) for snp_positions in multibase_position_dict.values()])
      if cgmlst_db is None:
          snp_cutoff = int(rmlst_gene_length/10000) + 1
      else:
          snp_cutoff = 10
      if multi_positions >= snp_cutoff:
          percent_contam, contam_stddev = estimate_percent_contamination(contamination_report_file=report_file)
      else:
          percent_contam = 0
          contam_stddev = 0
      logging.info('Done! Number of contaminating SNVs found: {}\n'.format(multi_positions))
      write_output(output_report=os.path.join(output_folder, 'confindr_report.csv'),
                   sample_name=sample_name,
                   multi_positions=multi_positions,
                   genus=genus,
                   percent_contam=percent_contam,
                   contam_stddev=contam_stddev,
                   total_gene_length=rmlst_gene_length,
                   snp_cutoff=snp_cutoff,
                   cgmlst=cgmlst_db,
                   database_download_date=database_download_date)
      if keep_files is False:
          shutil.rmtree(sample_tmp_dir)