# 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)

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


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

• 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!
– lzh
Apr 23 '20 at 14:05
• 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
Apr 23 '20 at 14:10
• 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
– lzh
Apr 23 '20 at 14:22
• 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
Apr 23 '20 at 14:29
• 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
– lzh
Apr 23 '20 at 14:33

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=[])]

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