Hi,
I have a bunch of EC numbers based on which I would like to download the corresponding genes from several species.
An example of a search would be "Bacillus[ORGN] AND 1.7.7.2[EC/RN Number]" and if I use:
esearch -db nuccore -query "Bacillus[ORGN] AND 1.7.7.2[EC/RN Number]"
I can see 7 matches (the same that I see using the search on the NCBI website). These are all genomic sequences including genes/CDS with the EC number I am looking for.
My question is: how can I retrieve the sequences of the CDS, not the full genome?
I have tried to use efetch -format gene_fasta but I get the sequences from all CDS and all genomes as a single fasta. It would be okay if the EC number was in the name of the sequences but this info is not displayed (only have the locus_tag, location, sometines gene).
Thanks for the help!
•
written
by
agathejouet •
1 answer
Try:
$ time esearch -db nuccore -query "Bacillus[ORGN] AND 1.7.7.2[EC/RN Number]" | efetch -format gbc | xtract -insd CDS gene product feat_location sub_sequence | sed 's/ /_/g' > file
1m53s
$ head -2 file
CP017247.1 - Chromosomal_replication_initiator_protein_DnaA 34..1341 TTGGGGCAGATCGAAAAAAAATTGAGCAAGCCCAGCTTTGAAACATGGATGAAATCGACAAAGGCCCATTCATTGCAGGGCGATACGCTGATCATCACCGCACCGAACGAGTTCGCCAGAGACTGGCTTGAATCAAGATACCTGCACCTGATCGCCGATACGATCTATGATCTGACAGGAGAAGAATTGAGCATTAAATTTGTCATTCCTCAGAATCAAAATGAAGAAGATTTTATGCCAAAGTCTCCAATCAAAAAAATGTCGAAAGAAGAACCGGCTGATTTTCCGCAAAACATGCTGAATCCCAAATATACATTTGATACGTTCGTTATCGGTTCAGGAAACCGATTCGCCCACGCAGCGTCTTTGGCAGTGGCTGAAGCCCCGGCGAAAGCTTACAATCCGCTGTTTATTTACGGGGGAGTCGGACTTGGAAAGACTCACTTAATGCATGCGATCGGGCACTATGTCATCGATCACAATCCATCTGCAAAAGTGGTTTATTTGTCATCTGAGAAATTTACAAATGAGTTCATTAACTCGATCCGTGACAATAAAGCTGTCGATTTTCGCAATCGCTATCGAAATGTTGACGTTCTTTTAATAGACGATATTCAATTTTTAGCCGGAAAAGAACAGACGCAAGAGGAATTTTTCCATACGTTTAATACGCTGCATGAAGAAACAAAGCAGATTGTCATTTCCAGCGACCGGCCTCCAAAAGAGATCCCAACGCTTGAAGACCGTTTGCGCTCCCGTTTTGAATGGGGATTGATCACTGACATCACGCCTCCTGATCTGGAAACAAGAATTGCGATTTTAAGAAAGAAAGCAAAAGCAGAAGGACTTGATATCCCGAATGAAGTCATGCTTTATATTGCCAATCAGATCGACAGCAATATCAGGGAGCTGGAAGGGGCATTAATCAGGGTTGTCGCATATTCTTCCTTGATCAATAAAGACATTAACGCCGATCTGGCTGCTGAAGCTTTGAAAGATATCATTCCTTCTTCAAAGCCGAAAGTCATTACGATCAAAGACATCCAAAGAATCGTCGGCCAGCAGTTCAATATCAAGCTGGAGGATTTCAAGGCGAAGAAACGGACAAAATCGGTGGCTTTTCCGCGGCAGATCGCTATGTATCTATCAAGAGAAATGACAGATTCTTCTCTTCCGAAGATCGGCGAAGAATTTGGCGGACGCGACCACACGACGGTCATCCATGCCCATGAGAAAATATCAAAACTGCTGAGCGATGATGAACAGCTTCAGCAGCAGATTAAAGAAATTAAAGAGCAGCTGAGATAA
CP017247.1 dnaN DNA-directed_DNA_polymerase 1518..2654 ATGAAGTTTACAATTCAAAAAGACCGCCTGGTCGAAAGTGTCCAAGATGTGTTAAAAGCCGTTTCTTCAAGAACGATGATTCCGATCTTAACCGGTATTAAAATCGTGGCCTCTGATGAAGGGGTCTCTCTGACAGGCAGCGATTCCGATATTTCGATTGAATCGTTTATCCCGAAAGAAGACGGCGATTTAGAGATCGTGACAATTGAACAGCCCGGCAGCATTGTGCTTCAAGCCCGTTTTTTCAGTGAAATTGTCAAAAAGCTGCCGATGTCAACAGTGGAAATCGAGGTTCAAAATCAATACTTAACGATCATCCGCTCCGGCAAAGCAGAGTTTAACTTAAACGGTTTGGATGCAAGCGAATATCCGCTTTTGCCGCAAATTGAAGAGCATCACGCTTTTCAAATTCCGACCGATCTGCTGAAAAACCTGATCCGCCAAACCGTTTTTGCAGTGTCCACCTCAGAAACACGCCCAATCTTGACAGGTGTAAACTGGAATGTCACTGGCGGTGAATTAATATGCACTGCAACGGATAGTCATCGTCTTGCGCTAAGGAAAGCTAAGCTCGACATTAACGAAGACAGTTCATACAATGTCGTCATCCCAGGAAAAAGCTTAACCGAGCTCAGCAAAATCCTTGATGACCATCAGGAGCTTGTCGATATTGTGATTACCGAAACACAAGTGCTGTTTAAAACAAAAAACGTTCTGTTTTTCTCCAGGCTTCTTGACGGAAACTATCCGGATACGAACCGCCTGATTCCTCAGGAAAGCAAAACGAACTTGATTGTCAATACTAAGGAATTTCTCCAGGCGATCGACAGGGCTTCGCTTTTGGCGAGAGAAGGACGAAACAACGTGGTGAAACTTTCCGCAGCCGCCAATGAGTCGATCGAAATTTCTTCAAACTCTCCGGAAATCGGAAAGGTTGTTGAAACGGTGAATGCCGAGCAGATCGAAGGGGAAGACTTAAAGATATCCTTTAGTCCGAAATATATGCTGGATGCCCTTAAAGTTCTTGAAGGAGAGGACATTCATGTAAGCTTCACAGGCGCTATGAGGCCTTTTCTGATCCGTACGCCGAATGACGATTCGATCGTCCAATTAATTCTTCCTGTCCGGACGTATTAA
You will get a table with 5 columns with -- one row per gene.
Gene
The command, indeed, retrieves all genes from Bacillus sequences that have the required EC# -- 1 gene sequence per line. You just need to grep the output, e.g., here I got 19 sequences:
time esearch -db nuccore -query "Bacillus[ORGN] AND 1.7.7.2[EC/RN Number]" | efetch -format gbc |
xtract -insd CDS gene product feat_location sub_sequence | sed 's/ /_/g' | grep -i ferredoxin > file**
real 1m55.001s
user 0m5.944s
sys 0m2.172s
$ wc -l file
19 file
$ head -1 file
CP017247.1 - Ferredoxin 2511566..2511814 ATGCCAAAATATACAATTGTAGACAAAGATACGTGCATCGCATGCGGAGCTTGTGGTGCTGCGGCTCCTGATATTTATGATTACGACGATGAGGGAATCGCATTTGTCACCCTTGACGACAATCAGGGTGTCGTCGAAGTCCCTGACGTCTTAGAAGAAGACATGATGGACGCGTTTGAAGGCTGTCCTACAGATTCGATCAAAGTTGCGGATGAGCCGTTCGAAGGCGACCCGCTTAAACACGAATAA
Log in to answer this question.
-
Finding the paper published from the SRA run riles
written by scholaraspect2008 •
-
Why are number of CDS smaller than corresponding genes [M. tb]
written by Alewa
-
How to download all accession ID of CDS from all viral genome
written by Kumar
-
How to retrieved protein Fasta sequence from accession number by Entrez
written by Nelo •
-
Download discontinued gene record from NCBI via E-utils
written by lmlukoseviciute •
-
How to limit sequences downloaded from NCBI to published data only
written by sbombin •
-
How do I loop over a list of gene symbols to extract sequences from NCBI?
written by MAPK
-
Extracting multiple features using NCBI's e-utilities?
written by ThePresident
-
E-utilities for obtain gene sequences from the gene database
written by dllopezr •
-
Fetching multiple FASTA from NCBI
written by mikael.lenz.strube •
I don't think this strategy is going work since the result from your query is not the specific gene but the genome accession. Trying to think about a way to do this.
The thing is that if I search the gene database, I do not get any results. I am now pulling the features from the genomes using efetch -format ft, retrieving the locus_tag associated with the EC number I am looking for and the using that locus_tag to get the fasta sequence from the sequences obtained using efetch -format gene_fasta. Probably not the most effective way but it appears to be working.