Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:05:55.027Z
Answers.Text: 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!
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:10:26.863Z
Answers.Text: 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`
Answers.UserDisplayName: mrhd
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:22:17.303Z
Answers.Text: 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
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:29:40.637Z
Answers.Text: 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?
Answers.UserDisplayName: mrhd
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:33:23.533Z
Answers.Text: 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
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:41:59.923Z
Answers.Text: 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)?
Answers.UserDisplayName: mrhd
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:43:20.420Z
Answers.Text: 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?
Answers.UserDisplayName: mrhd
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:48:36.957Z
Answers.Text: 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.
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:55:33.020Z
Answers.Text: Thank you, I try it with simpler regular expression, and it got results. Great!!! Thank you very much.
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T14:58:36.073Z
Answers.Text: 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.
Answers.UserDisplayName: mrhd
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T15:14:38.440Z
Answers.Text: Anyway, it saves a lot time to extract sequence, not just find, copy and paste, you make it more interesting. Thank you!
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T07:09:12.713Z
Answers.LastActivityDate: 2020-04-23T07:09:12.713Z
Answers.UserDisplayName: mrhd
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-05-11T16:02:17.857Z
Answers.Text: Thanks! I try it and it is also work.
Answers.UserDisplayName: user8692
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T18:04:36.913Z
Answers.Text: This is a very neat comprehension
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T18:13:35.853Z
Answers.Text: lol OP chose the mighty awk =p . Thank you for the kind comment
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T19:14:18.127Z
Answers.Text: The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records]))
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-23T11:16:44.747Z
Answers.LastActivityDate: 2020-04-23T11:16:44.747Z
Answers.UserDisplayName: StupidWolf
Body: <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>
  
Tags: <python><genome><sequence-analysis><perl>
Title: 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)
ThreadText: # 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.  
  
  ---
  
  