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