<commit_before>import os
import helpers
from helpers import dict_get, eval_surface_exposed_loop, \
                    chop_nterminal_peptide, log_stderr

def predict_surface_exposure_barrel(params, protein):
  # TODO: This is a placeholder for a function which will do something
  #       similar to predict_surface_exposure, but focussed on inferring 
  #       outer membrane beta barrel topology.
  #
  #       At the moment it's not implemented since I have not
  #       found an easily usable method which gives predictions
  #       accurate enough to be of any real use, particularly
  #       in terms of detecting barrel vs. non-barrel regions
  #       in multidomain BB OMPs (TMBETA-NET strand prediction is 
  #       implemented but it suffers from this problem). Probably 
  #       we will have to do domain detection independently, then 
  #       feed only the barrel domain to the strand predictor. 
  #       I've found PSI-PRED actually works reasonably well on OMP BBs - 
  #       this may be a simpler option over requiring a ProfTMB dependency 
  #       or interfacing with the transFold web service.
  #       
  #       Essentially, we should:
  #        * Move through the strand list in reverse.
  #        * Strand annotation alternates 'up' strand and 'down' strand
  #        * Loop annotation (starting with the C-terminal residue) alternates
  #          'inside' and 'outside'.
  #        * If everything is sane, we should finish on a down strand. If not,
  #          consider a rule to make an 'N-terminal up strand' become 
  #          an 'inside loop'
  #        * Sanity check on loop lengths ? 'Outside' loops should be on average
  #          longer than non-terminal 'inside' loops.
  #        * For alternative strand predictors (eg transFold, ProfTMB), which
  #          may specifically label inner and outer loops, we should obviously
  #          use those annotations directly.
  pass

def summary_table(params, proteins):
  """
  Returns a string representing a simple summary table of
  protein classifcations.
  """
  out = ""
  counts = {}
  for seqid in proteins:
    category = proteins[seqid]['category']
    
    if category not in counts:
      counts[category] = 0
    else:
      counts[category] += 1
      
  out += "\n\n# Number of proteins in each class:\n"
  for c in counts:
    out += "# %-15s %i\n" % (c, counts[c])
    
  return out
    

# TODO: Workflow:
#       These rules can be approached in two ways. Run every annotation plugin
#       on every sequence, then post-hoc apply these rules (potentially inefficient
#       and rude to web services). Alternatively the protocol can control the
#       annotation loop and only run certain analyses (eg BOMP) if a certain 
#       condition is met (eg is_signalp). 
#       * If is_lipop, 
#         check retention signal annotation (+2Asp) -> periplasmic inner or outer
#       * If is_signalp (or TAT) and no TMs and not lipop or barrel -> periplasmic
#       * If is_signalp, run OMPBB (BOMP) check. If OMPBB -> outer membrane bb
#       * If is_signalp, and has after the signal tranmembranes -> inner membrane
#       * If inner membrane, check for long cyto or peri loops: inner[+cyto][+peri]
#       * If not is_signalp and not is_lipop (and not TAT): -> cytoplasmic
def get_annotations(params):
  params['signalp4_organism'] = 'gram-'
  
  annotations = ['annotate_signalp4', 'annotate_lipop1', \
                 'annotate_tatfind_web']
  annotations += ['annotate_bomp_web']
  
  # TMBETA-NET knows to only run on predicted barrels
  # with the category 'OM(barrel)'
  if 'tmbeta' in dict_get(params, 'barrel_programs'):
    annotations.append('annotate_tmbeta_net_web')

  if dict_get(params, 'helix_programs'):
    if 'tmhmm' in params['helix_programs']:
      annotations.append('annotate_tmhmm')
    if 'memsat3' in params['helix_programs']:
      annotations.append('annotate_memsat3')

  # run some hmm profiles to detect features (eg Tat signal)
  annotations += ['annotate_hmmsearch3']
  params['hmm_profiles_dir'] = os.path.join(
      os.path.dirname(__file__), 'gram_neg_profiles')

  return annotations

def post_process_protein(params, protein, stringent=False):
    
  def has_tm_helix(protein):
    for program in params['helix_programs']:
      if dict_get(protein, '%s_helices' % program):
        return True
    return False

  # TODO: rather than this, lets classify each inner membrane protein as 
  # IM with optional large periplasmic or large cytoplasmic loops/domain
  # (eg possible classifications: IM+cyt, IM+peri, IM+cyt+peri )
  def has_periplasmic_loop(protein):
    for program in params['helix_programs']:
      if eval_surface_exposed_loop(
          protein['sequence_length'], 
          len(protein['%s_helices' % (program)]), 
          protein['%s_outer_loops' % (program)], 
          params['terminal_exposed_loop_min'], 
          params['internal_exposed_loop_min']):
        return True
    return False
  
  details = []
  category = "UNKNOWN"
  is_hmm_profile_match = dict_get(protein, 'hmmsearch')
  is_signalp = dict_get(protein, 'is_signalp')
  # TODO: in terms of logic, a Tat signal is essentially the same
  #       as a Sec (signalp) signal. Consider setting is_signal_pept = True
  #       if is_tatp = True, so all logic just looks at is_signal_pept
  is_tatfind = dict_get(protein, 'is_tatfind')
  is_lipop = dict_get(protein, 'is_lipop')
  
  is_signal_pept = False
  if is_signalp or is_tatfind or \
     ("Tat_PS51318" in dict_get(protein, 'hmmsearch')):
    is_signal_pept = True
  
  is_barrel = False
  if dict_get(protein, 'bomp'):
    is_barrel = True
    details += ['bomp']
  if dict_get(protein, 'tmbhunt'):
    is_barrel = True 
    details += ['tmbhunt']
  
  # if stringent, predicted OM barrels must also have a predicted
  # signal sequence 
  if stringent and is_signalp and is_barrel:
     category = 'OM(barrel)'
  elif is_barrel:
     category = 'OM(barrel)'
  protein['category'] = category

  # TODO: use thresholds for barrel annotation
  #if (dict_get(proteins[seqid], 'tmbhunt_prob') >= params['tmbhunt_cutoff']) or \
  #   (dict_get(proteins[seqid], 'bomp') >= params['bomp_cutoff']):
  #  protein['category'] = 'OM(barrel)'

  # set number of predicted OM barrel strands in details
  if (dict_get(protein, 'category') == 'OM(barrel)') and \<commit_msg>Added back in gram_neg barrel threshold rule.<commit_after>import os
import helpers
from helpers import dict_get, eval_surface_exposed_loop, \
                    chop_nterminal_peptide, log_stderr

def predict_surface_exposure_barrel(params, protein):
  # TODO: This is a placeholder for a function which will do something
  #       similar to predict_surface_exposure, but focussed on inferring 
  #       outer membrane beta barrel topology.
  #
  #       At the moment it's not implemented since I have not
  #       found an easily usable method which gives predictions
  #       accurate enough to be of any real use, particularly
  #       in terms of detecting barrel vs. non-barrel regions
  #       in multidomain BB OMPs (TMBETA-NET strand prediction is 
  #       implemented but it suffers from this problem). Probably 
  #       we will have to do domain detection independently, then 
  #       feed only the barrel domain to the strand predictor. 
  #       I've found PSI-PRED actually works reasonably well on OMP BBs - 
  #       this may be a simpler option over requiring a ProfTMB dependency 
  #       or interfacing with the transFold web service.
  #       
  #       Essentially, we should:
  #        * Move through the strand list in reverse.
  #        * Strand annotation alternates 'up' strand and 'down' strand
  #        * Loop annotation (starting with the C-terminal residue) alternates
  #          'inside' and 'outside'.
  #        * If everything is sane, we should finish on a down strand. If not,
  #          consider a rule to make an 'N-terminal up strand' become 
  #          an 'inside loop'
  #        * Sanity check on loop lengths ? 'Outside' loops should be on average
  #          longer than non-terminal 'inside' loops.
  #        * For alternative strand predictors (eg transFold, ProfTMB), which
  #          may specifically label inner and outer loops, we should obviously
  #          use those annotations directly.
  pass

def summary_table(params, proteins):
  """
  Returns a string representing a simple summary table of
  protein classifcations.
  """
  out = ""
  counts = {}
  for seqid in proteins:
    category = proteins[seqid]['category']
    
    if category not in counts:
      counts[category] = 0
    else:
      counts[category] += 1
      
  out += "\n\n# Number of proteins in each class:\n"
  for c in counts:
    out += "# %-15s %i\n" % (c, counts[c])
    
  return out
    

# TODO: Workflow:
#       These rules can be approached in two ways. Run every annotation plugin
#       on every sequence, then post-hoc apply these rules (potentially inefficient
#       and rude to web services). Alternatively the protocol can control the
#       annotation loop and only run certain analyses (eg BOMP) if a certain 
#       condition is met (eg is_signalp). 
#       * If is_lipop, 
#         check retention signal annotation (+2Asp) -> periplasmic inner or outer
#       * If is_signalp (or TAT) and no TMs and not lipop or barrel -> periplasmic
#       * If is_signalp, run OMPBB (BOMP) check. If OMPBB -> outer membrane bb
#       * If is_signalp, and has after the signal tranmembranes -> inner membrane
#       * If inner membrane, check for long cyto or peri loops: inner[+cyto][+peri]
#       * If not is_signalp and not is_lipop (and not TAT): -> cytoplasmic
def get_annotations(params):
  params['signalp4_organism'] = 'gram-'
  
  annotations = ['annotate_signalp4', 'annotate_lipop1', \
                 'annotate_tatfind_web']
  annotations += ['annotate_bomp_web']
  
  # TMBETA-NET knows to only run on predicted barrels
  # with the category 'OM(barrel)'
  if 'tmbeta' in dict_get(params, 'barrel_programs'):
    annotations.append('annotate_tmbeta_net_web')

  if dict_get(params, 'helix_programs'):
    if 'tmhmm' in params['helix_programs']:
      annotations.append('annotate_tmhmm')
    if 'memsat3' in params['helix_programs']:
      annotations.append('annotate_memsat3')

  # run some hmm profiles to detect features (eg Tat signal)
  annotations += ['annotate_hmmsearch3']
  params['hmm_profiles_dir'] = os.path.join(
      os.path.dirname(__file__), 'gram_neg_profiles')

  return annotations

def post_process_protein(params, protein, stringent=False):
    
  def has_tm_helix(protein):
    for program in params['helix_programs']:
      if dict_get(protein, '%s_helices' % program):
        return True
    return False

  # TODO: rather than this, lets classify each inner membrane protein as 
  # IM with optional large periplasmic or large cytoplasmic loops/domain
  # (eg possible classifications: IM+cyt, IM+peri, IM+cyt+peri )
  def has_periplasmic_loop(protein):
    for program in params['helix_programs']:
      if eval_surface_exposed_loop(
          protein['sequence_length'], 
          len(protein['%s_helices' % (program)]), 
          protein['%s_outer_loops' % (program)], 
          params['terminal_exposed_loop_min'], 
          params['internal_exposed_loop_min']):
        return True
    return False
  
  details = []
  category = "UNKNOWN"
  is_hmm_profile_match = dict_get(protein, 'hmmsearch')
  is_signalp = dict_get(protein, 'is_signalp')
  is_tatfind = dict_get(protein, 'is_tatfind')
  is_lipop = dict_get(protein, 'is_lipop')
  
  # in terms of most sublocalization logic, a Tat signal is essentially the 
  # same as a Sec (signalp) signal. We use is_signal_pept to denote that 
  # either is present.
  is_signal_pept = False
  if is_signalp or is_tatfind or \
     ("Tat_PS51318" in dict_get(protein, 'hmmsearch')):
    is_signal_pept = True
  
  is_barrel = False
  if dict_get(protein, 'bomp') >= params['bomp_cutoff']:
    is_barrel = True
    details += ['bomp']
  if dict_get(protein, 'tmbhunt_prob') >= params['tmbhunt_cutoff']:
    is_barrel = True 
    details += ['tmbhunt']
  
  # if stringent, predicted OM barrels must also have a predicted
  # signal sequence
  if stringent and is_signal_pept and is_barrel:
    category = 'OM(barrel)'
  elif is_barrel:
    category = 'OM(barrel)'
  protein['category'] = category

  # set number of predicted OM barrel strands in details
  if (dict_get(protein, 'category') == 'OM(barrel)') and \
{'fraction_non_alphanumeric': 0.08790508798471215, 'fraction_numerical': 0.0025479735647742654, 'mean_word_length': 3.0805717998700457, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4514853', 'n_tokens_mistral': 4031, 'n_tokens_neox': 3738, 'n_words': 1290}