利用biopython-NcbitblastnCommandline提取非同义代换

Question

我试图使用NcbitblastnCommandline对一个核苷酸序列进行蛋白质查询，然后报告结果。程序运行正常。但是，结果是，我的查询序列包含XXXXXXXX而不是输入序列。有人知道怎么解决这个问题吗？你知道吗

我使用的代码是：

output=NcbitblastnCommandline(query=QUERY, subject=ALL_SUBJECTS, evalue=0.001, outfmt=5)()[0]
blast_result_record = NCBIXML.read(StringIO(output))
print(output)

我的hsp\u qseq输出如下（很多XXXXX）： Dtligvaitdgnqqimlfsnegkairfaetdvramgrtakgvrgmrvsfasstlxxxxxxxxxxxxxxxxxxxxxxxxxxxxxpetgevlcasangygkrtpvndftkkrggggviaiktserngelvgavsidetkelllisdggtlvrtraaevamtgrnaqgvrlirlseetlvvsixxxxxxxxxxxxxxxxxxxxxxxx见avsnnedtsee

而我的查询实际上是这样的： dtligvaitdgnqqimlfsnegkairfaetdvramgrtakgvrvsfaststlseed公司 ADVENDDSDDNDDSADSSLVSRIVSLVPETGEVLCASANGYGKRTPVNDFPTKKRG公司 GKGVIAIKTSERNGELVGAVSIDETKELLISDGGTLVRTRAAEVAMTGRNAQGVRLI公司 RLSEEETLVGVSIEVEDEELLEGEVDTTETDSEEAVSNNEDTSEE公司

我搞砸了什么吗？你知道吗

Answer 1

我在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)