Использование биопиона NcbitblastnCommandline для извлечения несинонимичных замен - PullRequest
0 голосов
/ 07 февраля 2019

Я пытаюсь использовать NcbitblastnCommandline, чтобы взорвать запрос белка против нуклеотидной последовательности, а затем сообщить о попадании.Программа запустилась без ошибок.Однако в результате моя последовательность запросов оказалась содержит XXXXXXXX вместо моей входной последовательности.Кто-нибудь знает, как это решить?

1002 * Код я использовал:
output=NcbitblastnCommandline(query=QUERY, subject=ALL_SUBJECTS, evalue=0.001, outfmt=5)()[0]
blast_result_record = NCBIXML.read(StringIO(output))
print(output)

мой выход hsp_qseq выглядит так (много ХХХХХ): DTLIGVAITDGNQQIMLFSNEGKAIRFAETDVRAMGRTAKGVRGMRVSFASSTLXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXPETGEVLCASANGYGKRTPVNDFPTKKRGGKGVIAIKTSERNGELVGAVSIDETKELLLISDGGTLVRTRAAEVAMTGRNAQGVRLIRLSEEETLVGVVSIXXXXXXXXXXXXXXXXXXXXSEEAVSNNEDTSEE

1007 * В то время как мой QUERY на самом деле выглядит следующим образом: DTLIGVAITDGNQQIMLFSNEGKAIRFAETDVRAMGRTAKGVRGMRVSFASSTLSEEDADVENDDSDDNDDSADSSLVSRIVSLVVVPETGEVLCASANGYGKRTPVNDFPTKKRG GKGVIAIKTSERNGELVGAVSIDETKELLLISDGGTLVRTRAAEVAMTGRNAQGVRLI RLSEEETLVGVVSIEVE * 0 * * 8В * ВНУТРИ * * ВНУТРИ * * ВЗЯЛКА * 100В * ВЗЯТКА * * ВНУТРИ * * ВЗЯЛКА?

1 Ответ

0 голосов
/ 13 февраля 2019

Мой способ обойти несинонимную замену после TBLASTN ниже:


from Bio.Blast.Applications import NcbitblastnCommandline
from io import StringIO
from Bio.Blast import NCBIXML
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO
import os
from Bio.SeqRecord import SeqRecord
from Bio.Seq import MutableSeq

file="subject.fasta"
QUERY="query.fasta"
for seq_record in SeqIO.parse(open(QUERY, mode='r'), 'fasta'):
    query_list = seq_record.seq
output=NcbitblastnCommandline(query=QUERY, subject=file, max_target_seqs=1, evalue=0.001,outfmt=5)()[0]
blast_result_record = NCBIXML.parse(StringIO(output))
count = 0
for blast_record in blast_result_record:
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            if count > 0: continue 
            count+=1
            mutation_list = []
            for index in range(len(hsp.match)):
                match_aa = hsp.match[index]
                str_match = MutableSeq(hsp.query[0:index])
                number_gap = str_match.count("-")
                index_in_query = index
                match_query = hsp.query[index]
                if number_gap > 0:
                    index_in_query = len(str_match) - number_gap
                if match_aa == " " or match_aa == "+":
                    if match_query == "-":
                        mutation = str(index_in_query + 1) + query_list[index_in_query] + 'ins' +  hsp.sbjct[index]
                    else:
                        mutation = str(index_in_query + 1) + query_list[index_in_query] + '>' +  hsp.sbjct[index]
                    mutation_list.append(mutation)
            print(mutation_list)
...