This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get a FASTA file from an ID gene or Locus Tag list

Hello everyone. I have a list of genes after a Differencial expresion analysis and i would like to get a FASTA file to work on OmicsBox (ex Blast2Go). Those genes are ID or Locus tag. By the other hand i have a gff file and a fa genome file.

For example i have the EGR_04594 gen, manually i can go to https://www.ncbi.nlm.nih.gov/gene/?term=EGR_04594 and click on the "Go to nucleotide: FASTA" and go to the site https://www.ncbi.nlm.nih.gov/nuccore/NW_020170409.1?report=fasta&from=687977&to=691293&strand=true to get the sequence, save on a text file and save as *.fa and upload it to OmicsBox, then it works.

The problem is that i have so much genes to do it manually, so the question is:

Is there any way to do this kind of "automatic" with more genes and get only one .fa file?

Thank you all!

rna-seq
• 3,631 views • link
written by ovariohisterectomia • ovariohisterectomia

You can use Entrezdirect. An epost solution can be applied to many entries in one file to get them all.

To get nucleotide sequence

$ esearch -db gene -query "EGR_04594" | elink -target nuccore | efetch -format fasta_cds_na

Problem is that locus_tag returns 163 entries for this genome so I not sure which of these sequences are you interested in or all of them.

$ esearch -db gene -query "EGR_04594 AND Echinococcus granulosus [ORGN]" | elink -target nuccore | efetch -format fasta_cds_na | grep ">" | wc -l
     163

OR

$ esearch -db gene -query "EGR_04594" | elink -target protein | efetch -format fasta
>XP_024351771.1 Glutamate dehydrogenase [Echinococcus granulosus]
MVMCVRKMSKIAQMAPNETIEPSFHEMVQMFAKEGVPHVQKKLLAELPFKGSLTEKEHHIRGIMMSMEQC
DDVLEFNFPIKLDNGEFEMIQGWRAQHSHHVKPCKGGIRFASDVNMNEVKALASLMTYKCSVVDVPYGGA
KAGVKVDPTRYSLGELERICRRFALELAKKNFIGPTTDVPAPDVGTGEKEMAWIADTYANTVGLGDINAL
SCVTGKPVQYGGVEGRSAATGLGLFFGVSNFLNSEKNCKACGLDKPGIAGKTIIIQGFGNVGSYAFRFLQ
EGGAKIIGLIEIDTALYDATGLDFKELIEYRNAHGGVKGFPGAKEVDRKELMYHECDVLICAALQKQITA
ENVDKIKAKVIAEGANGPTTFYAHKQLLKRNIMVIPDLYLNSGGVTVSYFEWLKNLNHVSFGRLNFKYEK
ENAMCLLRNIELDLGKPVKLSKELSDRMKNASELDIVNSGLEYTMERAANQIMEVAERYNLGLDFRTAAY
ILAVEKIFYALHLSGNIF
>EUB60575.1 Glutamate dehydrogenase [Echinococcus granulosus]
MVMCVRKMSKIAQMAPNETIEPSFHEMVQMFAKEGVPHVQKKLLAELPFKGSLTEKEHHIRGIMMSMEQC
DDVLEFNFPIKLDNGEFEMIQGWRAQHSHHVKPCKGGIRFASDVNMNEVKALASLMTYKCSVVDVPYGGA
KAGVKVDPTRYSLGELERICRRFALELAKKNFIGPTTDVPAPDVGTGEKEMAWIADTYANTVGLGDINAL
SCVTGKPVQYGGVEGRSAATGLGLFFGVSNFLNSEKNCKACGLDKPGIAGKTIIIQGFGNVGSYAFRFLQ
EGGAKIIGLIEIDTALYDATGLDFKELIEYRNAHGGVKGFPGAKEVDRKELMYHECDVLICAALQKQITA
ENVDKIKAKVIAEGANGPTTFYAHKQLLKRNIMVIPDLYLNSGGVTVSYFEWLKNLNHVSFGRLNFKYEK
ENAMCLLRNIELDLGKPVKLSKELSDRMKNASELDIVNSGLEYTMERAANQIMEVAERYNLGLDFRTAAY
ILAVEKIFYALHLSGNIF
• 0 views • link
written by GenoMax GenoMax

Thank you for your answer, that is nearly about what im searching for. As you say, in this example "EGR_04594" the command returns 163 entries, but if you look at every entry, there is only one called "EGR_04594". Here are some examples (only first lines are showed):

The code:

$ esearch -db gene -query "EGR_04594" | elink -target nuccore | efetch -format fasta_cds_na

Some examples of the entries:

This one returns the locus tag for EGR_04541

lcl|NW_020170409.1_cds_XP_024351718.1_2 [locus_tag=EGR_04541] [db_xref=GeneID:36340256] [protein=hypothetical protein] [protein_id=XP_024351718.1] [location=complement(join(7473..7625,8373..8450,8927..8970,9039..9154,9309..9514,12530..12816,13863..13992,15596..15616))] [gbkey=CDS]ATGTTGCGTGGGAGGGGGGGGGGTAAACCAGGCGGCTTCCTTGACAAGGGGGCCGATTTGGACGACCTTT TGGATGGATTTGATTTTAATGAAAAGGGCAGTTCTAAAACAGGGACCGCGGCAGCTAGCTCTACGTTGGA GACAAGCTTGGATTCTTTAAAGTTGTCTGTTGAAGCGAAATCTCTGGGCGCCTCGTCGTTGAACGAGTCG TCCAACACTTCGTCCACCACTACTTCCGCCATGCGTGGTGTGCAGAAGAGGAGTGGTGTCGCGGCCTCAT

This returns the locus tag for EGR_04543

>lcl|NW_020170409.1_cds_XP_024351720.1_4 [locus_tag=EGR_04543] [db_xref=GeneID:36340258] [protein=Zinc finger protein 316] [protein_id=XP_024351720.1] [location=complement(join(23776..24284,24356..24439,24510..24675,24871..25305,25377..25456,25523..25595,28083..28151))] [gbkey=CDS]
ATGGAAAGTCCTGGATCAGACAAATCGTCGATTCTAACATCGGACCCTCATATGCCTAATCCTAGAGAGT
CCTTCCCCGAGTCCACACGTCTCCCCCTTGATCAGATGCCCAAAACGATAGAGGATAGCAGTGAAGTTCA
AGACAGTGGCATTGGAGTCTCTCAGCTCTCTTCATCCCTATCGGGACTTCCTTCCTCTCGAAGGAGTCCT
GAATCGACAGGGGCCCACTTACCGAATAACTTTTCCACCCCACTGAATACTCGTGATGCGGTAGTGCCAT
TGGTAGTGCCGCAAGAGAGTATTTCGAACGTTGCGGAGGCTACGAGTTTTGTCACCTCCACCTCAGCGTC

And this entry is for EGR_04494 (this is the entry that im looking for)

>lcl|APAU02000029.1_cds_EUB60575.1_55 [locus_tag=EGR_04594] [protein=Glutamate dehydrogenase] [protein_id=EUB60575.1] [location=complement(join(687977..688221,688300..688645,688736..688975,689052..689710,691257..691293))] [gbkey=CDS]
ATGGTTATGTGCGTGCGTAAAATGTCAAAAATAGCCCAAATGGCACCTAATGAGACTATTGAGCCTTCGT
TCCATGAAATGGTTCAAATGTTCGCCAAAGAGGGCGTTCCTCACGTCCAAAAAAAACTACTGGCTGAGTT
GCCTTTCAAAGGCTCACTAACCGAGAAAGAACACCACATCAGGGGCATTATGATGTCTATGGAACAATGT
GATGATGTTTTGGAGTTCAATTTCCCCATTAAACTGGACAACGGTGAATTTGAGATGATTCAAGGATGGC
GTGCACAACACTCTCATCACGTTAAACCATGCAAGGGGGGCATTCGCTTTGCTTCGGACGTCAATATGAA
CGAAGTCAAGGCTCTGGCTAGTCTGATGACCTATAAGTGCAGTGTCGTTGATGTGCCATATGGTGGCGCT

I dont know why it returns many entries, if i only asked for "EGR_04594", i really appreciate your help, i was watching for something like this many time ago, hope you may still help me with this last question, thank you!

• 0 views • link
written by ovariohisterectomia • ovariohisterectomia

See the new answer below.

• 0 views • link
written by GenoMax GenoMax

Thank you for your kind help! Now that i found that, im trying to perform the same command, but with many genes from a list. I Have a file called "file.txt" with my genes, and its like this (only few genes are shown, for space reasons):

EGR_01804
EGR_05374
EGR_07251
EGR_10508
EGR_07266
EGR_05599
EGR_01271
EGR_11080
EGR_10767
EGR_03982

Then i tried a command, but did not work...i will show:

INPUT

epost -db gene -input file.txt | elink -target protein -name gene_protein_refseq | efetch -format fasta_cds_na

OUTPUT

Non-numeric value found in post input:EGR_01804,EGR_05374,EGR_07251,EGR_10508,EGR_07266,EGR_05599,EGR_01271,EGR_11080,EGR_10767,EGR_03982,EGR_10757,EGR_10836,EGR_09530,EGR_08771,EGR_04836,EGR_03663,EGR_02967,EGR_08116,EGR_10423,EGR_10105,EGR_04953,EGR_00981,EGR_10846,EGR_00908,EGR_08444,EGR_00872,EGR_01543,EGR_09863,EGR_04115,EGR_00982,EGR_06499,EGR_10572,EGR_10875,EGR_06890,EGR_05116,EGR_00641,EGR_04982,EGR_07298,EGR_07536,EGR_05347,EGR_00253,EGR_08711,EGR_10082,EGR_05414,EGR_04568,EGR_10759,EGR_05146,EGR_06230,EGR_08680,EGR_10237,EGR_03250,EGR_10428,EGR_06812,EGR_09445,EGR_06594,EGR_07181,EGR_06729,EGR_00930,EGR_00326,EGR_00236,EGR_05086,EGR_08356,EGR_00203,EGR_11152,EGR_09718,EGR_07253,EGR_08977(...)EGR_02453,EGR_00040,EGR_04995,EGR_04353,EGR_08156,EGR_00825,EGR_06875,EGR_10986,EGR_02017,EGR_06716,EGR_00740,EGR_02558,EGR_01328,EGR_02144,EGR_01878,EGR_07268,EGR_10137,EGR_00892,EGR_00929,EGR_02623,EGR_09132,EGR_01860,EGR_08946,EGR_07998,EGR_06877,EGR_03641,EGR_08264,EGR_01616,EGR_06240

Db value not found in link input
Db value not found in fetch input

Hope you may help me again, im really newbie in this things. Thank you again.

• 0 views • link
written by ovariohisterectomia • ovariohisterectomia

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized. SUBMIT ANSWER is for new answers to original question.

• 0 views • link
written by GenoMax GenoMax

2 answers

I refined previous answer (moved to a comment to provide context). There are two entries that get pulled up. One is from genome and other is from RefSeq. See the two queries below. [Sequence truncated for space reasons]

Using Entrezdirect:

$  esearch -db gene -query "EGR_04594 [LOCUS_TAG]" | elink -target protein -name gene_protein_refseq | efetch -format fasta_cds_na
>lcl|XM_024493843.1_cds_XP_024351771.1_1 [locus_tag=EGR_04594] [db_xref=GeneID:36340309] [protein=Glutamate dehydrogenase] [protein_id=XP_024351771.1] [location=1..1527] [gbkey=CDS]
ATGGTTATGTGCGTGCGTAAAATGTCAAAAATAGCCCAAATGGCACCTAATGAGACTATTGAGCCTTCGT
TCCATGAAATGGTTCAAATGTTCGCCAAAGAGGGCGTTCCTCACGTCCAAAAAAAACTACTGGCTGAGTT

OR

$ esearch -db gene -query "EGR_04594 [LOCUS_TAG]" | elink -target protein  | efetch -format fasta_cds_na
>lcl|XM_024493843.1_cds_XP_024351771.1_1 [locus_tag=EGR_04594] [db_xref=GeneID:36340309] [protein=Glutamate dehydrogenase] [protein_id=XP_024351771.1] [location=1..1527] [gbkey=CDS]
ATGGTTATGTGCGTGCGTAAAATGTCAAAAATAGCCCAAATGGCACCTAATGAGACTATTGAGCCTTCGT
TCCATGAAATGGTTCAAATGTTCGCCAAAGAGGGCGTTCCTCACGTCCAAAAAAAACTACTGGCTGAGTT
>lcl|APAU02000029.1_cds_EUB60575.1_1 [locus_tag=EGR_04594] [protein=Glutamate dehydrogenase] [protein_id=EUB60575.1] [location=complement(join(687977..688221,688300..688645,688736..688975,689052..689710,691257..691293))] [gbkey=CDS]
ATGGTTATGTGCGTGCGTAAAATGTCAAAAATAGCCCAAATGGCACCTAATGAGACTATTGAGCCTTCGT
TCCATGAAATGGTTCAAATGTTCGCCAAAGAGGGCGTTCCTCACGTCCAAAAAAAACTACTGGCTGAGTT
GCCTTTCAAAGGCTCACTAACCGAGAAAGAACACCACATCAGGGGCATTATGATGTCTATGGAACAATGT

Edited: Based on new request I am amending this answer.

Since your identifiers are not accession numbers you can't use the epost method. Try

$ for i in $(cat your_identifier_file); do esearch -db gene -query "${i} [LOCUS_TAG]" | elink -target protein -name gene_protein_refseq | efetch -format fasta_cds_na >> sequence.fa; done
• 0 views • link
written by GenoMax GenoMax

Updating my previous post, i made my gene list with GeneID (no more with Locus Tag), then i have a list like this

36344776
36346103
36336986
36341279
36339830
36346795
36343976
36345820
36338828
36337519
36344362
36344959
36346561
36342521
36346551
36340258
36341314
36342966
36337111
36336623

Then i performed this command

epost -db gene -input "file2.txt" | elink -target protein -name gene_protein_refseq | efetch -format fasta_cds_na

and i finally get what in need, here is (part of) the output:

>lcl|XM_024500329.1_cds_XP_024345260.1_1 [locus_tag=EGR_11080] [db_xref=GeneID:36346795] [protein=Non-capsid protein NS-1] [protein_id=XP_024345260.1] [location=1..588] [gbkey=CDS]
ATGTTTAACCTACGCGCGGTTCCGACGGCCGGCCGTCGCGCAGGCAATCCGGATTTGCGAACCAATCAGC
GTAGAGCACACGCCCGCACCACGGAGACGCACGCATCCCACCCCCCCGGAAGCGGAAGTGAGTCCATTCA
ATGGCTGGAAGACATGTTTTTAGCCAATGAAATCGCCCTTGTCGATTTTGCAATTACGCTTCGCATTATA
ATGAATTGCGAAGATGAGAAAATCAACACTCTTGTGTTGTACGGCCCGACCAATACGGGCAAATCGCTTA
TTTGTAAGCTGACAACGACCTTCCTTGAGCATGGCAGTGTCATGCGCAGGCAGGGGGCATCAGCCTTCGC
TTACGAGAACCTTCTTAATAGGAAGGTTGCGTTAATGGAGGAGCCTGGGATCTGCGCTGCTAACCAGCAG
GATCTGAAGCAGATCCTAGGAGGCGAGACATTTAAGGGCCCCAAAGACATGCAGACGACCCAACAGGCTC
CACAGCCTGTTCAGAGCACTGCCCAACACCCGCCTGCACCTGCGACTCCTACACCACATACATGGCTTAA
AGGTGACACCACCACCTTGTCGGCTTGA
>lcl|XM_024500095.1_cds_XP_024345491.1_1 [locus_tag=EGR_10846] [db_xref=GeneID:36346561] [protein=hypothetical protein] [protein_id=XP_024345491.1] [location=1..693] [gbkey=CDS]
ATGTCGCGACACTTCATGCACTCCGCCAAGGGTACGGATTCCCTACAGCAATCACTGGAAGTGAGGTCGG
TTCGAAGAAGACTAGATGGTCGCAATGTTGCCGTGGACTTGCGTCGTCGTCGGTGCTCCCCTCGGTCAAC
CTCCTGTATCGGCCATCACAACGTAGTCGGTGTTGTAAGAAAGGCACCCAAGAACCAACCGTGCCAATTA
TGGAAATGGCTTGTAACTGAGAAGATATTTGACCACAACGTCAGTGACGAGTTGCTGCAGATCAATGGTT
TCATGACGGCGGGAATGTCGTATCGTCGTGCGGTGGAGATAATCCGGGCGGGTGGCAACTTGGTCCGCCC
CGCTGCGCTGTCGCCGCTGCCTCTACCTGCGCTCTCCTTCGCACGCTTCCATCCACTGCCGCCTATACCG
TCTTCAGCCCCCGTGGTGCCGCAGCCACCACCGGCGCCTCCACTACCACCGCCATATCTCCGCTCACTGC
AGTTCTTTGCCCCTTCATCCGTGCACAAATCATCCATTGGATCCCATCCCCTACTCTCCTCGGCCTCCGC
CGCAGCCGTTGCCTCCACCGCACTACCAGAATGGGGGAATGAGGATGCCTCTGGGTTTCCCCTTCTCTCC
ACCACATTTCCGCCACCCGTTTCTCTCCTCCTCTCTCCACCTCCCACCATCACTGAAATCTAA

So i think this is the way to do this, please correct me if im wrong.

My last question is, is there any way to print the output to a file? Thank you again and hope for the last time!

• 0 views • link
written by ovariohisterectomia • ovariohisterectomia

I amended my answer above with a new solution that will address your questions.

• 0 views • link
written by GenoMax GenoMax

Log in to answer this question.

More posts like this