Package csb :: Package bio :: Package hmm :: Module pseudocounts

[frames] | no frames]

# Source Code for Module csb.bio.hmm.pseudocounts

      1  import sys
      2  import csb.core
      3  import csb.bio.sequence as sequence
      4
      5  from csb.bio.hmm import States, ScoreUnits, Transition, State
      6
      7
      8  GONNET = [10227, 3430, 2875, 3869, 1625, 2393, 4590, 6500, 2352, 3225, 5819, 4172, 1435,
      9            1579, 3728, 4610, 6264, 418, 1824, 5709, 3430, 7780, 2209, 2589, 584, 2369,
     10            3368, 3080, 2173, 1493, 3093, 5701, 763, 859, 1893, 2287, 3487, 444, 1338, 2356,
     11            2875, 2209, 3868, 3601, 501, 1541, 2956, 3325, 1951, 1065, 2012, 2879, 532, 688,
     12            1480, 2304, 3204, 219, 1148, 1759, 3869, 2589, 3601, 8618, 488, 2172, 6021, 4176,
     13            2184, 1139, 2151, 3616, 595, 670, 2086, 2828, 3843, 204, 1119, 2015, 1625, 584,
     14            501, 488, 5034, 355, 566, 900, 516, 741, 1336, 591, 337, 549, 419, 901,
     15            1197, 187, 664, 1373, 2393, 2369, 1541, 2172, 355, 1987, 2891, 1959, 1587, 1066,
     16            2260, 2751, 570, 628, 1415, 1595, 2323, 219, 871, 1682, 4590, 3368, 2956, 6021,
     17            566, 2891, 8201, 3758, 2418, 1624, 3140, 4704, 830, 852, 2418, 2923, 4159, 278, 1268, 2809,
     18            6500, 3080, 3325, 4176, 900, 1959, 3758, 26066, 2016, 1354, 2741, 3496, 741, 797, 2369,
     19            3863, 4169, 375, 1186, 2569, 2352, 2173, 1951, 2184, 516, 1587, 2418, 2016, 5409,
     20            1123, 2380, 2524, 600, 1259, 1298, 1642, 2446, 383, 876, 1691, 3225, 1493, 1065,
     21            1139, 741, 1066, 1624, 1354, 1123, 6417, 9630, 1858, 1975, 2225, 1260, 1558, 3131,
     22            417, 1697, 7504, 5819, 3093, 2012, 2151, 1336, 2260, 3140, 2741, 2380, 9630, 25113,
     23            3677, 4187, 5540, 2670, 2876, 5272, 1063, 3945, 11005, 4172, 5701, 2879, 3616, 591,
     24            2751, 4704, 3496, 2524, 1858, 3677, 7430, 949, 975, 2355, 2847, 4340, 333, 1451, 2932,
     25            1435, 763, 532, 595, 337, 570, 830, 741, 600, 1975, 4187, 949, 1300, 1111, 573,
     26            743, 1361, 218, 828, 2310, 1579, 859, 688, 670, 549, 628, 852, 797, 1259, 2225,
     27            5540, 975, 1111, 6126, 661, 856, 1498, 1000, 4464, 2602, 3728, 1893, 1480, 2086, 419,
     28            1415, 2418, 2369, 1298, 1260, 2670, 2355, 573, 661, 11834, 2320, 3300, 179, 876, 2179,
     29            4610, 2287, 2304, 2828, 901, 1595, 2923, 3863, 1642, 1558, 2876, 2847, 743, 856, 2320,
     30            3611, 4686, 272, 1188, 2695, 6264, 3487, 3204, 3843, 1197, 2323, 4159, 4169, 2446, 3131,
     31            5272, 4340, 1361, 1498, 3300, 4686, 8995, 397, 1812, 5172, 418, 444, 219, 204, 187,
     32            219, 278, 375, 383, 417, 1063, 333, 218, 1000, 179, 272, 397, 4101, 1266, 499,
     33            1824, 1338, 1148, 1119, 664, 871, 1268, 1186, 876, 1697, 3945, 1451, 828, 4464, 876,
     34            1188, 1812, 1266, 9380, 2227, 5709, 2356, 1759, 2015, 1373, 1682, 2809, 2569, 1691, 7504,
     35            11005, 2932, 2310, 2602, 2179, 2695, 5172, 499, 2227.0, 11569.0]
     36  """
     37  Gonnet matrix frequencies taken from HHpred
     38  """
     39
     40
     41 -class PseudocountBuilder(object):
     42      """
     43      Constructs profile HMMs with pseudocounts.
     44      """
     45
     46 -    def __init__(self, hmm):
     47          self._hmm = hmm
     48
     49      @property
     50 -    def hmm(self):
     51          return self._hmm
     52
     53 -    def add_emission_pseudocounts(self, tau=0.1, pca=2.5, pcb=0.5, pcc=1.0):
     54          """
     55          Port from HHpred, it uses the conditional background probabilities,
     56          inferred from the Gonnet matrix.
     57
     58          @param tau: admission weight, i.e how much of the final score is
     59                      determined by the background probabilities.
     60                      0.0=no pseudocounts.
     61          @type tau: float
     62          """
     63          from numpy import array, dot, transpose, clip
     64
     65          if self.hmm.pseudocounts or self.hmm.emission_pseudocounts:
     66              return
     67          if abs(tau) < 1e-6:
     68              return
     69
     70          # Assume probabilities
     71          if not self.hmm.score_units == ScoreUnits.Probability:
     72              self.hmm.convert_scores(units=ScoreUnits.Probability)
     73
     74          alphabet = csb.core.Enum.values(sequence.StdProteinAlphabet)
     75
     76          ## S = SubstitutionMatrix(substitution_matrix)
     77          s_mat = array(GONNET)
     78          #Normalize
     79          s_mat /= s_mat.sum()
     80          s_mat = s_mat.reshape((len(alphabet), len(alphabet)))
     81          # Marginalize matrix
     82          s_marginal = s_mat.sum(-1)
     83          s_conditional = s_mat / s_marginal
     84          # Get data and info from hmm
     85          em = array([ [layer[States.Match].emission[aa] or 0.0 for aa in alphabet]
     86                      for layer in self.hmm.layers])
     87
     88          em = clip(em, sys.float_info.min, 1.)
     89
     90          neff_m = array([l.effective_matches for l in self.hmm.layers])
     91
     92          g = dot(em, transpose(s_conditional))
     93
     94          if neff_m is not None:
     95              tau = clip(pca / (1. + (neff_m / pcb) ** pcc), 0.0, pcc)
     96              e = transpose((1. - tau) * transpose(em) + tau * transpose(g))
     97          else:
     98              e = (1. - tau) * em + tau * g
     99
    100          # Renormalize e
    101          e = transpose(transpose(e) / e.sum(-1))
    102
    103          for i, layer in enumerate(self.hmm.layers):
    104              layer[States.Match].emission.set(dict(zip(alphabet, e[i])))
    105
    106          self.hmm.emission_pseudocounts = True
    107          return
    108
    109
    110
    111 -    def add_transition_pseudocounts(self, gapb=1., gapd=0.15, gape=1.0, gapf=0.6, gapg=0.6, gapi=0.6):
    112          """
    113          Add pseudocounts to the transitions. A port from hhsearch
    114          -gapb 1.0 -gapd 0.15 -gape 1.0 -gapf 0.6 -gapg 0.6 -gapi 0.6
    115          """
    116
    117          from numpy import array
    118
    119          if not self.hmm._score_units == ScoreUnits.Probability:
    120              self.hmm.convert_scores(units=ScoreUnits.Probability)
    121
    122          if self.hmm.pseudocounts or self.hmm.transition_pseudocounts:
    123              return
    124
    125          # We need a fully populated HMM so first add all missing states
    126          states = [States.Match, States.Insertion, States.Deletion]
    127          background = self.hmm.layers[1][States.Match].background
    128          for layer in self.hmm.layers:
    129              rank = layer.rank
    130              for state in states:
    131                  if state not in layer:
    132
    133                      if state is States.Deletion:
    134                          # Add a new Deletion state
    135                          deletion = State(States.Deletion)
    136                          deletion.rank = rank
    137                          layer.append(deletion)
    138
    139                      elif state is States.Insertion:
    140                          # Add a new Deletion state
    141                          insertion = State(States.Insertion,
    142                                            emit=csb.core.Enum.members(
    143                                            sequence.SequenceAlphabets.Protein))
    144                          insertion.background.set(background)
    145                          insertion.emission.set(background)
    146                          insertion.rank = rank
    147                          layer.append(insertion)
    148
    149          if not self.hmm.start_insertion:
    150              insertion = State(States.Insertion,
    151                                            emit=csb.core.Enum.members(
    152                                            sequence.SequenceAlphabets.Protein))
    153              insertion.background.set(background)
    154              insertion.emission.set(background)
    155              insertion.rank = 0
    156              self.hmm.start_insertion = insertion
    157
    158          # make hmm completly connected
    159          for i in range(1, self.hmm.layers.length):
    160              layer = self.hmm.layers[i]
    161              #Start with match state
    162              state = layer[States.Match]
    163              if not States.Insertion in state.transitions:
    164                  state.transitions.append(Transition(state,
    165                                                      self.hmm.layers[i][States.Insertion],
    166                                                      0.0))
    167              if not States.Deletion in state.transitions:
    168                  state.transitions.append(Transition(state,
    169                                                      self.hmm.layers[i + 1][States.Deletion],
    170                                                      0.0))
    171              state = layer[States.Insertion]
    172              if not States.Insertion in state.transitions:
    173                  state.transitions.append(Transition(state,
    174                                                      self.hmm.layers[i][States.Insertion],
    175                                                      0.0))
    176              if not States.Match in state.transitions:
    177                  state.transitions.append(Transition(state,
    178                                                      self.hmm.layers[i + 1][States.Match],
    179                                                      0.0))
    180              state = layer[States.Deletion]
    181              if not States.Deletion in state.transitions:
    182                  state.transitions.append(Transition(state,
    183                                                      self.hmm.layers[i + 1][States.Deletion],
    184                                                      0.0))
    185              if not States.Match in state.transitions:
    186                  state.transitions.append(Transition(state,
    187                                                      self.hmm.layers[i + 1][States.Match],
    188                                                      0.0))
    189          # start layer
    190          state = self.hmm.start
    191          if not States.Insertion in self.hmm.start.transitions:
    192              state.transitions.append(Transition(self.hmm.start,
    193                                                  self.hmm.start_insertion,
    194                                                  0.0))
    195          if not States.Deletion in self.hmm.start.transitions:
    196              state.transitions.append(Transition(self.hmm.start,
    197                                                  self.hmm.layers[1][States.Deletion],
    198                                                  0.0))
    199
    200          state = self.hmm.start_insertion
    201          if not States.Insertion in self.hmm.start_insertion.transitions:
    202              state.transitions.append(Transition(self.hmm.start_insertion,
    203                                                  self.hmm.start_insertion,
    204                                                  0.0))
    205          if not States.Match in self.hmm.start_insertion.transitions:
    206              state.transitions.append(Transition(self.hmm.start_insertion,
    207                                                  self.hmm.layers[1][States.Match],
    208                                                  0.0))
    209
    210          # last layer
    211          state = self.hmm.layers[-1][States.Match]
    212          if not States.Insertion in state.transitions:
    213              state.transitions.append(Transition(state,
    214                                                  self.hmm.layers[-1][States.Insertion],
    215                                                  0.0))
    216          state = self.hmm.layers[-1][States.Insertion]
    217          if not States.Insertion in state.transitions:
    218              state.transitions.append(Transition(state,
    219                                                  self.hmm.layers[-1][States.Insertion],
    220                                                  0.0))
    221
    222          if not States.End in state.transitions:
    223              state.transitions.append(Transition(state,
    224                                                  self.hmm.end,
    225                                                  0.0))
    226          state = self.hmm.layers[-1][States.Deletion]
    227          if not States.End in state.transitions:
    228              state.transitions.append(Transition(state,
    229                                                  self.hmm.end,
    230                                                  0.0))
    231
    232
    233
    234          # Now we have created a fully connected HMM
    235          # Lates add pseuod counts
    236          # Calculate pseudo counts
    237
    238          # to be honest I really do not know how they came up with this
    239          pc_MD = pc_MI = 0.0286 * gapd
    240          pc_MM = 1. - 2 * pc_MD
    241          pc_DD = pc_II = gape / (gape - 1 + 1 / 0.75)
    242          pc_DM = pc_IM = 1. - pc_II
    243
    244
    245          # Get current transtion probabilities
    246          t_mm = self.hmm.start.transitions[States.Match].probability
    247          t_mi = self.hmm.start.transitions[States.Insertion].probability
    248          t_md = self.hmm.start.transitions[States.Deletion].probability
    249
    250          # Transitions from Match state
    251          n_eff = self.hmm.effective_matches
    252
    253          t = array([(n_eff - 1) * t_mm + gapb * pc_MM,
    254                     (n_eff - 1) * t_mi + gapb * pc_MI,
    255                     (n_eff - 1) * t_md + gapb * pc_MD])
    256          # normalize to one
    257          t /= t.sum()
    258          # Set
    259          self.hmm.start.transitions[States.Match].probability = t[0]
    260          self.hmm.start.transitions[States.Insertion].probability = t[1]
    261          self.hmm.start.transitions[States.Deletion].probability = t[2]
    262
    263          # Rinse and repeat
    264          t_im = self.hmm.start_insertion.transitions[States.Match].probability
    265          t_ii = self.hmm.start_insertion.transitions[States.Insertion].probability
    266
    267          t = array([t_im + gapb * pc_IM, t_ii + gapb * pc_II])
    268          t /= t.sum()
    269
    270          self.hmm.start_insertion.transitions[States.Match].probability = t[0]
    271          t_ii = self.hmm.start_insertion.transitions[States.Insertion].probability = t[1]
    272
    273          # And now for all layers
    274          for layer in self.hmm.layers[:-1]:
    275              # Get current transtion probabilities
    276              t_mm = layer[States.Match].transitions[States.Match].probability
    277              t_mi = layer[States.Match].transitions[States.Insertion].probability
    278              t_md = layer[States.Match].transitions[States.Deletion].probability
    279              n_eff = layer.effective_matches
    280              t = array([(n_eff - 1) * t_mm + gapb * pc_MM,
    281                         (n_eff - 1) * t_mi + gapb * pc_MI,
    282                         (n_eff - 1) * t_md + gapb * pc_MD])
    283              # normalize to one
    284              t /= t.sum()
    285              layer[States.Match].transitions[States.Match].probability = t[0]
    286              layer[States.Match].transitions[States.Insertion].probability = t[1]
    287              layer[States.Match].transitions[States.Deletion].probability = t[2]
    288
    289              # Transitions from insert state
    290              t_im = layer[States.Insertion].transitions[States.Match].probability
    291              t_ii = layer[States.Insertion].transitions[States.Insertion].probability
    292              n_eff = layer.effective_insertions
    293              t = array([t_im * n_eff + gapb * pc_IM,
    294                         t_im * n_eff + gapb * pc_II])
    295              # normalize to one
    296              t /= t.sum()
    297              layer[States.Insertion].transitions[States.Match].probability = t[0]
    298              layer[States.Insertion].transitions[States.Insertion].probability = t[1]
    299
    300              # Transitions form deletion state
    301              t_dm = layer[States.Deletion].transitions[States.Match].probability
    302              t_dd = layer[States.Deletion].transitions[States.Deletion].probability
    303              n_eff = layer.effective_deletions
    304              t = array([t_dm * n_eff + gapb * pc_DM,
    305                         t_dd * n_eff + gapb * pc_DD])
    306              # normalize to one
    307              t /= t.sum()
    308              layer[States.Deletion].transitions[States.Match].probability = t[0]
    309              layer[States.Deletion].transitions[States.Deletion].probability = t[1]
    310
    311          #Last layer
    312
    313          layer = self.hmm.layers[-1]
    314          t_mm = layer[States.Match].transitions[States.End].probability
    315          t_mi = layer[States.Match].transitions[States.Insertion].probability
    316          n_eff = layer.effective_matches
    317          # No deletion
    318          t = array([(n_eff - 1) * t_mm + gapb * pc_MM,
    319                     (n_eff - 1) * t_mi + gapb * pc_MI])
    320          # normalize to one
    321          t /= t.sum()
    322          layer[States.Match].transitions[States.End].probability = t[0]
    323          layer[States.Match].transitions[States.Insertion].probability = t[1]
    324
    325          # Transitions from insert state
    326          t_im = layer[States.Insertion].transitions[States.End].probability
    327          t_ii = layer[States.Insertion].transitions[States.Insertion].probability
    328          n_eff = layer.effective_insertions
    329          t = array([t_im * n_eff + gapb * pc_IM,
    330                     t_im * n_eff + gapb * pc_II])
    331          # normalize to one
    332          t /= t.sum()
    333          layer[States.Insertion].transitions[States.End].probability = t[0]
    334          layer[States.Insertion].transitions[States.Insertion].probability = t[1]
    335
    336          layer[States.Deletion].transitions[States.End].probability = 1.
    337
    338          self.hmm.transition_pseudocounts = True
    339          return
    340