text: user9641: I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that "...the object of type 'FastaIterator' has no len()...". The script however worked when only a single sequence was in the fasta file.
  In essence, I need some help figuring out how to:
  
  Prepare an instances object(is this correct?) with microsatellite
  motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10
  and (AG)10;
  Parse my 10,000 long sequence file using this object(?); and
  Print an output specifying which of my sequences contain microsatellite motifs.
  
  Appreciate any help I can get on this.
  #Seq Utils to find Motifs
  >>> pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG") 
  >>> sequence = Seq(a)
  >>> results = SeqUtils.nt_search(str(a),pattern)
  >>> print(results)
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
   56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
  8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
  8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines
  
  #BrutForce to find Motifs in a fasta file containing 10,000 sequences
  >>> def BruteForce(s, t):
  ...     occurrences = []
  ...     for i in range(len(s)-len(t)+1):
  ...         match = True
  ...         for j in range(len(t)):
  ...             if s[i+j] != t[j]:
  ...                 match = False
  ...                 break
  ...         if match:
  ...             occurrences.append(i)
  ...             
  ...     print(occurrences)
  ...     
  ... 
  
  >>> t = "TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG"
  >>> len(t)
  34
  >>> s = SeqIO.parse("group23.fas", "fasta")
  
  >>> BruteForce(s,t)
  Traceback (most recent call last):
    File "<console>", line 1, in <module>
    File "<console>", line 3, in BruteForce
  TypeError: object of type 'FastaIterator' has no len()
  
  #BrutForce to find Motifs in a fasta file containing a single sequence
  >>> s = SeqIO.read("group23v.fas", "fasta")
  >>> BruteForce(s,t)
  [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]
  
  
  user3541: Bio.SeqIO.parse() returns a SeqRecord iterator. To access the sequence of this sequence record, you need to use the .seq attribute so you should update your s with s.seq.
  
  user4980: To expand on haci's answer (and since you mentioned being new to Python) you can loop over the iterator object and use the seq attribute like this:
  pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
  for record in SeqIO.parse("group23.fas", "fasta"):
      results = SeqUtils.nt_search(str(record.seq), pattern)
      print(results)
  
  (More info on iterators here.)  That will work the same whatever the number of sequences in group23.fas.
  I'd be cautious about looking for microsatellites with fixed strings.  If the sequences you're looking for are short and have imperfections among the repeats, you might miss them, depending on how many perfect repeats your search pattern has.  For example with these stub example sequences:
  >seq1
  TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
  >seq2
  ACGATGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
  >seq3
  TGTGTGTGTGTGTGTGTGTGTCTGTGTGTGTGTGTGTGTG
  
  The nt_search output gives you:
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 0]
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 4]
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG']
  
  That is, one match at position 0 of the first sequence, one at position 4 of the second, and no matches for the third (there's a C I put in the middle).  This might not be relevant for you if you're looking microsatellites much longer than your search patterns, but for shorter ones it might be a problem.
  
  user9641: Here's a solution that I came up with after help from both haci and Jesse. It involves using the following modules: os, glob, SeqIO, Seq, SeqUtils.
  import os, glob
  
  os.chdir("C:\python38\Lib\site-packages\Bio\SeqIO")
  
  from Bio import SeqIO
  
  from Bio.Seq import Seq
  
  from Bio import SeqUtils
  
  patttern = Seq("ACACACACACACACACACAC")
  
  folder_path = 'C:\python38\Lib\site-packages\Bio\SeqIO' #specified the folder where my sequences were stored
  
  fasta_paths = glob.glob(os.path.join(folder_path, '*.fasta')) #glob module finds all the pathnames matching a specified pattern i.e. '*.fasta'
  
  with open("motif4.txt", "w") as f:
  ... for fasta_path in fasta_paths:
  ...     print(fasta_paths)
  ...     for seq_record in SeqIO.parse(fasta_path, "fasta"):
  ...         results = SeqUtils.nt_search(str(seq_record.seq),pattern)
  ...         f.write(str(results)) #replace with seq_record.seq to save the list of node names of all sequences that were parsed
  ...         f.write("\n") #ensures results of the parse through all the sequences in the fasta_path is written to a new line
  ... f.close()
  
  The nt_search output gave the following output (similar to what Jesse showed):
  1 ['ACACACACACACACACACAC']
  2 ['ACACACACACACACACACAC']
  3 ['ACACACACACACACACACAC']
  4 ['ACACACACACACACACACAC']
  5 ['ACACACACACACACACACAC']
  6 ['ACACACACACACACACACAC', 16636, 16638, 16640, 16642, 16644, 16646, 16648, 16650, 16652, 16654, 16656, 16658, 16660, 16662, 16664]
  7 ['ACACACACACACACACACAC']
  8 ['ACACACACACACACACACAC']
  9 ['ACACACACACACACACACAC']
  10 ['ACACACACACACACACACAC']
  11 ['ACACACACACACACACACAC', 9174]
  
  It would miss sequences that did not match (AC)10 or those with imperfect repeat motifs. However, in cases where the motif was > (AC)10, nt_search listed down all locations where the motif was found sequentially. For example, in line 6, all locations were listed and showed the motif was in fact 36 nucleotides long in total (pos #16,636 to #16,673).
  You can then import the txt. files into a databasing app and analyse the results.
  
original_text: 9641: <p>I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that &quot;...the object of type 'FastaIterator' has no len()...&quot;. The script however worked when only a single sequence was in the fasta file.</p>
  <p>In essence, I need some help figuring out how to:</p>
  <ul>
  <li>Prepare an instances object(is this correct?) with microsatellite
  motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10
  and (AG)10;</li>
  <li>Parse my 10,000 long sequence file using this object(?); and</li>
  <li>Print an output specifying which of my sequences contain microsatellite motifs.</li>
  </ul>
  <p>Appreciate any help I can get on this.</p>
  <p>#Seq Utils to find Motifs</p>
  <pre><code>&gt;&gt;&gt; pattern = Seq(&quot;TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG&quot;) 
  &gt;&gt;&gt; sequence = Seq(a)
  &gt;&gt;&gt; results = SeqUtils.nt_search(str(a),pattern)
  &gt;&gt;&gt; print(results)
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
   56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
  8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
  8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines
  </code></pre>
  <p>#BrutForce to find Motifs in a fasta file containing 10,000 sequences</p>
  <pre><code>&gt;&gt;&gt; def BruteForce(s, t):
  ...     occurrences = []
  ...     for i in range(len(s)-len(t)+1):
  ...         match = True
  ...         for j in range(len(t)):
  ...             if s[i+j] != t[j]:
  ...                 match = False
  ...                 break
  ...         if match:
  ...             occurrences.append(i)
  ...             
  ...     print(occurrences)
  ...     
  ... 
  
  &gt;&gt;&gt; t = &quot;TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG&quot;
  &gt;&gt;&gt; len(t)
  34
  &gt;&gt;&gt; s = SeqIO.parse(&quot;group23.fas&quot;, &quot;fasta&quot;)
  
  &gt;&gt;&gt; BruteForce(s,t)
  Traceback (most recent call last):
    File &quot;&lt;console&gt;&quot;, line 1, in &lt;module&gt;
    File &quot;&lt;console&gt;&quot;, line 3, in BruteForce
  TypeError: object of type 'FastaIterator' has no len()
  </code></pre>
  <p>#BrutForce to find Motifs in a fasta file containing a single sequence</p>
  <pre><code>&gt;&gt;&gt; s = SeqIO.read(&quot;group23v.fas&quot;, &quot;fasta&quot;)
  &gt;&gt;&gt; BruteForce(s,t)
  [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]
  </code></pre>
   3541: <p><code>Bio.SeqIO.parse()</code> returns a <code>SeqRecord iterator</code>. To access the sequence of this sequence record, you need to use the <code>.seq</code> attribute so you should update your <code>s</code> with <code>s.seq</code>.</p>
   4980: <p>To expand on haci's answer (and since you mentioned being new to Python) you can loop over the iterator object and use the <code>seq</code> attribute like this:</p>
  <pre class="lang-python prettyprint-override"><code>pattern = Seq(&quot;TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG&quot;)
  for record in SeqIO.parse(&quot;group23.fas&quot;, &quot;fasta&quot;):
      results = SeqUtils.nt_search(str(record.seq), pattern)
      print(results)
  </code></pre>
  <p>(More info on iterators <a href="https://wiki.python.org/moin/Iterator" rel="nofollow noreferrer">here</a>.)  That will work the same whatever the number of sequences in group23.fas.</p>
  <p>I'd be cautious about looking for microsatellites with fixed strings.  If the sequences you're looking for are short and have imperfections among the repeats, you might miss them, depending on how many perfect repeats your search pattern has.  For example with these stub example sequences:</p>
  <pre><code>&gt;seq1
  TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
  &gt;seq2
  ACGATGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
  &gt;seq3
  TGTGTGTGTGTGTGTGTGTGTCTGTGTGTGTGTGTGTGTG
  </code></pre>
  <p>The nt_search output gives you:</p>
  <pre><code>['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 0]
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 4]
  ['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG']
  </code></pre>
  <p>That is, one match at position 0 of the first sequence, one at position 4 of the second, and no matches for the third (there's a C I put in the middle).  This might not be relevant for you if you're looking microsatellites much longer than your search patterns, but for shorter ones it might be a problem.</p>
   9641: <p>Here's a solution that I came up with after help from both haci and Jesse. It involves using the following modules: os, glob, SeqIO, Seq, SeqUtils.</p>
  <pre><code>import os, glob
  
  os.chdir(&quot;C:\python38\Lib\site-packages\Bio\SeqIO&quot;)
  
  from Bio import SeqIO
  
  from Bio.Seq import Seq
  
  from Bio import SeqUtils
  
  patttern = Seq(&quot;ACACACACACACACACACAC&quot;)
  
  folder_path = 'C:\python38\Lib\site-packages\Bio\SeqIO' #specified the folder where my sequences were stored
  
  fasta_paths = glob.glob(os.path.join(folder_path, '*.fasta')) #glob module finds all the pathnames matching a specified pattern i.e. '*.fasta'
  
  with open(&quot;motif4.txt&quot;, &quot;w&quot;) as f:
  ... for fasta_path in fasta_paths:
  ...     print(fasta_paths)
  ...     for seq_record in SeqIO.parse(fasta_path, &quot;fasta&quot;):
  ...         results = SeqUtils.nt_search(str(seq_record.seq),pattern)
  ...         f.write(str(results)) #replace with seq_record.seq to save the list of node names of all sequences that were parsed
  ...         f.write(&quot;\n&quot;) #ensures results of the parse through all the sequences in the fasta_path is written to a new line
  ... f.close()
  </code></pre>
  <p>The nt_search output gave the following output (similar to what Jesse showed):</p>
  <pre><code>1 ['ACACACACACACACACACAC']
  2 ['ACACACACACACACACACAC']
  3 ['ACACACACACACACACACAC']
  4 ['ACACACACACACACACACAC']
  5 ['ACACACACACACACACACAC']
  6 ['ACACACACACACACACACAC', 16636, 16638, 16640, 16642, 16644, 16646, 16648, 16650, 16652, 16654, 16656, 16658, 16660, 16662, 16664]
  7 ['ACACACACACACACACACAC']
  8 ['ACACACACACACACACACAC']
  9 ['ACACACACACACACACACAC']
  10 ['ACACACACACACACACACAC']
  11 ['ACACACACACACACACACAC', 9174]
  </code></pre>
  <p>It would miss sequences that did not match (AC)10 or those with imperfect repeat motifs. However, in cases where the motif was &gt; (AC)10, nt_search listed down all locations where the motif was found sequentially. For example, in line 6, all locations were listed and showed the motif was in fact 36 nucleotides long in total (pos #16,636 to #16,673).</p>
  <p>You can then import the txt. files into a databasing app and analyse the results.</p>
  