извлечь последовательность ДНК из файла fasta с идентификаторами генов в другом месте - PullRequest
Новая
       9

извлечь последовательность ДНК из файла fasta с идентификаторами генов в другом месте

0 голосов
/ 24 мая 2018

Я создал небольшую программу для извлечения выбранных идентификаторов + последовательности из файла fasta.Интересующие идентификаторы - это имена файлов, которые содержат несколько последовательностей для этого гена.Это программа: 

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

У меня есть файл фаста с 900 генами, и я хочу извлечь и записать в другой файл 700 из них.Программа работает нормально, но записывает только один ген и его последовательность в новый файл.Я не понимаю, почему этот цикл работает только один раз.Изначально может кто-нибудь помочь мне исправить проблему, пожалуйста?Заранее спасибо.

Ответы [ 2 ]

0 голосов
/ 04 июня 2018

У меня была проблема с отступом.Если вы посмотрите приведенный выше код, то увидите, что def record_extraction() не был вызван в цикле for в основной функции ( def main () ).Я изменил этот отступ, и теперь он работает отлично.Смотрите выше новый скрипт: 

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()
0 голосов
/ 24 мая 2018

Вопрос немного неясен, но, возможно, это то, что вы ищете: 

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

Если это не то, что вы ищете.Пожалуйста, дайте более подробную информацию о том, какие проблемы и решения вам нужны.

...