python - Selectively choosing nucleotide sequences from a fasta file?
Question: Using biopython how can I snip genes of my interest from a fasta file if the
gene names are stored in a text file?

    
    
    #extract genes                
    f1 = open('ortholog1.txt','r')
    f2 = open('all.fasta','r')
    f3 = open('ortholog1.fasta','w')
    
    genes = [line.rstrip('\n') for line in f1.readlines()]
    
    
    i=0
    for seq_record in SeqIO.parse(f2, "fasta"):
        if genes[i] == seq_record.id:              
            print genes[i]
            f3.write('>'+genes[i])
            i=i+1
            if i==18:
                break
            f3.write('\n')
            f3.write(str(seq_record.seq))
            f3.write('\n')
    
    
    f2.close()
    f3.close()
    

I was trying the above code. But it has some mistakes and is not generic,
since like `ortholog1.txt` (which contain gene names) there are 5 more similar
files. Also the number of genes in each file varies (not 18 always as here).
Here `all.fasta` is the file which contains all genes. `ortholog1.fasta` must
contain the snipped nucleotide sequence.


Answer: Basically, you can make Biopython do all the work.

I'm going to guess that the gene names in "ortholog1.txt" are exactly the same
as in the fasta file, and that there is one gene name per line. If not, you'd
need to tweak them as necessary to make them line up.

    
    
    from Bio import SeqIO
    
    with open('ortholog1.txt','r') as f:
        orthologs_txt = f.read()
    orthologs = orthologs_txt.splitlines()
    
    genes_to_keep = []
    for record in SeqIO.parse(open('all.fasta','r'), 'fasta'):
        if record.description in orthologs:
            genes_to_keep.append(record)
    
    with open('ortholog1.fasta','w') as f:
        SeqIO.write(genes_to_keep, f, 'fasta')
    

Edit: Here is one way to keep the output genes in the same order as in the
orthologs file:

    
    
    from Bio import SeqIO
    
    with open('all.fasta','r') as fasta_file:
        record_dict = SeqIO.to_dict(open(SeqIO.parse(fasta_file, 'fasta')
    
    with open('ortholog1.txt','r') as text_file:
        orthologs_txt = text_file.read()
    
    genes_to_keep = []     
    for ortholog in orthologs_txt.splitlines():
        try:
            genes_to_keep.append( record_dict[ortholog] )
        except KeyError:
            pass
    
    with open('ortholog1.fasta','w') as output_file:
        SeqIO.write(genes_to_keep, output_file, 'fasta')
    

