<p>If you can use <code>awk</code>, you could use this simple script:</p>

<pre><code>awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^&gt;/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta
</code></pre>

<p>You can also put the code inside a script file:</p>

<p>extract_genes.awk</p>

<pre><code>#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^&gt;/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?
</code></pre>

<p>And call it like this:</p>

<pre><code>awk -f extract_genes.awk gene_list cds.fasta
</code></pre>

CC BY-SA 4.0
2020-04-23T14:05:55.027Z
hi, thank you for your answer. Where to find my output file? After I run it, I don't know where to find the output file when finished. Thank you!
user8692
CC BY-SA 4.0
2020-04-23T14:10:26.863Z
Happy to help. The output is written to stdout. You can forward it into a file like this: `awk -f extract_genes.awk gene_list cds.fasta > output.fasta`
mrhd
CC BY-SA 4.0
2020-04-23T14:22:17.303Z
Thank you. However, I don't know why the output.fasta file is empty. I use command like this: awk -f gene_cds_extract.awk gene_list_red_green_cpsf100.txt all_species_longset_cds.fa > example_cds.fa
user8692
CC BY-SA 4.0
2020-04-23T14:29:40.637Z
That's odd. Are the input files formatted the same as in your question? And are there overlaps between the gene names and the fasta headers?
mrhd
CC BY-SA 4.0
2020-04-23T14:33:23.533Z
Thank you for your quick reply. input cds.fast format like this: >PNS24249 cds supercontig:Brachypodium_distachyon_v3.0:Bd1_centromere_containing_Bradi1g41430:24784:29531:-1 gene:BRADI_41430s00200v3 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGCGGCCGCGCGGCAGAGGAGACGACGGGGACGACGACGAAGACGAAGAGCCCTCGCGT
CAACGGAGCGCGCCGGAAGCCGACCCTGATCCCGAGCCCGCGCCAGCGCAGCAGCCCACA
CGAGCGCCCCTCAGCAGCCTAGTCGTCAGGCCGCCGCCGCCGCAGGACAACGGCGCCTCC     input list.txt: BRADI_41430s00200v3 like this. no gene overlap, I extract longest cds in cds.fa file
user8692
CC BY-SA 4.0
2020-04-23T14:41:59.923Z
I just made a test with your input and the script above and it outputs the sequence. Which operating system do you use, and which awk version (I use Lubuntu 18.04 and GNU Awk 4.1.4)?
mrhd
CC BY-SA 4.0
2020-04-23T14:43:20.420Z
One thing you could try is to use a simpler (less specific) regular expression. If you replace `re = "\\sgene:"g"\\s"` with `re = "gene:"g` does it output anything?
mrhd
CC BY-SA 4.0
2020-04-23T14:48:36.957Z
I use mac terminal. I don't know it already have awk or not. Sorry, It's my first time to use awk in terminal.
user8692
CC BY-SA 4.0
2020-04-23T14:55:33.020Z
Thank you, I try it with simpler regular expression, and it got results. Great!!! Thank you very much.
user8692
CC BY-SA 4.0
2020-04-23T14:58:36.073Z
Aha! Great that it worked. Maybe awk on Mac understands special characters in regular expressions (in this case "\s" for a space) differently than awk on Linux.
mrhd
CC BY-SA 4.0
2020-04-23T15:14:38.440Z
Anyway, it saves a lot time to extract sequence, not just find, copy and paste, you make it more interesting. Thank you!
user8692
CC BY-SA 4.0
2020-04-23T07:09:12.713Z
2020-04-23T07:09:12.713Z
mrhd
<p>If you read it in python, they would take care of that:</p>

<pre><code>from Bio import SeqIO                                                                 
records = list(SeqIO.parse("test.fa","fasta")) 

records

[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm002', name='Zm002', description='Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm001', name='Zm001', description='Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding', dbxrefs=[])]
</code></pre>

<p>Get the gene and sequence.. maybe there's a easier to do this in biopython, this is like my 10th time using it:</p>

<pre><code>def format(rec):
    des = rec.description.split()
    gene = [entry[5:] for entry in des if "gene:" in entry][0]
    return("&gt;"+gene+"\n"+str(rec.seq))

[format(rec) for rec in records]

['&gt;Zm1d035916\nATCGGCAT', '&gt;Zm1d048\nATGCGGCA']
</code></pre>

<p>Write it out:</p>

<pre><code>f=open('example.fasta','w')
f.writelines("\n".join([format(rec) for rec in records]))
f.close()

list(SeqIO.parse("example.fasta","fasta"))                                           
[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm1d035916', name='Zm1d035916', description='Zm1d035916', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm1d048', name='Zm1d048', description='Zm1d048', dbxrefs=[])]
</code></pre>

CC BY-SA 4.0
2020-05-11T16:02:17.857Z
Thanks! I try it and it is also work.
user8692
CC BY-SA 4.0
2020-04-23T18:04:36.913Z
This is a very neat comprehension
M__
CC BY-SA 4.0
2020-04-23T18:13:35.853Z
lol OP chose the mighty awk =p . Thank you for the kind comment
StupidWolf
CC BY-SA 4.0
2020-04-23T19:14:18.127Z
The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records]))
M__
CC BY-SA 4.0
2020-04-23T11:16:44.747Z
2020-04-23T11:16:44.747Z
StupidWolf
<p>i have a question when i want to extract cds sequence using gene id. but cds file is not just start with >gene is, it has many other annotation. the only same  is star with gene:</p>

<p>cds fasta:</p>

<pre><code>&gt;Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding 
ATCGGCAT
&gt;Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA
</code></pre>

<p>gene_list</p>

<pre><code>Zm1d035916
Zm1d048 
</code></pre>

<p>how to get result like</p>

<pre><code>&gt;Zm1d035916
ATCGGCAT
&gt;Zm1d048 
ATGCGGCA
</code></pre>

<python><genome><sequence-analysis><perl>
how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)
# how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)

(Asked by: username_0 on 2020-04-23)

i have a question when i want to extract cds sequence using gene id. but cds file is not just start with >gene is, it has many other annotation. the only same  is star with gene:
cds fasta:
>Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding 
ATCGGCAT
>Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA

gene_list
Zm1d035916
Zm1d048 

how to get result like
>Zm1d035916
ATCGGCAT
>Zm1d048 
ATGCGGCA

**Tags:** python, genome, sequence-analysis, perl

---

## 2 answers

### Answer by username_1 on 2020-04-23. Score: 2 | Accepted answer

If you can use awk, you could use this simple script:
awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

You can also put the code inside a script file:
extract_genes.awk
#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

And call it like this:
awk -f extract_genes.awk gene_list cds.fasta

> **Comments:**
> **username_0**: hi, thank you for your answer. Where to find my output file? After I run it, I don't know where to find the output file when finished. Thank you!  
> **username_1**: Happy to help. The output is written to stdout. You can forward it into a file like this: `awk -f extract_genes.awk gene_list cds.fasta > output.fasta`  
> **username_0**: Thank you. However, I don't know why the output.fasta file is empty. I use command like this: awk -f gene_cds_extract.awk gene_list_red_green_cpsf100.txt all_species_longset_cds.fa > example_cds.fa  
> **username_1**: That's odd. Are the input files formatted the same as in your question? And are there overlaps between the gene names and the fasta headers?  
> **username_0**: Thank you for your quick reply. input cds.fast format like this: >PNS24249 cds supercontig:Brachypodium_distachyon_v3.0:Bd1_centromere_containing_Bradi1g41430:24784:29531:-1 gene:BRADI_41430s00200v3 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein ATGCGGCCGCGCGGCAGAGGAGACGACGGGGACGACGACGAAGACGAAGAGCCCTCGCGT CAACGGAGCGCGCCGGAAGCCGACCCTGATCCCGAGCCCGCGCCAGCGCAGCAGCCCACA CGAGCGCCCCTCAGCAGCCTAGTCGTCAGGCCGCCGCCGCCGCAGGACAACGGCGCCTCC     input list.txt: BRADI_41430s00200v3 like this. no gene overlap, I extract longest cds in cds.fa file  
> **username_1**: I just made a test with your input and the script above and it outputs the sequence. Which operating system do you use, and which awk version (I use Lubuntu 18.04 and GNU Awk 4.1.4)?  
> **username_1**: One thing you could try is to use a simpler (less specific) regular expression. If you replace `re = "\\sgene:"g"\\s"` with `re = "gene:"g` does it output anything?  
> **username_0**: I use mac terminal. I don't know it already have awk or not. Sorry, It's my first time to use awk in terminal.  
> **username_0**: Thank you, I try it with simpler regular expression, and it got results. Great!!! Thank you very much.  
> **username_1**: Aha! Great that it worked. Maybe awk on Mac understands special characters in regular expressions (in this case "\s" for a space) differently than awk on Linux.  
> **username_0**: Anyway, it saves a lot time to extract sequence, not just find, copy and paste, you make it more interesting. Thank you!  

---

### Answer by username_2 on 2020-04-23. Score: 1

If you read it in python, they would take care of that:
from Bio import SeqIO                                                                 
records = list(SeqIO.parse("test.fa","fasta")) 

records

[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm002', name='Zm002', description='Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm001', name='Zm001', description='Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding', dbxrefs=[])]

Get the gene and sequence.. maybe there's a easier to do this in biopython, this is like my 10th time using it:
def format(rec):
    des = rec.description.split()
    gene = [entry[5:] for entry in des if "gene:" in entry][0]
    return(">"+gene+"\n"+str(rec.seq))

[format(rec) for rec in records]

['>Zm1d035916\nATCGGCAT', '>Zm1d048\nATGCGGCA']

Write it out:
f=open('example.fasta','w')
f.writelines("\n".join([format(rec) for rec in records]))
f.close()

list(SeqIO.parse("example.fasta","fasta"))                                           
[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm1d035916', name='Zm1d035916', description='Zm1d035916', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm1d048', name='Zm1d048', description='Zm1d048', dbxrefs=[])]

> **Comments:**
> **username_3**: This is a very neat comprehension  
> **username_2**: lol OP chose the mighty awk =p . Thank you for the kind comment  
> **username_3**: The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records]))  
> **username_0**: Thanks! I try it and it is also work.  

---

