1
$\begingroup$

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
$\endgroup$
2
$\begingroup$

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
$\endgroup$
11
  • $\begingroup$ 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! $\endgroup$
    – lzh
    Apr 23 '20 at 14:05
  • $\begingroup$ 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 $\endgroup$
    – mrhd
    Apr 23 '20 at 14:10
  • $\begingroup$ 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 $\endgroup$
    – lzh
    Apr 23 '20 at 14:22
  • $\begingroup$ 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? $\endgroup$
    – mrhd
    Apr 23 '20 at 14:29
  • $\begingroup$ 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 $\endgroup$
    – lzh
    Apr 23 '20 at 14:33
1
$\begingroup$

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=[])]
$\endgroup$
4
  • $\begingroup$ This is a very neat comprehension $\endgroup$
    – M__
    Apr 23 '20 at 18:04
  • 1
    $\begingroup$ lol OP chose the mighty awk =p . Thank you for the kind comment $\endgroup$
    – StupidWolf
    Apr 23 '20 at 18:13
  • $\begingroup$ The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records])) $\endgroup$
    – M__
    Apr 23 '20 at 19:14
  • $\begingroup$ Thanks! I try it and it is also work. $\endgroup$
    – lzh
    May 11 '20 at 16:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.