Body: <p>I want to map TFBS (transcription factor binding sites) on to the sequences from dataset containing 8448 (210 base) sequences, but it is  not giving all the hits I need</p>
  <pre><code>import pickle
  from pickle import load, dump
  from collections import defaultdict
  from Bio.Seq import Seq
  
  #first element is BREu = SSRCGCC, length = 7 bases, size of matrix = 4*7, max score = 3.33
  # first row = A, second row = T, third row = G and forth row = C
  BREu = [
      [0.0, 0.0, 0.5, 0, 0, 0, 0],
      [0.0, 0.0, 0.0, 0, 0, 0, 0],
      [0.5, 0.5, 0.5, 0, 1, 0, 0],
      [0.5, 0.5, 0.0, 1, 0, 1, 1]]
  
  #second element is ETS = MGGAART, length = 7 bases, size of matrix = 4*7, max score = 6
  ETS = [
      [0.5, 0, 0, 1, 1, 0.5, 0],
      [0.0, 0, 0, 0, 0, 0.0, 1],
      [0.0, 1, 1, 0, 0, 0.5, 0],
      [0.5, 0, 0, 0, 0, 0.0, 0]]
  
  #Third element is INR = YYANWYY, length = 7 bases, size of matrix = 4*7, max score = 3.75
  INR = [
      [0.0, 0.0, 1, 0.25, 0.5, 0.0, 0.0],
      [0.5, 0.5, 0, 0.25, 0.5, 0.5, 0.5],
      [0.0, 0.0, 0, 0.25, 0.0, 0.0, 0.0],
      [0.5, 0.5, 0, 0.25, 0.0, 0.5, 0.5]]
  
  #forth element is NFY = CCAAT, length = 5 bases, size of matrix = 4*5, max score = 5
  NFY = [
      [0, 0, 1, 1, 0],
      [0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0],
      [1, 1, 0, 0, 0]]
  
  #seventh element is SP1_NC = GGGNGGG, length = 7 bases, size of matrix = 4*7, max score = 7
  #Gives some possibilities of having any bases at position 4 while still favoring G
  #threshold score &gt;= 6.5 : keep all degenerated motif with degenerated = 6.5 and canonical = 7
  
  SP1 = [
      [0, 0, 0, .5, 0, 0, 0],
      [0, 0, 0, .5, 0, 0, 0],
      [1, 1, 1, .5, 1, 1, 1],
      [0, 0, 0, 1., 0, 0, 0]]
  
  #sixth element is TATA = TATAWAAG, length = 8 bases, size of matrix = 4*8, max score = 7.5
  #threshold score &gt; 20 : keep all degenerated motif with degenerated = 20 and canonical = 23
  TATA = [
      [0, 4, 0, 4, 4, 1, 1, 0],
      [4, 0, 4, 0, 4, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0]]
  
  # define matrix and thresholds as tuple
  motifs = dict(
      BREu = (BREu, 5.5),
      ETS = (ETS, 6),
      INR = (INR, 3.75),
      NFY = (NFY, 5),
      SP1 = (SP1, 6.25),
      TATA = (TATA, 20))
  
  def get_motif_score(seq, motif):
  
      motif_size = len(motif[0])
      bases_encodings = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
      coded_seq = [
          bases_encodings[b] if b in 'ATGC' else None
          for b in seq]
  
      return [
          sum(
              0 if coded_b is None else motif[coded_b][i]
              for i, coded_b in enumerate(coded_seq[start:start+motif_size])
          ) for start in range(len(seq) - motif_size)]
  
  
  def get_motifs(seq):
      &quot;&quot;&quot;find motif occurences in given sequence
      
      return list of (motif sequence, start, end, score)
  
      &quot;&quot;&quot;
      seq_size = len(seq)
      motif_hits = {}
  
      for motif_name, (motif_matrix, threshold) in motifs.items():
          motif_size = len(motif_matrix[0])
          rcseq = str(Seq(seq).reverse_complement())
          scores = get_motif_score(seq, motif_matrix)
          rc_scores = get_motif_score(rcseq, motif_matrix)
          
          motif_hits[motif_name] = [
              (seq[i:i+motif_size], i, i+motif_size, score,'+')
              for i, score in enumerate(scores)
              if score &gt;= threshold
          ] + [
              (rcseq[i:i+motif_size], seq_size-(i+motif_size), (seq_size-i), score,'-')
              for i, score in enumerate(rc_scores)
              if score &gt;= threshold
          ]
  
      return motif_hits
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-08T17:36:08.913Z
Comments.Text: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T13:17:35.363Z
Comments.Text: What is the use and how we should interpret the matrix in your code ???
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T13:25:32.217Z
Comments.Text: Hi pippo1930, The matrices are self designed to map the motifs on basis of arbitary score, with the function of getting score, I am calculating the score on the window size equal to the length of the motif (length of matrix row), then if the score is above the threshold as indicated, it means that the motif exist. note that the sequence are non - consensus, which means that some basis are not very particular (that is why the score to map that particular base is less). I am not sure if that answer your question and thank you for taking interest to solve this problem.
Comments.UserDisplayName: Bioinfo_nerd
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T13:28:29.803Z
Comments.Text: Its a really hard alghiritm to understand, here usually you should provide a test input with the expected results vs actual ones to help people understand the issue
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T13:34:36.460Z
Comments.Text: How of curiosity could you show me how a BREu motif can get a > 5,5 score
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T14:04:29.990Z
Comments.Text: print(get_motifs('TATAWAAG')) gives {'BREu': [], 'ETS': [], 'INR': [], 'NFY': [], 'SP1': [], 'TATA': []}
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T14:52:40.327Z
Comments.Text: lower the threshold ???
Comments.UserDisplayName: pippo1980
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-09-09T14:56:55.110Z
Comments.Text: or change the matrix values
Comments.UserDisplayName: pippo1980
Tags: <python><phylogenetics>
Title: I want to map TFBS on to the sequences from dataset containing 8448 (210 base) sequences. its not giving all the hits I need
ThreadText: # I want to map TFBS on to the sequences from dataset containing 8448 (210 base) sequences. its not giving all the hits I need
  
  (Asked by: username_0 on 2022-09-08)
  
  I want to map TFBS (transcription factor binding sites) on to the sequences from dataset containing 8448 (210 base) sequences, but it is not giving all the hits I need
  
  ```
  import pickle
  from pickle import load, dump
  from collections import defaultdict
  from Bio.Seq import Seq
  
  #first element is BREu = SSRCGCC, length = 7 bases, size of matrix = 4*7, max score = 3.33
  # first row = A, second row = T, third row = G and forth row = C
  BREu = [
      [0.0, 0.0, 0.5, 0, 0, 0, 0],
      [0.0, 0.0, 0.0, 0, 0, 0, 0],
      [0.5, 0.5, 0.5, 0, 1, 0, 0],
      [0.5, 0.5, 0.0, 1, 0, 1, 1]]
  
  #second element is ETS = MGGAART, length = 7 bases, size of matrix = 4*7, max score = 6
  ETS = [
      [0.5, 0, 0, 1, 1, 0.5, 0],
      [0.0, 0, 0, 0, 0, 0.0, 1],
      [0.0, 1, 1, 0, 0, 0.5, 0],
      [0.5, 0, 0, 0, 0, 0.0, 0]]
  
  #Third element is INR = YYANWYY, length = 7 bases, size of matrix = 4*7, max score = 3.75
  INR = [
      [0.0, 0.0, 1, 0.25, 0.5, 0.0, 0.0],
      [0.5, 0.5, 0, 0.25, 0.5, 0.5, 0.5],
      [0.0, 0.0, 0, 0.25, 0.0, 0.0, 0.0],
      [0.5, 0.5, 0, 0.25, 0.0, 0.5, 0.5]]
  
  #forth element is NFY = CCAAT, length = 5 bases, size of matrix = 4*5, max score = 5
  NFY = [
      [0, 0, 1, 1, 0],
      [0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0],
      [1, 1, 0, 0, 0]]
  
  #seventh element is SP1_NC = GGGNGGG, length = 7 bases, size of matrix = 4*7, max score = 7
  #Gives some possibilities of having any bases at position 4 while still favoring G
  #threshold score >= 6.5 : keep all degenerated motif with degenerated = 6.5 and canonical = 7
  
  SP1 = [
      [0, 0, 0, .5, 0, 0, 0],
      [0, 0, 0, .5, 0, 0, 0],
      [1, 1, 1, .5, 1, 1, 1],
      [0, 0, 0, 1., 0, 0, 0]]
  
  #sixth element is TATA = TATAWAAG, length = 8 bases, size of matrix = 4*8, max score = 7.5
  #threshold score > 20 : keep all degenerated motif with degenerated = 20 and canonical = 23
  TATA = [
      [0, 4, 0, 4, 4, 1, 1, 0],
      [4, 0, 4, 0, 4, 0, 0, 0],
      [0, 0, 0, 0, 0, 0, 0, 1],
      [0, 0, 0, 0, 0, 0, 0, 0]]
  
  # define matrix and thresholds as tuple
  motifs = dict(
      BREu = (BREu, 5.5),
      ETS = (ETS, 6),
      INR = (INR, 3.75),
      NFY = (NFY, 5),
      SP1 = (SP1, 6.25),
      TATA = (TATA, 20))
  
  def get_motif_score(seq, motif):
  
      motif_size = len(motif[0])
      bases_encodings = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
      coded_seq = [
          bases_encodings[b] if b in 'ATGC' else None
          for b in seq]
  
      return [
          sum(
              0 if coded_b is None else motif[coded_b][i]
              for i, coded_b in enumerate(coded_seq[start:start+motif_size])
          ) for start in range(len(seq) - motif_size)]
  
  def get_motifs(seq):
      """find motif occurences in given sequence
      
      return list of (motif sequence, start, end, score)
  
      """
      seq_size = len(seq)
      motif_hits = {}
  
      for motif_name, (motif_matrix, threshold) in motifs.items():
          motif_size = len(motif_matrix[0])
          rcseq = str(Seq(seq).reverse_complement())
          scores = get_motif_score(seq, motif_matrix)
          rc_scores = get_motif_score(rcseq, motif_matrix)
          
          motif_hits[motif_name] = [
              (seq[i:i+motif_size], i, i+motif_size, score,'+')
              for i, score in enumerate(scores)
              if score >= threshold
          ] + [
              (rcseq[i:i+motif_size], seq_size-(i+motif_size), (seq_size-i), score,'-')
              for i, score in enumerate(rc_scores)
              if score >= threshold
          ]
  
      return motif_hits
  
  ```
  
  **Tags:** python, phylogenetics
  
  > **Comments:**
  > **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
  > **username_2**: What is the use and how we should interpret the matrix in your code ???  
  > **username_0**: Hi pippo1930, The matrices are self designed to map the motifs on basis of arbitary score, with the function of getting score, I am calculating the score on the window size equal to the length of the motif (length of matrix row), then if the score is above the threshold as indicated, it means that the motif exist. note that the sequence are non - consensus, which means that some basis are not very particular (that is why the score to map that particular base is less). I am not sure if that answer your question and thank you for taking interest to solve this problem.  
  > **username_2**: Its a really hard alghiritm to understand, here usually you should provide a test input with the expected results vs actual ones to help people understand the issue  
  > **username_2**: How of curiosity could you show me how a BREu motif can get a > 5,5 score  
  > **username_2**: print(get_motifs('TATAWAAG')) gives {'BREu': [], 'ETS': [], 'INR': [], 'NFY': [], 'SP1': [], 'TATA': []}  
  > **username_2**: lower the threshold ???  
  > **username_2**: or change the matrix values  
  
  ---
  
  ## 0 answers
  
  