Answers.Body: <p>Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-13T09:49:18.730Z
Answers.Text: Right I see, so I need to preprocess the sequence alignments first? So that I match the gaps properly?
Answers.UserDisplayName: space_gladiator
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-13T10:54:30.307Z
Answers.Text: Yes, to verify this step for logical errors
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-12T22:41:14.647Z
Answers.LastActivityDate: 2020-04-12T22:41:14.647Z
Answers.UserDisplayName: M__
Body: <p>I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.</p>
  
  <p>So, here is my code for the Monte Carlo method:</p>
  
  <pre><code>def MonteCarlo(S, max_iter, threshold):
      """
      Computes the discriminating position for a sequence column
  
      Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
      Output: 2 classes
      """
      n = len(S)
      S_copy = S.copy()
      num_class1 = int(0.30*n) #number in class1
      print(num_class1)
      for _ in range(num_class1): #intialisation
          rand_position = random.randint(0,len(S)-1)
          if S_copy[rand_position] == 'X':
              rand_position = random.randint(0,len(S)) #reselect if used position
          S_copy[rand_position] = 'X' 
  
      #Monte Carlo with Insertions and Deletions
  
      s_entropy = entropy(S)
      n_iter = 0 #counter for iterations
      class_1 = dict() #class 1 determined by X's
      class_2 = dict() #class 2 is rest of S
      for i in range(len(S)):
          if S_copy[i] == 'X':
              class_1[i] =S[i]
          else:
              class_2[i] =S[i]
  
      res = score(S, class_1, class_2)
  
      while n_iter &lt; max_iter:
          r = random.random()
  
          if r &lt; 0.25: #25% chance to do an insertion
            (c1,c2,pos) = Insertion(S, class_1, class_2)
  
          elif r &lt; 0.5: #25% chance to do a deletion
            (c1,c2,pos) = Deletion(S, class_1, class_2)
  
          else: # 50% chance to do a swap
            (c1,c2,pos) = Swap(S, class_1, class_2)
  
          res2 = score(S, c1, c2)
  
          #Boltzmann Distribution coefficient 
          delta_S = abs(res2 - res)
          lambda_, N = 0.2, 19
          size = 6000
          ind1 = random.randint(0,size - 1)
          rv = boltzmann.rvs(lambda_, N, size=size)/N 
          beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
          if  (res2 &gt; res) or (beta*delta_S &gt; threshold): #if better score i.e. higher entropy gain
              res = res2
              class_1 = c1
              class_2 = c2
  
          n_iter += 1
  
  
      return class_1, class_2, pos
  
  </code></pre>
  
  <p>Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)</p>
  
  <p>Let me explain what my code is trying to do: Given a sequence alignment (as a list) say ['M','M','M',...,'L','L','L'] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.</p>
  
  <p>While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results. </p>
  
  <p>Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?</p>
  
  <p>Thanks!</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-12T22:09:08.977Z
Comments.Text: It's not very clear what you mean by "discriminating position". Can you explain more or link to a source that describes in a little more detail what you are talking about?
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-13T09:48:25.540Z
Comments.Text: Hi! By discriminating position, I mean I'm trying to find a position within the sequence alignment that splits the columns into two classes such that the two have the lowest entropy possible. So, in my sample case with a clear division between Ms and Ls, I should've gotten the the position at which the Ls begin. Does that make sense?
  
  Here is a link: https://pubs.rsc.org/en/content/articlelanding/2019/sc/c8sc03669a#!divAbstract
Comments.UserDisplayName: space_gladiator
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-13T15:01:06.250Z
Comments.Text: Hi, the paper that you've linked to is talking about a discrimination between two structural states of a protein by a kinase; "discriminating position" does not appear to occur in its full text. It's not clear how this generalizes to alignments. From your explanation though I think I have a better sense of what you mean: for seq alignment $S$ with positions i = 1,2,3...n find the position $S_i$ such that a measure of entropy across aligned positions $S_1 : S_{i-1}$ and $S_{i+1} : S_n$ is minimized for both slices across each column of the aligned slice. (words are easier than code for me).
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-04-13T18:37:46.100Z
Comments.Text: Yes! no problem (: I also misunderstood what I was trying to solve actually, I have managed to fix the problem! But for future reference: the code is supposed to split the sequence column into two classes such that each distinctly has the highest entropy possible (or lowest for that matter)
Comments.UserDisplayName: space_gladiator
Tags: <python><phylogenetics><proteins>
Title: Finding discriminating positions for a sequence alignment column
ThreadText: # Finding discriminating positions for a sequence alignment column
  
  (Asked by: username_0 on 2020-04-12)
  
  I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.
  
  So, here is my code for the Monte Carlo method:
  
  ```
  def MonteCarlo(S, max_iter, threshold):
      """
      Computes the discriminating position for a sequence column
  
      Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
      Output: 2 classes
      """
      n = len(S)
      S_copy = S.copy()
      num_class1 = int(0.30*n) #number in class1
      print(num_class1)
      for _ in range(num_class1): #intialisation
          rand_position = random.randint(0,len(S)-1)
          if S_copy[rand_position] == 'X':
              rand_position = random.randint(0,len(S)) #reselect if used position
          S_copy[rand_position] = 'X' 
  
      #Monte Carlo with Insertions and Deletions
  
      s_entropy = entropy(S)
      n_iter = 0 #counter for iterations
      class_1 = dict() #class 1 determined by X's
      class_2 = dict() #class 2 is rest of S
      for i in range(len(S)):
          if S_copy[i] == 'X':
              class_1[i] =S[i]
          else:
              class_2[i] =S[i]
  
      res = score(S, class_1, class_2)
  
      while n_iter < max_iter:
          r = random.random()
  
          if r < 0.25: #25% chance to do an insertion
            (c1,c2,pos) = Insertion(S, class_1, class_2)
  
          elif r < 0.5: #25% chance to do a deletion
            (c1,c2,pos) = Deletion(S, class_1, class_2)
  
          else: # 50% chance to do a swap
            (c1,c2,pos) = Swap(S, class_1, class_2)
  
          res2 = score(S, c1, c2)
  
          #Boltzmann Distribution coefficient 
          delta_S = abs(res2 - res)
          lambda_, N = 0.2, 19
          size = 6000
          ind1 = random.randint(0,size - 1)
          rv = boltzmann.rvs(lambda_, N, size=size)/N 
          beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
          if  (res2 > res) or (beta*delta_S > threshold): #if better score i.e. higher entropy gain
              res = res2
              class_1 = c1
              class_2 = c2
  
          n_iter += 1
  
      return class_1, class_2, pos
  
  ```
  
  Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)
  
  Let me explain what my code is trying to do: Given a sequence alignment (as a list) say ['M','M','M',...,'L','L','L'] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num\_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.
  
  While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results.
  
  Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?
  
  Thanks!
  
  **Tags:** python, phylogenetics, proteins
  
  > **Comments:**
  > **username_1**: It's not very clear what you mean by "discriminating position". Can you explain more or link to a source that describes in a little more detail what you are talking about?  
  > **username_0**: Hi! By discriminating position, I mean I'm trying to find a position within the sequence alignment that splits the columns into two classes such that the two have the lowest entropy possible. So, in my sample case with a clear division between Ms and Ls, I should've gotten the the position at which the Ls begin. Does that make sense?  Here is a link: https://pubs.rsc.org/en/content/articlelanding/2019/sc/c8sc03669a#!divAbstract  
  > **username_1**: Hi, the paper that you've linked to is talking about a discrimination between two structural states of a protein by a kinase; "discriminating position" does not appear to occur in its full text. It's not clear how this generalizes to alignments. From your explanation though I think I have a better sense of what you mean: for seq alignment $S$ with positions i = 1,2,3...n find the position $S_i$ such that a measure of entropy across aligned positions $S_1 : S_{i-1}$ and $S_{i+1} : S_n$ is minimized for both slices across each column of the aligned slice. (words are easier than code for me).  
  > **username_0**: Yes! no problem (: I also misunderstood what I was trying to solve actually, I have managed to fix the problem! But for future reference: the code is supposed to split the sequence column into two classes such that each distinctly has the highest entropy possible (or lowest for that matter)  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2020-04-12. Score: 1
  
  Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.
  
  > **Comments:**
  > **username_0**: Right I see, so I need to preprocess the sequence alignments first? So that I match the gaps properly?  
  > **username_2**: Yes, to verify this step for logical errors  
  
  ---
  
  