1

我正在尝试使用 NCBIWWW 获得特定蛋白质的 blastP 结果。问题是返回的不是我识别为对齐数据的,我得到的是这个(这是源代码中'Blast_record'的内容);我正在使用从“BioPython 教程和食谱”中获得的代码,我已经在它和互联网上搜索了错误的来源,但我找不到它。我的源代码是这样的;

# biopython
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

# first get the sequence we want to parse from a FASTA file
# f_record = next(SeqIO.parse('m_cold.fasta', 'fasta'))

print('Doing the BLAST and retrieving the results...')
result_handle = NCBIWWW.qblast('blastp', 'tsa', '365176198')

# save the results for later, in case we want to look at it
save_file = open('m_cold_blast.out', 'w')
blast_results = result_handle.read()
save_file.write(blast_results)
save_file.close()

结果文件是:

<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput>
<BlastOutput_program>blastp</BlastOutput_program>
<BlastOutput_version>BLASTP 2.2.31+</BlastOutput_version>
<BlastOutput_reference>Stephen F. Altschul, Thomas L. Madden, Alejandro A. Sch&amp;auml;ffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), &quot;Gapped BLAST and PSI-BLAST: a new generation of protein database search programs&quot;, Nucleic Acids Res. 25:3389-3402.</BlastOutput_reference>
<BlastOutput_db>tsa</BlastOutput_db>
<BlastOutput_query-ID>gi|365176198|gb|AEW67975.1|</BlastOutput_query-ID>
<BlastOutput_query-def>polyprotein [Black queen cell virus]    </BlastOutput_query-def>
<BlastOutput_query-len>171</BlastOutput_query-len>
<BlastOutput_param>
<Parameters>
  <Parameters_matrix>BLOSUM62</Parameters_matrix>
  <Parameters_expect>10</Parameters_expect>
  <Parameters_gap-open>11</Parameters_gap-open>
  <Parameters_gap-extend>1</Parameters_gap-extend>
  <Parameters_filter>F</Parameters_filter>
</Parameters>
</BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
<Iteration_iter-num>1</Iteration_iter-num>
<Iteration_query-ID>gi|365176198|gb|AEW67975.1|</Iteration_query-ID>
<Iteration_query-def>polyprotein [Black queen cell virus]</Iteration_query-def>
<Iteration_query-len>171</Iteration_query-len>
<Iteration_hits>
</Iteration_hits>
<Iteration_stat>
<Statistics>
  <Statistics_db-num>0</Statistics_db-num>
  <Statistics_db-len>0</Statistics_db-len>
  <Statistics_hsp-len>0</Statistics_hsp-len>
  <Statistics_eff-space>0</Statistics_eff-space>
  <Statistics_kappa>-1</Statistics_kappa>
  <Statistics_lambda>-1</Statistics_lambda>
  <Statistics_entropy>-1</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
</BlastOutput_iterations>
</BlastOutput>

现在,如果我使用 BlastN 和上面使用的蛋白质的核苷酸序列进行搜索,我会得到所有匹配的序列、它们的 E 值和分数等。那么为什么在使用 BlastP 时不是这种情况呢?

我对 Python 和 Biopython 都很陌生,对于我的生活,我无法弄清楚我做错了什么。

4

2 回答 2

1

QUERY 365176198是一种蛋白质

数据库是核糖体

什么是转录组 Shotgun 组装 (TSA) 数据库?

TSA 是从原始数据(例如 EST、踪迹和下一代测序技术)计算组装序列的存档。来自完整转录组的重叠序列读数通过计算方法而不是通过克隆 cDNA 的传统克隆和测序组装成转录本。

是否可以通过 BLAST 搜索获得 TSA 序列?

Transcriptome Shotgun Assembly (TSA) BLAST 数据库现已推出。这些序列最初包含在 nt 中,但现在已被分离到一个单独的数据库中。TSA 数据库可从 BLAST 主页的 Basic BLAST 下的核苷酸、tblastn 和 tblastx 链接获得。这些序列在 nt 中不可用。

爆炸风味

爆炸:
     将氨基酸查询序列与蛋白质序列进行比较
     数据库

blastn:将核苷酸查询序列与核苷酸进行比较
     序列数据库

blastx:比较一个核苷酸查询序列
     针对蛋白质序列数据库的阅读框

tblastn:将蛋白质查询序列与核苷酸进行比较
     在所有阅读框中动态翻译的序列数据库

tblastx:比较核苷酸查询的六帧翻译
     针对核苷酸的六帧翻译的序列
     序列数据库。请注意,tblastx 程序不能
     与 BLAST 网页上的 nr 数据库一起使用。

手动应答

手动输入爆破

手动输出爆炸

生物蟒答案

您必须使用“tsa_nt”而不是“tsa”,并且必须使用“tblastn”而不是“blastp”

query = '365176198'
#note: this may take several minutes
result_handle = NCBIWWW.qblast('tblastn', 'tsa_nt', query, format_type="Text")

你得到:

…………

                                                                   分数 E
产生重要比对的序列:(位)值

国标|GAZV01037943.1| Apis mellifera comp13466_c0_seq1 转录... 342 4e-105
国标|GAZF01116856.1| Essigella californica C629542 转录... 179 1e-54
国标|GBYB01008381.1| Fopius arisanus c20283_g1_i1 转录 R... 149 2e-37
国标|GAUO01000423.1| Velia caprai s423_L_1942_0 转录 RNA... 58.9 3e-07
国标|GAXG01028220.1| Gynaikothrips ficorum s28263_L_292921_0 tr... 57.0 9e-07
国标|GAWP01023404.1| Grylloblatta bifratrilecta s23438_L_295244... 52.8 4e-05
国标|GAWZ01143177.1| Gryllotalpa sp。AD-2013 C589197 转录... 45.4 0.002
国标|GAXW01013938.1| Euroleon nostras s13984_L_116369_0 transcr... 45.8 0.006
国标|GAXC01050700.1| 棕榈蓟马 C235436 转录 RNA 序列 42.0 0.017
国标|GAXH01037906.1| Parides eurimedes C235744 转录 RNA ... 40.4 0.069
国标|GBES01007135.1| Dichelops melacanthus Locus_17334_Transcri... 39.7 0.18  
国标|GBXI01014067.1| 瓜实蝇 c16593_g1_i1 transcr... 40.4 0.49  
国标|GAMC01001920.1| Ceratitis capitata comp55379_c0_seq1 mRNA ... 40.4 0.49  
国标|GARL01030594.1| 草地贪夜蛾SEUC25635_TC01转录... 39.7 0.88  
国标|GAZS01034153.1| Acanthoscurria geniculata L2169_T1/2_Turan... 38.5 2.1   
国标|GAZS01034154.1| Acanthoscurria geniculata L2170_T2/2_Turan... 38.5 2.1   
国标|GAYD01030921.1| Blaberus atropos s30958_L_499964_0 transcr... 38.5 2.2   
国标|GAZR01021123.1| Stegodyphus mimosarum L19863_T1/1_Velvet_W... 36.2 3.0   
国标|GBCX01022664.1| Dastarcus helophoroides Unigene14575 trans... 37.7 3.8   
国标|GAYF01148415.1| Nilaparvata lugens C730037 转录 RNA... 37.7 4.5   
国标|GAMK01054259.1| 菜豆 Ref_259_comp8866_c0_seq... 37.0 8.5   
国标|EZ343106.1| 青蒿菌株乌干达 Contig10322.Uhm ... 35.8 8.9   

对齐方式
>gb|GAZV01037943.1| TSA:Apis mellifera comp13466_c0_seq1 转录的 RNA 序列
长度=6998

 分数 = 342 位 (878),期望 = 4e-105,方法:组合矩阵调整。
 同一性 = 164/168 (98%),正数 = 167/168 (99%),差距 = 0/168 (0%)
 帧 = +2

查询 4 ​​YALYRGGVRVKVVTGRGVDFVRATVSPQQTYGSEVAPTTHISTPLAIEQIPIKGVAEFQI 63
             YALYRGGVRVKVVT +GVDFVRATVSPQQTYGS+VAPTTHISTPLAIEQIPIKGVAEFQI
Sbjct 6317 YALYRGGVRVKVVTEKGVDFVRATVSPQQTYGSDVAPTTHISTPLAIEQIPIKGVAEFQI 6496

查询 64 PYYAPCLSSSFRANSETFYYSSGRNNLDIATSPPSINRYYAVGAGDDMDFSIFIGTPPCI 123
             PYYAPCLSSSFRANSETFYYSSGRNNLDI+TSPPSINRYYAVGAGDDMDFSIFIGTPPCI
Sbjct 6497 PYYAPCLSSSFRANSETFYYSSGRNNLDISTSPPSINRYYAVGAGDDMDFSIFIGTPPCI 6676

…………
于 2015-05-26T15:27:59.370 回答
0

我注意到您的输出文件以“.out”结尾。尝试将其保存到 XML 文件中,它会将它们映射到漂亮的列。在输出文件的第一行,您会看到 '?xml. qBLAST 函数默认使用 XML,并且还有一个可选的“文本”参数,尽管格式是一场噩梦。

确实,blastp 和 tsa 是不同的数据库。qblast 模块有一些内置帮助,可以帮助处理不同的参数,可以通过这个来访问。

>>> from Bio.Blast import NCBIWWW
>>> help(NCBIWWW.qblast)

(我会运行您的代码,但 qBLAST 在撰写本文时遇到了问题)

于 2016-11-17T12:44:08.847 回答