Answers.Body: <p>I overlooked the fact that (for my data) about 1:50,000 results returned from my blasting yielded a ";" within the hit denoting two taxids and did not account for these in my script.  </p>
  
  <p>How a given researcher will respond to their own results having ";"-hits will probably vary (you may want to do extra analysis, etc.).  In my case, I am only interested in the first taxid before the ";".  So I added a simply python .split() to my script before appending the taxid to the list that I query:</p>
  
  <pre><code>if __name__ == '__main__':
      file = open(sys.argv[1], "r")    
      taxids = []
      contigs = []
      for line in file:
          line = line.split("\n")[0]
          tID = line.split(",")[0]
          if ';' in tID: 
              tID = tID.split(";")[0]
          taxids.append(tID)
          contigs.append(line.split(",")[1])
  
      desired_ranks = ['superkingdom', 'phylum']
      results = list()
      for taxid in taxids:
          results.append(list())
          results[-1].append(str(taxid))
          ranks = get_desired_ranks(taxid, desired_ranks)
          for key, rank in ranks.items():
              if rank != '&lt;not present&gt;':
                  results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
              else:
                  results[-1].append(rank)
  
      i = 0
      for result in results:
          print(contigs[i] + ','),
          print(','.join(result))
          i += 1
  
      file.close()
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-01T08:06:11.057Z
Answers.LastActivityDate: 2018-12-01T08:06:11.057Z
Answers.UserDisplayName: ljs
Body: <p>I am using this script:</p>
  
  <pre><code>import csv
  import time
  import sys
  from ete3 import NCBITaxa
  
  ncbi = NCBITaxa()
  
  def get_desired_ranks(taxid, desired_ranks):
      lineage = ncbi.get_lineage(taxid)   
      names = ncbi.get_taxid_translator(lineage)
      lineage2ranks = ncbi.get_rank(names)
      ranks2lineage = dict((rank,taxid) for (taxid, rank) in lineage2ranks.items())
      return{'{}_id'.format(rank): ranks2lineage.get(rank, '&lt;not present&gt;') for rank in desired_ranks}
  
  if __name__ == '__main__':
      file = open(sys.argv[1], "r")    
      taxids = []
      contigs = []
      for line in file:
          line = line.split("\n")[0]
          taxids.append(line.split(",")[0])
          contigs.append(line.split(",")[1])
  
      desired_ranks = ['superkingdom', 'phylum']
      results = list()
      for taxid in taxids:
          results.append(list())
          results[-1].append(str(taxid))
          ranks = get_desired_ranks(taxid, desired_ranks)
          for key, rank in ranks.items():
              if rank != '&lt;not present&gt;':
                  results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
              else:
                  results[-1].append(rank)
  
      i = 0
      for result in results:
          print(contigs[i] + ','),
          print(','.join(result))
          i += 1
  
      file.close()
  </code></pre>
  
  <p>The script takes taxids from a file and fetches their respective lineages from a local copy of NCBI's Taxonomy database.  Strangely, this script works fine when I run it on small sets of taxids (~70, ~100), but most of my datasets are upwards of 280k taxids and these break the script.</p>
  
  <p>I get this complete error:</p>
  
  <pre><code>Traceback (most recent call last):
    File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 31, in &lt;module&gt;
      ranks = get_desired_ranks(taxid, desired_ranks)
    File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 11, in get_desired_ranks
      lineage = ncbi.get_lineage(taxid)   
    File "/data1/lstout/.local/lib/python2.7/site-packages/ete3/ncbi_taxonomy/ncbiquery.py", line 227, in get_lineage
      result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
  sqlite3.Warning: You can only execute one statement at a time.
  </code></pre>
  
  <p>The first two files from the traceback are simply the script I referenced above, the third file is one of ete3's.  And as I stated, the script works fine with small datasets.</p>
  
  <p><strong>What I have tried:</strong></p>
  
  <ul>
  <li>Importing the time module and sleeping for a few milliseconds/hundredths of a second before/after my offending lines of code on lines 11 and 31.  No effect.</li>
  <li><p>Went to line 227 in ete3's code...</p>
  
  <pre><code>result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
  </code></pre></li>
  </ul>
  
  <p>and changed the "execute" function to "executescript" in order to be able to handle multiple queries at once (as that seems to be the problem).  This produced a new error and led to a rabbit hole of me changing minor things in their script trying to fudge this to work.  No result.  This is the complete offending function:</p>
  
  <pre><code>    def get_lineage(self, taxid):
      """Given a valid taxid number, return its corresponding lineage track as a
      hierarchically sorted list of parent taxids.
      """
      if not taxid:
          return None
      result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
      raw_track = result.fetchone()
      if not raw_track:
          #perhaps is an obsolete taxid
          _, merged_conversion = self._translate_merged([taxid])
          if taxid in merged_conversion:
              result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
              raw_track = result.fetchone()
          # if not raise error
          if not raw_track:
              #raw_track = ["1"]
              raise ValueError("%s taxid not found" %taxid)
          else:
              warnings.warn("taxid %s was translated into %s" %(taxid, merged_conversion[taxid]))
  
      track = list(map(int, raw_track[0].split(",")))
      return list(reversed(track))
  </code></pre>
  
  <p>What bothers me so much is that this works on small amounts of data!  I'm running these scripts from my school's high performance computer and have tried running on their head node and in an interactive moab scheduler.  Nothing has helped.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-29T11:36:03.120Z
Comments.Text: Please [edit] your question and include the script you are using. Both so we don't need to go back and forth between different tabs, and because the linked script may be edited. Also so we can see your modifications, understand what lines the error messages refer to etc.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-29T12:28:17.297Z
Comments.Text: Since you have 280k taxids i think querying is not the way to go; instead download a flat file from NCBI taxonomy with all the data and process that
Comments.UserDisplayName: Chris_Rands
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-11-29T15:25:00.833Z
Comments.Text: @Chris_Rands I am using a local copy of taxonomy,  ncbi=ncbi.Taxa() gets the DB.
Comments.UserDisplayName: ljs
Tags: <python><phylogenetics><database><taxonomy><ete3>
Title: Python scripting with ete3 to query NCBI's Taxonomy: "sqlite3 Warning (can only execute one statement at a time)"
ThreadText: # Python scripting with ete3 to query NCBI's Taxonomy: "sqlite3 Warning (can only execute one statement at a time)"
  
  (Asked by: username_0 on 2018-11-29)
  
  I am using this script:
  
  ```
  import csv
  import time
  import sys
  from ete3 import NCBITaxa
  
  ncbi = NCBITaxa()
  
  def get_desired_ranks(taxid, desired_ranks):
      lineage = ncbi.get_lineage(taxid)   
      names = ncbi.get_taxid_translator(lineage)
      lineage2ranks = ncbi.get_rank(names)
      ranks2lineage = dict((rank,taxid) for (taxid, rank) in lineage2ranks.items())
      return{'{}_id'.format(rank): ranks2lineage.get(rank, '<not present>') for rank in desired_ranks}
  
  if __name__ == '__main__':
      file = open(sys.argv[1], "r")    
      taxids = []
      contigs = []
      for line in file:
          line = line.split("\n")[0]
          taxids.append(line.split(",")[0])
          contigs.append(line.split(",")[1])
  
      desired_ranks = ['superkingdom', 'phylum']
      results = list()
      for taxid in taxids:
          results.append(list())
          results[-1].append(str(taxid))
          ranks = get_desired_ranks(taxid, desired_ranks)
          for key, rank in ranks.items():
              if rank != '<not present>':
                  results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
              else:
                  results[-1].append(rank)
  
      i = 0
      for result in results:
          print(contigs[i] + ','),
          print(','.join(result))
          i += 1
  
      file.close()
  
  ```
  
  The script takes taxids from a file and fetches their respective lineages from a local copy of NCBI's Taxonomy database. Strangely, this script works fine when I run it on small sets of taxids (~70, ~100), but most of my datasets are upwards of 280k taxids and these break the script.
  
  I get this complete error:
  
  ```
  Traceback (most recent call last):
    File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 31, in <module>
      ranks = get_desired_ranks(taxid, desired_ranks)
    File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 11, in get_desired_ranks
      lineage = ncbi.get_lineage(taxid)   
    File "/data1/lstout/.local/lib/python2.7/site-packages/ete3/ncbi_taxonomy/ncbiquery.py", line 227, in get_lineage
      result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
  sqlite3.Warning: You can only execute one statement at a time.
  
  ```
  
  The first two files from the traceback are simply the script I referenced above, the third file is one of ete3's. And as I stated, the script works fine with small datasets.
  
  **What I have tried:**
  
  * Importing the time module and sleeping for a few milliseconds/hundredths of a second before/after my offending lines of code on lines 11 and 31. No effect.
  * Went to line 227 in ete3's code...
  
    ```
    result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
  
    ```
  
  and changed the "execute" function to "executescript" in order to be able to handle multiple queries at once (as that seems to be the problem). This produced a new error and led to a rabbit hole of me changing minor things in their script trying to fudge this to work. No result. This is the complete offending function:
  
  ```
      def get_lineage(self, taxid):
      """Given a valid taxid number, return its corresponding lineage track as a
      hierarchically sorted list of parent taxids.
      """
      if not taxid:
          return None
      result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
      raw_track = result.fetchone()
      if not raw_track:
          #perhaps is an obsolete taxid
          _, merged_conversion = self._translate_merged([taxid])
          if taxid in merged_conversion:
              result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
              raw_track = result.fetchone()
          # if not raise error
          if not raw_track:
              #raw_track = ["1"]
              raise ValueError("%s taxid not found" %taxid)
          else:
              warnings.warn("taxid %s was translated into %s" %(taxid, merged_conversion[taxid]))
  
      track = list(map(int, raw_track[0].split(",")))
      return list(reversed(track))
  
  ```
  
  What bothers me so much is that this works on small amounts of data! I'm running these scripts from my school's high performance computer and have tried running on their head node and in an interactive moab scheduler. Nothing has helped.
  
  **Tags:** python, phylogenetics, database, taxonomy, ete3
  
  > **Comments:**
  > **username_1**: Please [edit] your question and include the script you are using. Both so we don't need to go back and forth between different tabs, and because the linked script may be edited. Also so we can see your modifications, understand what lines the error messages refer to etc.  
  > **username_2**: Since you have 280k taxids i think querying is not the way to go; instead download a flat file from NCBI taxonomy with all the data and process that  
  > **username_0**: @username_2 I am using a local copy of taxonomy,  ncbi=ncbi.Taxa() gets the DB.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2018-12-01. Score: 1 | Accepted answer
  
  I overlooked the fact that (for my data) about 1:50,000 results returned from my blasting yielded a ";" within the hit denoting two taxids and did not account for these in my script.
  
  How a given researcher will respond to their own results having ";"-hits will probably vary (you may want to do extra analysis, etc.). In my case, I am only interested in the first taxid before the ";". So I added a simply python .split() to my script before appending the taxid to the list that I query:
  
  ```
  if __name__ == '__main__':
      file = open(sys.argv[1], "r")    
      taxids = []
      contigs = []
      for line in file:
          line = line.split("\n")[0]
          tID = line.split(",")[0]
          if ';' in tID: 
              tID = tID.split(";")[0]
          taxids.append(tID)
          contigs.append(line.split(",")[1])
  
      desired_ranks = ['superkingdom', 'phylum']
      results = list()
      for taxid in taxids:
          results.append(list())
          results[-1].append(str(taxid))
          ranks = get_desired_ranks(taxid, desired_ranks)
          for key, rank in ranks.items():
              if rank != '<not present>':
                  results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
              else:
                  results[-1].append(rank)
  
      i = 0
      for result in results:
          print(contigs[i] + ','),
          print(','.join(result))
          i += 1
  
      file.close()
  
  ```
  
  
  ---
  
  