How to get the nucleotide sequence for protein NCBI accession code, programatically?

Question

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
...

Answer 1

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 = "[email protected]"
        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 2

You can use R Rentrez, in the link you can find a tutorial for the package usage.