Get a FASTA file from an ID gene or Locus Tag list
2
0
Entering edit mode
4.5 years ago

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 • 2.5k views
ADD COMMENT •
2
Entering edit mode

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
ADD REPLY •
0
Entering edit mode

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!

ADD REPLY •
0
Entering edit mode

See the new answer below.

ADD REPLY •
0
Entering edit mode

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.

ADD REPLY •
0
Entering edit mode

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.

ADD REPLY •
2
Entering edit mode
4.5 years ago
GenoMax 147k

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
ADD COMMENT •
0
Entering edit mode
4.5 years ago

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!
ADD COMMENT •
0
Entering edit mode

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

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 2577 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6