docs: def get_orthology_matrix(self, pid_cutoff=None, bitscore_cutoff=None, evalue_cutoff=None, filter_condition='OR',
                               remove_strains_with_no_orthology=True,
                               remove_strains_with_no_differences=False,
                               remove_genes_not_in_base_model=True):
          """Create the orthology matrix by finding best bidirectional BLAST hits. Genes = rows, strains = columns
  
          Runs run_makeblastdb, run_bidirectional_blast, and calculate_bbh for protein sequences.
  
          Args:
              pid_cutoff (float): Minimum percent identity between BLAST hits to filter for in the range [0, 100]
              bitscore_cutoff (float): Minimum bitscore allowed between BLAST hits
              evalue_cutoff (float): Maximum E-value allowed between BLAST hits
              filter_condition (str): 'OR' or 'AND', how to combine cutoff filters. 'OR' gives more results since it
              is less stringent, as you will be filtering for hits with (>80% PID or >30 bitscore or <0.0001 evalue).
              remove_strains_with_no_orthology (bool): Remove strains which have no orthologous genes found
              remove_strains_with_no_differences (bool): Remove strains which have all the same genes as the base model.
                  Default is False because since orthology is found using a PID cutoff, all genes may be present but
                  differences may be on the sequence level.
              remove_genes_not_in_base_model (bool): Remove genes from the orthology matrix which are not present in our
                  base model. This happens if we use a genome file for our model that has other genes in it.
  
          Returns:
              DataFrame: Orthology matrix calculated from best bidirectional BLAST hits.
  
          """
          # TODO: document and test other cutoffs
  
          # Get the path to the reference genome
          r_file = self.reference_gempro.genome_path
  
          bbh_files = {}
  
          log.info('Running bidirectional BLAST and finding best bidirectional hits (BBH)...')
          for strain_gempro in tqdm(self.strains):
              g_file = strain_gempro.genome_path
  
              # Run bidirectional BLAST
              log.debug('{} vs {}: Running bidirectional BLAST'.format(self.reference_gempro.id, strain_gempro.id))
              r_vs_g, g_vs_r = ssbio.protein.sequence.utils.blast.run_bidirectional_blast(reference=r_file,
                                                                                          other_genome=g_file,
                                                                                          dbtype='prot',
                                                                                          outdir=self.sequences_by_organism_dir)
  
              # Using the BLAST files, find the BBH
              log.debug('{} vs {}: Finding BBHs'.format(self.reference_gempro.id, strain_gempro.id))
              bbh = ssbio.protein.sequence.utils.blast.calculate_bbh(blast_results_1=r_vs_g, blast_results_2=g_vs_r,
                                                                     outdir=self.sequences_by_organism_dir)
              bbh_files[strain_gempro.id] = bbh
  
          # Make the orthologous genes matrix
          log.info('Creating orthology matrix from BBHs...')
          ortho_matrix = ssbio.protein.sequence.utils.blast.create_orthology_matrix(r_name=self.reference_gempro.id,
                                                                                    genome_to_bbh_files=bbh_files,
                                                                                    pid_cutoff=pid_cutoff,
                                                                                    bitscore_cutoff=bitscore_cutoff,
                                                                                    evalue_cutoff=evalue_cutoff,
                                                                                    filter_condition=filter_condition,
                                                                                    outname='{}_{}_orthology.csv'.format(self.reference_gempro.id, 'prot'),
                                                                                    outdir=self.data_dir)
  
          log.info('Saved orthology matrix at {}. See the "df_orthology_matrix" attribute.'.format(ortho_matrix))
          self.df_orthology_matrix = pd.read_csv(ortho_matrix, index_col=0)
  
          # Filter the matrix to genes only in our analysis, and also check for strains with no differences or no orthologous genes
          self._filter_orthology_matrix(remove_strains_with_no_orthology=remove_strains_with_no_orthology,
                                             remove_strains_with_no_differences=remove_strains_with_no_differences,
                                             remove_genes_not_in_base_model=remove_genes_not_in_base_model)