docs: def get_orthology_matrix(self, outfile, sc, outdir=None,
                               pid_cutoff=None, bitscore_cutoff=None, evalue_cutoff=None,
                               force_rerun=False):
          """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:
              outfile (str): Filename with extension of the orthology matrix (ie. df_orthology.csv)
              outdir (str): Path to output of orthology matrix, default is ATLAS data_dir
              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.
  
          """
          # Output directory for orthology matrix
          if not outdir:
              outdir = self.data_dir
  
          ortho_matrix = op.join(outdir, outfile)
          if ssbio.utils.force_rerun(flag=force_rerun, outfile=ortho_matrix):
              if not sc:
                  raise ValueError('Please initialize SparkContext')
              ################################################################################################################
              # BIDIRECTIONAL BLAST
              def run_bidirectional_blast(strain_id, strain_genome_path,
                                          r_file=self.reference_gempro.genome_path,
                                          outdir2=self.sequences_by_organism_dir):
                  import ssbio.protein.sequence.utils.blast
  
                  # Run bidirectional BLAST
                  r_vs_g, g_vs_r = ssbio.protein.sequence.utils.blast.run_bidirectional_blast(reference=r_file,
                                                                                              other_genome=strain_genome_path,
                                                                                              dbtype='prot',
                                                                                              outdir=outdir2)
  
                  # Using the BLAST files, find the BBH
                  bbh = ssbio.protein.sequence.utils.blast.calculate_bbh(blast_results_1=r_vs_g, blast_results_2=g_vs_r,
                                                                         outdir=outdir2)
                  return strain_id, bbh
  
              log.info('Running bidirectional BLAST and finding best bidirectional hits (BBH)...')
              for_strains_rdd = [(k, v['genome_path']) for k,v in self.strain_infodict.items()]
              from random import shuffle
              shuffle(for_strains_rdd)
              strains_rdd = sc.parallelize(for_strains_rdd)
              result = strains_rdd.map(lambda x: run_bidirectional_blast(strain_id=x[0], strain_genome_path=x[1])).collect()
              bbh_files = dict(result)
              ################################################################################################################
  
              ################################################################################################################
              # ORTHOLOGY 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,
                                                                                        outname=outfile,
                                                                                        outdir=outdir,
                                                                                        force_rerun=force_rerun)
  
          log.info('Orthology matrix at {}. See the "df_orthology_matrix" attribute.'.format(ortho_matrix))
          self.df_orthology_matrix = pd.read_csv(ortho_matrix, index_col=0)
          self.df_orthology_matrix = self.df_orthology_matrix.rename_axis('gene').rename_axis("strain", axis="columns")