text: user8692: 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
  
  
  user5077: 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=[])]
  
  
  user8629: 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
  
  
original_text: 8692: <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>
   5077: <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>
   8629: <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>
  