docs: def filter_genes_and_strains(self, remove_genes_not_in_reference_model=True,
                                   remove_strains_with_no_orthology=True, remove_strains_with_no_differences=False,
                                   custom_keep_strains=None, custom_keep_genes=None):
          """Filters the analysis by keeping a subset of strains or genes based on certain criteria.
  
          Args:
              remove_genes_not_in_reference_model (bool): Remove genes from reference model not in orthology matrix
              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.
              custom_keep_genes (list): List of gene IDs to keep in analysis
              custom_keep_strains (list): List of strain IDs to keep in analysis
  
          """
  
          if len(self.df_orthology_matrix) == 0:
              raise RuntimeError('Empty orthology matrix, please calculate first!')
  
          reference_strain_gene_ids = [x.id for x in self.reference_gempro.genes]
          initial_num_genes = len(reference_strain_gene_ids)
  
          initial_num_strains = len(self.strain_ids)
  
          # Gene filtering
          to_remove_genes = []
          if custom_keep_genes:
              to_remove_genes.extend([x for x in reference_strain_gene_ids if x not in custom_keep_genes])
          if remove_genes_not_in_reference_model:
              to_remove_genes.extend([x for x in reference_strain_gene_ids if x not in self.df_orthology_matrix.index.tolist()])
  
          to_remove_genes = list(set(to_remove_genes))
          if self.reference_gempro.model:
              cobra.manipulation.delete_model_genes(self.reference_gempro.model, to_remove_genes)
          else:
              for g_id in to_remove_genes:
                  self.reference_gempro.genes.get_by_id(g_id).functional = False
  
          # Create new orthology matrix with only our genes of interest
          new_gene_subset = [x.id for x in self.reference_gempro.functional_genes]
          tmp_new_orthology_matrix = self.df_orthology_matrix[self.df_orthology_matrix.index.isin(new_gene_subset)]
  
          # Strain filtering
          if custom_keep_strains or remove_strains_with_no_orthology or remove_strains_with_no_differences:
              for strain_id in self.strain_ids:
                  if custom_keep_strains:
                      if strain_id not in custom_keep_strains:
                          self.strain_ids.remove(strain_id)
                          continue
  
                  if remove_strains_with_no_orthology:
                      if strain_id not in tmp_new_orthology_matrix.columns:
                          self.strain_ids.remove(strain_id)
                          log.info('{}: no orthologous genes found for this strain, removed from analysis.'.format(strain_id))
                          continue
                      elif tmp_new_orthology_matrix[strain_id].isnull().all():
                          self.strain_ids.remove(strain_id)
                          log.info('{}: no orthologous genes found for this strain, removed from analysis.'.format(strain_id))
                          continue
  
                  if remove_strains_with_no_differences:
                      not_in_strain = tmp_new_orthology_matrix[pd.isnull(tmp_new_orthology_matrix[strain_id])][strain_id].index.tolist()
                      if len(not_in_strain) == 0:
                          self.strain_ids.remove(strain_id)
                          log.info('{}: strain has no differences from the base, removed from analysis.')
                          continue
  
          log.info('{} genes to be analyzed, originally {}'.format(len(self.reference_gempro.functional_genes), initial_num_genes))
          log.info('{} strains to be analyzed, originally {}'.format(len(self.strain_ids), initial_num_strains))