Answers.Body: <p>Yes its easy using Biopython.</p>
  <p>What I'm requesting is the sequence must be &gt;9500 otherwise it's discarded. In your case this value will be MUCH bigger, because these are bacterial genomes.</p>
  <pre><code>from Bio import Entrez
  from Bio import SeqIO, AlignIO
  
  def getgenbank (smallGenDict):
      bacteriadb = []
      for bacteria in smallGenDict.keys():
          Entrez.email = &quot;name@email.com&quot;
          with Entrez.efetch(
                  db=&quot;nucleotide&quot;, rettype=&quot;gb&quot;, retmode=&quot;text&quot;, id=bacteria
                  ) as handle:
              seq_record = SeqIO.read(handle, &quot;gb&quot;)
              # record = get_CDS(seq_record, bacteria, smallGenDict.keys)
              record = seq_record # better to remove the above function because its outside the question
              if len(record.seq) &gt; 9500:
                  bacteriadb.append(record)
      return bacteriadb
  
  def get_CDS(record, bacteria, smallGenDict):
      for feature in record.features:
          # custom manipulation of input sequences 
          # better to forget this bit of code
          # this whole part of the code has been skipped
          # its capturing the coding sequence (not relevant)
  return record
  
  def main(smallGenDict):
      db = getgenbank(smallGenDict)
      SeqIO.write(db, '/Users/username/Desktop/myoutput.fasta', &quot;fasta&quot;)
  
  if __name__ == '__main__':
      smallGenDict = {'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'} # better to automate this bit (different question)
      main(smallGenDict)
  </code></pre>
  <p>Thats it.</p>
  <p><strong>Notes</strong> In the code below, the file smallGenDict is a dictionary of Genbank accession No.s and their sequence.
  <code>{'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'}</code>
  This better written here as a list
  <code>['ABxxxxxx','ABxxxxx2' ...]</code> say call it <code>smallGenList</code>.</p>
  <p>Where it says <code>smallGenDict.keys():</code> just replace it with <code>smallGenList</code> and remove the method <code>.keys()</code>, otherwise the code will through a bug.
  The function <code>get_CDS</code> is processing the sequence file to remove the non-coding regions. I've not included it because its specific stuff.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-11T16:14:19.860Z
Answers.LastActivityDate: 2023-05-11T16:19:50.040Z
Answers.LastEditDate: 2023-05-11T16:19:50.040Z
Answers.UserDisplayName: M__
Answers.Body: <p>You can use R <a href="https://cran.r-project.org/web/packages/rentrez/vignettes/rentrez_tutorial.html" rel="nofollow noreferrer">Rentrez</a>, in the link you can find a tutorial for the package usage.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-11T16:03:31.080Z
Answers.Text: Thank you! I've tried Entrez (R is not something I am familiar with), but I got stuck.
Answers.UserDisplayName: Laura
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-05-11T15:35:29.673Z
Answers.LastActivityDate: 2023-05-11T15:35:29.673Z
Answers.UserDisplayName: flavinsky
Body: <p>It is a lot easier in my work to do blast searches against protein databases when I can, since I work with distantly related organisms. But sometimes, I need to work with the actual nucleotide sequences. Until now, I've often done it manually, but it is getting unfeasible due to increasing data amounts.</p>
  <p>The NCBI website is structured such that I suspect automating the process (e.g. input a list of protein accession numbers, get out a fasta file with the corresponding nucleotide seqs) should be possible. At the same time, I am not quite sure what tools to use (biopython?), and if a tool already exists, I would rather not reinvent the wheel. Would you happen to know of anything like this? Thank you for your time!</p>
  <p><strong>Update and WIP:</strong></p>
  <p>I have tried Entrez Utils. In part it works, but I am having some issues scaling it up. For example,</p>
  <p><code>elink -db protein -target nuccore -id WP_244350345.1 | efetch -format fasta_cds_na | awk '/^&gt;/ {printf(&quot;\n%s\n&quot;,$0);next; } { printf(&quot;%s&quot;,$0);}  END {printf(&quot;\n&quot;);}' | grep WP_244350345.1 -A 1</code><br />
  works. (elinks to connect the protein accession WP_244350345.1 with the nucleotide database | efetch to grab the nucleotide sequence, which is a bit of a pain because it grabs the entire genome rather than the cds I want | so I need a bit of cleaning, including using awk to turn the resulting multi-line fasta to single-line fasta | and grep to actually get the cds sequence associated with WP_244350345.1. As I said, it works, and it results in:</p>
  <pre><code>&gt;lcl|NZ_JAFIRA010000020.1_cds_WP_244350345.1_35 [gene=bchM] [locus_tag=JX360_RS09120] [protein=magnesium protoporphyrin IX methyltransferase] [protein_id=WP_244350345.1] [location=35695..36360] [gbkey=CDS]
  ATGCAAGAGAAAAGCATCGTCCGCGACTATTTCAACTCCATTGGCTTCGACCGTTGGCGGCGCATCTACGGAGAGGATGAGGTCAATTTTGTCCAAAAAGACATTCGCTCCGGCCATGCCCGTACCGTAGCCACTGTGTTGGAATGGCTGGGGGATCCCTCAGGTTTGCAGATTTGTGATGCGGGGTGTGGGGTCGGTAGTCTTAGTTTGCCGTTGGCTGCCAAAGGGGCCCGCGTCTTCGCCAGCGACATTTCCGAACAAATGGTGAACGAAGCCCGCCGTCGCCAGCAAATCCAGTTGGGATCCACCGAGAACCCCCAATTTCGGGTCTCGGACTTAGAGGAACTAAGCGGAGAATATGATGCCGTCATTTGCCTGGATGTGATGATCCACTACCCAGAAGCCGATGCCCTGCGCATGCTGGAACATCTAACCCGCCTGGCCCGCTCTCGCTTGATCTTCAGTTTTGCCCCCAAAACCCCTCTGTTGACGTTGCTAAAAAAAGTGGGAGAGTTCTTCCCCGGCCCCAGCAAAGCCACCCGTGCCTACCAACACCGAGAAGCGATTTTGGTGAGCAAATTAACCGAATTGGGCTGGAACGTTCAGCAACGACAGACCATTCGTAGCCGCTTTTACTTTGCCTGGATTCTGGATTTGGCCCGGTGA
  </code></pre>
  <p><em>However</em>, trying to read multiple accession numbers from a file, in a loop, is not working. My current code (essentially the same as before, just replacing WP_244350345.1 with ${LINE}:</p>
  <pre><code>while IFS='' read -r LINE || [ -n &quot;${LINE}&quot; ]; do
  echo &quot;${LINE}&quot;
  elink -db protein -target nuccore -id &quot;${LINE}&quot; | efetch -format fasta_cds_na | awk '/^&gt;/ {printf(&quot;\n%s\n&quot;,$0);next; } { printf(&quot;%s&quot;,$0);}  END {printf(&quot;\n&quot;);}' | grep &quot;${LINE}&quot; -A 1
  done &lt; &quot;results.txt&quot; | more
  </code></pre>
  <p>The <code>echo &quot;${LINE}&quot;</code> part works fine, the rest does not. It has also worked perfectly fine with other parts of Entrez Utils (see below), so I am confused and I would appreciate some help.</p>
  <pre><code>while IFS='' read -r LINE || [ -n &quot;${LINE}&quot; ]; do
  organism=$(elink -db protein -target taxonomy -id &quot;${LINE}&quot; | efetch | tr --delete '\n' | sed 's/    species, /_/g'| sed 's/^.*\. //g' | sed 's/ /_/g')
  echo &quot;${LINE}&quot;
  echo &quot;$organism&quot;
  done &lt; results.txt | more
  </code></pre>
  <p>results in:</p>
  <pre><code>WP_244350345.1
  Thermostichus_vulcanus_cyanobacteria
  WP_073547903.1
  Chroogloeocystis_siderophila_cyanobacteria
  WP_088429095.1
  Halomicronema_hongdechloris_cyanobacteria
  ...
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-05-12T21:21:35.717Z
Comments.Text: Can you please add any solutions you find as an answer instead of as additions to questions? This keeps the Stack Exchange automated algorithms happy, so they don't try to get more answers for an already-answered question.
Comments.UserDisplayName: gringer
Tags: <proteins><gene><ncbi>
Title: How to get the nucleotide sequence for protein NCBI accession code, programatically?
ThreadText: # How to get the nucleotide sequence for protein NCBI accession code, programatically?
  
  (Asked by: username_0 on 2023-05-11)
  
  It is a lot easier in my work to do blast searches against protein databases when I can, since I work with distantly related organisms. But sometimes, I need to work with the actual nucleotide sequences. Until now, I've often done it manually, but it is getting unfeasible due to increasing data amounts.
  
  The NCBI website is structured such that I suspect automating the process (e.g. input a list of protein accession numbers, get out a fasta file with the corresponding nucleotide seqs) should be possible. At the same time, I am not quite sure what tools to use (biopython?), and if a tool already exists, I would rather not reinvent the wheel. Would you happen to know of anything like this? Thank you for your time!
  
  **Update and WIP:**
  
  I have tried Entrez Utils. In part it works, but I am having some issues scaling it up. For example,
  
  `elink -db protein -target nuccore -id WP_244350345.1 | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' | grep WP_244350345.1 -A 1`  
  works. (elinks to connect the protein accession WP\_244350345.1 with the nucleotide database | efetch to grab the nucleotide sequence, which is a bit of a pain because it grabs the entire genome rather than the cds I want | so I need a bit of cleaning, including using awk to turn the resulting multi-line fasta to single-line fasta | and grep to actually get the cds sequence associated with WP\_244350345.1. As I said, it works, and it results in:
  
  ```
  >lcl|NZ_JAFIRA010000020.1_cds_WP_244350345.1_35 [gene=bchM] [locus_tag=JX360_RS09120] [protein=magnesium protoporphyrin IX methyltransferase] [protein_id=WP_244350345.1] [location=35695..36360] [gbkey=CDS]
  ATGCAAGAGAAAAGCATCGTCCGCGACTATTTCAACTCCATTGGCTTCGACCGTTGGCGGCGCATCTACGGAGAGGATGAGGTCAATTTTGTCCAAAAAGACATTCGCTCCGGCCATGCCCGTACCGTAGCCACTGTGTTGGAATGGCTGGGGGATCCCTCAGGTTTGCAGATTTGTGATGCGGGGTGTGGGGTCGGTAGTCTTAGTTTGCCGTTGGCTGCCAAAGGGGCCCGCGTCTTCGCCAGCGACATTTCCGAACAAATGGTGAACGAAGCCCGCCGTCGCCAGCAAATCCAGTTGGGATCCACCGAGAACCCCCAATTTCGGGTCTCGGACTTAGAGGAACTAAGCGGAGAATATGATGCCGTCATTTGCCTGGATGTGATGATCCACTACCCAGAAGCCGATGCCCTGCGCATGCTGGAACATCTAACCCGCCTGGCCCGCTCTCGCTTGATCTTCAGTTTTGCCCCCAAAACCCCTCTGTTGACGTTGCTAAAAAAAGTGGGAGAGTTCTTCCCCGGCCCCAGCAAAGCCACCCGTGCCTACCAACACCGAGAAGCGATTTTGGTGAGCAAATTAACCGAATTGGGCTGGAACGTTCAGCAACGACAGACCATTCGTAGCCGCTTTTACTTTGCCTGGATTCTGGATTTGGCCCGGTGA
  
  ```
  
  *However*, trying to read multiple accession numbers from a file, in a loop, is not working. My current code (essentially the same as before, just replacing WP\_244350345.1 with ${LINE}:
  
  ```
  while IFS='' read -r LINE || [ -n "${LINE}" ]; do
  echo "${LINE}"
  elink -db protein -target nuccore -id "${LINE}" | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);}  END {printf("\n");}' | grep "${LINE}" -A 1
  done < "results.txt" | more
  
  ```
  
  The `echo "${LINE}"` part works fine, the rest does not. It has also worked perfectly fine with other parts of Entrez Utils (see below), so I am confused and I would appreciate some help.
  
  ```
  while IFS='' read -r LINE || [ -n "${LINE}" ]; do
  organism=$(elink -db protein -target taxonomy -id "${LINE}" | efetch | tr --delete '\n' | sed 's/    species, /_/g'| sed 's/^.*\. //g' | sed 's/ /_/g')
  echo "${LINE}"
  echo "$organism"
  done < results.txt | more
  
  ```
  
  results in:
  
  ```
  WP_244350345.1
  Thermostichus_vulcanus_cyanobacteria
  WP_073547903.1
  Chroogloeocystis_siderophila_cyanobacteria
  WP_088429095.1
  Halomicronema_hongdechloris_cyanobacteria
  ...
  
  ```
  
  **Tags:** proteins, gene, ncbi
  
  > **Comments:**
  > **username_1**: Can you please add any solutions you find as an answer instead of as additions to questions? This keeps the Stack Exchange automated algorithms happy, so they don't try to get more answers for an already-answered question.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2023-05-11. Score: 1 | Accepted answer
  
  Yes its easy using Biopython.
  
  What I'm requesting is the sequence must be >9500 otherwise it's discarded. In your case this value will be MUCH bigger, because these are bacterial genomes.
  
  ```
  from Bio import Entrez
  from Bio import SeqIO, AlignIO
  
  def getgenbank (smallGenDict):
      bacteriadb = []
      for bacteria in smallGenDict.keys():
          Entrez.email = "name@email.com"
          with Entrez.efetch(
                  db="nucleotide", rettype="gb", retmode="text", id=bacteria
                  ) as handle:
              seq_record = SeqIO.read(handle, "gb")
              # record = get_CDS(seq_record, bacteria, smallGenDict.keys)
              record = seq_record # better to remove the above function because its outside the question
              if len(record.seq) > 9500:
                  bacteriadb.append(record)
      return bacteriadb
  
  def get_CDS(record, bacteria, smallGenDict):
      for feature in record.features:
          # custom manipulation of input sequences 
          # better to forget this bit of code
          # this whole part of the code has been skipped
          # its capturing the coding sequence (not relevant)
  return record
  
  def main(smallGenDict):
      db = getgenbank(smallGenDict)
      SeqIO.write(db, '/Users/username/Desktop/myoutput.fasta', "fasta")
  
  if __name__ == '__main__':
      smallGenDict = {'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'} # better to automate this bit (different question)
      main(smallGenDict)
  
  ```
  
  Thats it.
  
  **Notes** In the code below, the file smallGenDict is a dictionary of Genbank accession No.s and their sequence.
  `{'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'}`
  This better written here as a list
  `['ABxxxxxx','ABxxxxx2' ...]` say call it `smallGenList`.
  
  Where it says `smallGenDict.keys():` just replace it with `smallGenList` and remove the method `.keys()`, otherwise the code will through a bug.
  The function `get_CDS` is processing the sequence file to remove the non-coding regions. I've not included it because its specific stuff.
  
  
  ---
  
  ### Answer by username_3 on 2023-05-11. Score: 0
  
  You can use R [Rentrez](https://cran.r-project.org/web/packages/rentrez/vignettes/rentrez_tutorial.html), in the link you can find a tutorial for the package usage.
  
  > **Comments:**
  > **username_0**: Thank you! I've tried Entrez (R is not something I am familiar with), but I got stuck.  
  
  ---
  
  