docs: def pairwise_sequence_alignment(a_seq, b_seq, engine, a_seq_id=None, b_seq_id=None,
                                  gapopen=10, gapextend=0.5,
                                  outfile=None, outdir=None, force_rerun=False):
      """Run a global pairwise sequence alignment between two sequence strings.
  
      Args:
          a_seq (str, Seq, SeqRecord, SeqProp): Reference sequence
          b_seq (str, Seq, SeqRecord, SeqProp): Sequence to be aligned to reference
          engine (str): `biopython` or `needle` - which pairwise alignment program to use
          a_seq_id (str): Reference sequence ID. If not set, is "a_seq"
          b_seq_id (str): Sequence to be aligned ID. If not set, is "b_seq"
          gapopen (int): Only for `needle` - Gap open penalty is the score taken away when a gap is created
          gapextend (float): Only for `needle` - Gap extension penalty is added to the standard gap penalty for each 
              base or residue in the gap
          outfile (str): Only for `needle` - name of output file. If not set, is {id_a}_{id_b}_align.txt
          outdir (str): Only for `needle` - Path to output directory. Default is the current directory.
          force_rerun (bool): Only for `needle` - Default False, set to True if you want to rerun the alignment 
              if outfile exists.
  
      Returns:
          MultipleSeqAlignment: Biopython object to represent an alignment
  
      """
      engine = engine.lower()
  
      if engine not in ['biopython', 'needle']:
          raise ValueError('{}: invalid engine'.format(engine))
  
      if not a_seq_id:
          a_seq_id = 'a_seq'
      if not b_seq_id:
          b_seq_id = 'b_seq'
  
      a_seq = ssbio.protein.sequence.utils.cast_to_str(a_seq)
      b_seq = ssbio.protein.sequence.utils.cast_to_str(b_seq)
  
      if engine == 'biopython':
          # TODO: allow different matrices? needle uses blosum62 by default, how to change that?
          # TODO: how to define gap open/extend when using matrix in biopython global alignment?
          log.warning('Gap penalties not implemented in Biopython yet')
  
          blosum62 = matlist.blosum62
          alignments = pairwise2.align.globaldx(a_seq, b_seq, blosum62)  # TODO: add gap penalties
          best_alignment = alignments[0]
  
          a = ssbio.protein.sequence.utils.cast_to_seq_record(best_alignment[0], id=a_seq_id)
          b = ssbio.protein.sequence.utils.cast_to_seq_record(best_alignment[1], id=b_seq_id)
          alignment = MultipleSeqAlignment([a, b], annotations={"score": best_alignment[2],
                                                                "start": best_alignment[3],
                                                                "end"  : best_alignment[4]})
          alignment.annotations['percent_identity'] = get_percent_identity(best_alignment[0], best_alignment[1]) * 100
  
          return alignment
  
      if engine == 'needle':
          alignment_file = run_needle_alignment(seq_a=a_seq, seq_b=b_seq, gapopen=gapopen, gapextend=gapextend,
                                                write_outfile=True,  # Has to be true, AlignIO parses files on disk
                                                outdir=outdir, outfile=outfile, force_rerun=force_rerun)
          log.debug('Needle alignment at {}'.format(alignment_file))
  
          if not op.exists(alignment_file):
              raise ValueError('{}: needle alignment file does not exist'.format(alignment_file))
  
          # Use AlignIO to parse the needle alignment, alignments[0] is the first alignment (the only one in pairwise)
          # alignments = list(AlignIO.parse(alignment_file, "emboss"))
          # alignment = alignments[0]
          alignment = needle_statistics_alignio(alignment_file)
  
          # Rename the sequence IDs
          alignment[0].id = a_seq_id
          alignment[1].id = b_seq_id
  
          # # Add needle statistics as annotations in the alignment object
          # stats = needle_statistics(alignment_file)
          # alignment_ids = list(stats.keys())
          # if len(alignment_ids) > 1:
          #     raise ValueError('Needle alignment file contains more than one pairwise alignment')
          # needle_id = alignment_ids[0]
          # alignment.annotations['percent_identity'] = stats[needle_id]['percent_identity']
          # alignment.annotations['percent_similarity'] = stats[needle_id]['percent_similarity']
          # alignment.annotations['percent_gaps'] = stats[needle_id]['percent_gaps']
          # alignment.annotations['score'] = stats[needle_id]['score']
  
          return alignment