to extract dna sequence from a fasta file with gene ids in another location

Question

I have created a little programe to extract selected ids + sequence from a fasta file. The ids of interest are file names that contains several seq for that gene. This is the programe:

import  glob, sys, os
from Bio import SeqIO, SearchIO
from Bio.SeqRecord import SeqRecord
import argparse




def help_function():
    print """
    usage: to_extract_seq_and_id.py [-h] [-i input_file:path to data]  [-r reference file: path_to_file ]
    [-o output_directory: path_to_store_new_file] """
parser = argparse.ArgumentParser()
parser.add_argument('-input_files', '-i',type=str,help='path_to_data')
parser.add_argument('-reference', '-r',type=str,help='path_to_the_fasta_reference_file')
parser.add_argument('-output_directory','-o', type=str, help='path_to_store_new_file')

opts = parser.parse_args()
#first function to check if the files exits.
def check_file_exists(filepath, file_description):
    if not os.path.exists(filepath):
        print("The " + file_description + " (" + filepath + ") does not exist")
        sys.exit(1)
    else:
        print file_description + " detected"
def record_extraction(geneID,reference,output):

    records=list(SeqIO.parse(opts.reference,'fasta'))
    new_reference=output + '/new_reference_common_genes.fa'
    output_handle=open(new_reference, 'a')
    with open (opts.reference, 'rU') as input_handle:
        for record in records:
            recordID=record.id 

            if recordID == geneID:              
                SeqIO.write(record, output_handle, 'fasta')
            else:
                continue
            return new_reference
def main():
    if len(sys.argv) <=2:
        parser.print_help()
        sys.exit()
    else:
        check_file_exists(opts.input_files, 'input_files')
        check_file_exists(opts.reference, 'reference_file')
        check_file_exists(opts.output_directory, 'output_directory')

        files=(glob.glob(opts.input_files + '/*.fa'))

        for f in files:
            database_files=glob.glob(f)[0]
            file_name=os.path.basename(f)
            gene_id=file_name.split('.')
            gene_name=gene_id[0].split('_')
            geneID=gene_name[1] + '_' + gene_name[2]
        print 'The new reference fasta file has been create'




        new_reference=record_extraction(geneID,opts.reference,opts.output_directory)
main()

I have a fasta file with 900 genes, and I want to extract and write in another file 700 of them. The programe runs OK but only write one gene and its sequence into the new file. I don't understand why that loop is running only once. Initially Can anyone help me to correct the problem, please? Thanks in advance.

A [Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve) that allows the reproduction of your problem would be nice. We don't have your files, so how would we know, what your data structure is. And you didn't even clean the code before posting. — Mr. T, May 24 '18 at 16:07

score 1 · Answer 1 · answered May 24 '18 at 16:00

The question is abit unclear but maybe this is what you are looking for:

results = []
for record in records:
    recordID=record.id  #.split('_')

    if recordID == geneID:
        results.append(record)
    else:
        continue           
SeqIO.write(" ".join(text for text in results), output_handle, 'fasta')
return new_reference

If this is not what you are looking for. Please give more details for what are the problem and solution you need.

Ana · Accepted Answer · 2018-06-04T15:13:57.400

The problem I had was a problem of indentation. If you have a look the code above, you can see that the def record_extraction() was not called within the for loop in the main function (def main()). I have changed that indentation, and now it does work perfectly. See above the new script:

import  glob, sys, os
from Bio import SeqIO, SearchIO
from Bio.SeqRecord import SeqRecord
import argparse

def help_function():
    print """
    usage: to_extract_seq_and_id.py [-h] [-i input_file:path to data]  [-r reference file: path_to_file ]
    [-o output_directory: path_to_store_new_file] """
parser = argparse.ArgumentParser()
parser.add_argument('-input_files', '-i',type=str,help='path_to_data')
parser.add_argument('-reference', '-r',type=str,help='path_to_the_fasta_reference_file')
parser.add_argument('-output_directory','-o', type=str, help='path_to_store_new_file')

opts = parser.parse_args()
#first function to check if the files exits.
def check_file_exists(filepath, file_description):
    if not os.path.exists(filepath):
        print("The " + file_description + " (" + filepath + ") does not exist")
        sys.exit(1)
    else:
        print file_description + " detected"
def record_extraction(geneID,reference,output,genelist):

    records=list(SeqIO.parse(opts.reference,'fasta'))
    new_reference=output + '/new_reference_common_genes.fa'
    output_handle=open(new_reference, 'a')
    with open (opts.reference, 'rU') as input_handle:
        for record in records:
            recordid=record.id.split('_')
            recordID=recordid[0]+'_'+recordid[1]                
            if recordID in genelist: 

                SeqIO.write(record, output_handle, 'fasta')
            else:
                continue
        return new_reference    
def main():
    if len(sys.argv) <=2:
        parser.print_help()
        sys.exit()
    else:
        check_file_exists(opts.input_files, 'input_files')
        check_file_exists(opts.reference, 'reference_file')
        check_file_exists(opts.output_directory, 'output_directory')
        #run the programme
        files=(glob.glob(opts.input_files + '/*.fa'))

        for f in files:
            database_files=glob.glob(f)[0]
            file_name=os.path.basename(f)
            gene_id=file_name.split('.')
            gene_name=gene_id[0].split('_')
            geneID=gene_name[1]+'_'+gene_name[2]
            genelist=[]
            if geneID not in genelist:
                genelist.append(geneID)
            new_reference=record_extraction(geneID,opts.reference,opts.output_directory,genelist)

    print 'The new reference fasta file has been create'        

main()

to extract dna sequence from a fasta file with gene ids in another location

2 Answers2