BioHackersNet(<-) |
. |
[../] |
(->) |
Chap 3.1 BLAST
이것봐요. 모두 BLAST 좋아하죠? 내가 얘기하는건, 당신의 sequence 중의 하나와 다른 알려진 sequence 사이의 비교를 어떻게 더 쉽게 더 쉽게 할 수 있는 가에 관한 것입니다. 빌어먹게도, 만약 내가 그 code를 쓴다면 아마도 하루 반 정도가 걸릴 것이지만, 결과는 그다지 만족스럽지 않을 겁니다. 그러나, 물론 이 section 은 BLAST 가 얼마나 멋진가에 관한 것이 아닙니다. 왜냐하면 우리는 이미 잘 알고 있으니까요. BLAST의 문제에 관한 것입니다. 많은 run 에 의해서 생성된 data 들을 처리하고 BLAST run을 자동화시키는 것은 일반적으로 정말 어려울 수도 있습니다.
다행히, BioPython을 만드는 친구들은 이것을 너무나 잘 압니다. 그래서 그들은 BLAST를 더 쉽게 다루는 많은 tool 들을 개발했습니다. 이 section 은 이러한 tool을 어떻게 사용해서 유용한 것들을 할 수 있는지를 설명합니다.
3.1.1 BLAST를 internet을 통해서 돌리기
BLAST를 자동으로 돌리는 첫 번째 step 은 모든 것을 python script 로부터 accessible 하도록 하는 것입니다. 그래서, biopython 은 당신이 WWW version 의 BLAST(http://www.ncbi.nlm.nih.gov/BLAST/) 를 당신의 python script 로부터 direct로 돌릴 수 있도록 하는 code를 가지고 있습니다.
특히 BLAST를 whole BLAST queue과 separate results page를 가지고 script로부터 다루는 것은 아주 힘들 수 있기 때문에, 이것은 아주 멋집니다. BioPythonCodeAnalysis를 NCBI 의 변화에 up to date 하게 유지하는 수고로 충분합니다!
WWW version of BLAST와 blast function을 다루는 code는 Bio.Blast.NCBIWWW module에서 찾을 수 있습니다. 우리가 BLAST의 정보를 FASTA format 으로 원한다고 해 봅시다. 첫째, 우리는 FASTA file 의 정보를 필요로 합니다. 가장 쉬운 방법은 Fasta iterator를 사용해서 record를 얻는 것입니다. 이 예에서 우리는 file로부터 Iterator (section 2.4.3에서 FASTA iterator를 보다 자세히 설명함) 를 사용해서 단지 처음의 FASTA record를 잡아낼 것입니다.
from Bio import Fasta
file_for_blast = open('m_cold.fasta', 'r')
f_iterator = Fasta.Iterator(file_for_blast)
f_record = f_iterator.next()이제 우리는 FASTA record를 얻었습니다. 우리는 이것은 blast 할 준비가 되었습니다. 가장 간단히 가능한 BLAST(모든 non-redundant database에 대해서 FASTA file 의 simple blastn)를 하는 code 는 다음과 같습니다.
from Bio.Blast import NCBIWWW
b_results = NCBIWWW.blast('blastn', 'nr', f_record)
역자주 : 이 작업이 끝나고 나서 윈도우에서는 화면을 최대화했을 때에만 IDLE 이 정상적으로 보이는 문제점이 있었습니다.
Blast function 은 역시 기본적으로 basic BLAST page(http://www.ncbi.nlm.nih.gov/blast/blast.cgi?Jform=0)에서 설정될 수 있는 parameter 들과 상응하는 몇 개의 다른 option argument를 가질 수 있습니다.
그러나 이를 위해서 나는 단지 가장 중요한(필수적인) 몇 개의 argument 만을 언급할 것입니다.
첫 번째 argument 는 search 에 사용되는 blast program입니다. Option 과 program 에 대한 서술은 http://www.ncbi.nlm.nih.gov/BLAST/blast_program.html 에서 볼 수 있습니다. 두 번째 argument 는 검색을 할 database를 지정합니다. The second argument specifies the databases to search against. 이를 위한 option 은 NCBI web pages http://www.ncbi.nlm.nih.gov/BLAST/blast_databases.html 에서 볼 수 있습니다.
Search option을 설정한 다음에는, 당신의 fasta sequence를 세 번째 argument인 string 으로 넣어야 합니다. 그리고 나면 BLAST를 돌리기 위한 모든 준비가 된 것입니다. Biopython 은 result를 얻어서 돌려줄수 있을때까지 pause 됩니다.
Result를 parsing 하기 전에, 때로 그것들을 나중에 다시 돌아가서 모든 것을 또한번 blast로 돌리지 않고 사용할 수 있도록 file 로 저장하는 것이 유용합니다. BLAST file 로부터 정보를 extract 하는 code를 debugging 할 때 이것은 특히 유용합니다. 그러나, 단지 당신이 한 것에 대해 backup을 만드는 것으로도 역시 유용할 수 있습니다.
Return 된 result를 save하고 parsing 하길 원한다면, result를 file로 저장하기 위해서 handle.read()를 사용하는 것은 handle 안에 있는 정보는 더 이상 가지지 못하기 때문에, 주의가 필요합니다.
- 역자주 : read를 하고 나면 handle 안에 있는 정보는 없어져 버립니다. 다른 변수에다가 다시 한번 똑같은 read를 해 보시면 압니다.
이것을 handle 하는 가장 좋은 방법은 처음에 read()를 사용해서 handle 로부터의 모든 정보를 string 으로 얻고, 그 다음에는 이 string을 다루는 것입니다. (주의 : 전에는 그렇게 쓰여졌음에도 불구하고 (으악!) handle을 copy 하는 것은 가능하지 않기 때문에 copy module 은 작동하지 않을 것입니다.)
우리는 아래와 같이 함으로써 BLAST result에서 string을 얻어서 그것을 file 로 쓸 수 있습니다.
# save the results for later, in case we want to look at it
save_file = open('my_blast.out', 'w')
blast_results = result_handle.read()
save_file.write(blast_results)
save_file.close()이렇게 한 다음에, result 는 my_blast.out 이라는 file 에 저장되고, blast_results 라는 변수는 blast result를 string form 으로 간직합니다.
- 역자주 : 위의 예제들을 그대로 따라 왔다면 여기서는 blast_results = result_handle.read() 라는 line을 blast_results = b_results.read() 로 바꿔줘야 제대로 동작합니다. 이때 blast_results 는 HTML 형식의 문서가 됩니다. 그리고 read method를 사용하면 handle에 있는 정보는 없어지기 때문에 만약 똑같은 blast_results = b_results.read()를 한번 더 쓴다면 공백만 들어가게 됩니다.)
앗.. 그랬군요.. 어쩐지... 왜이렇게 코딩이 되어있을까.. 쩝.. --yong27
BLAST parser 의 parse function 은 3.1.2 에 서술된 것과 마찬가지로, parse 되어야 할 file-handle-like object를 받습니다.
우리는 handle-like object를 python standard library module 인 cStringIO를 사용하여 우리의 BLAST result 의 string 으로부터 얻을 수 있습니다. 아래의 code 는 바로 parser 로 보낼 수 있는 handle을 우리에게 줍니다.
import cStringIO string_result_handle = cStringIO.StringIO(blast_results)
이제 우리는 BLAST result를 얻었습니다. 우리는 그것을 가지고 뭔가를 하려는 준비가 되었고, 그래서 우리는 parsing section 으로 넘어갑니다.
3.1.2 WWW version of BLAST 으로부터의 결과물의 parsing
WWW version 의 BLAST는 지속적으로 결과물을 표현하는 방법을 변경하는 것으로 악명높습니다. 이것은 parser 가 그것을 처리하지 못하도록 합니다.
이 때문에, Biopython 과 같은 centralized parser를 갖는 것이 매우 중요합니다. 왜냐하면, 그렇게 되면 많은 사람들이 test 하고 사용하여 그것을 up to date 하게 유지할 수 있기 때문입니다. 이것은 혼자서 이러한 tool을 만들려고 애쓰는 것보다 더 나을 것입니다.
당신은 WWW BLAST parser 로 parsing 할 정보를 두가지 방법 중 하나에서 얻을 수 있습니다.
첫번째는 blast를 다루는 biopython function을 사용해서 정보를 얻는 것입니다. 이것은 section 3.1.1 에 서술되어 있습니다. 두 번째는 수동으로 NCBI site에서 BLAST를 하고, 그 결과를 save 하는 것입니다. 결과물을 원하는 format 으로 얻기 위해서, final BLAST page (당신이 원하는 결과를 담고 있는 바로 그 페이지) 를 HTML 로 저장해야 합니다. (Text 로 save 하면 안됩니다.) 중요한 것은 그것을 parsing 하기 위해서 data를 fetching 하는 biopython script를 사용할 필요가 없다는 것입니다.
다음의 두가지 중 한가지로 하면 됩니다. 그 다음에는 result 에 대한 handle을 필요로 합니다. Doing things either of these two ways, you then need to get a handle to the results. Python에서 handle 은 어떠한 정보의 source라도 기술하는 아주 멋진 일반적인 방법입니다. 따라서 정보는 read() 와 readline() 함수를 사용해서 retrieval 될 수 있습니다. 이것이 BLAST parser 의 input 형식입니다. (또한 다른 biopython parser 가 취하는 방식이기도 합니다.)
만약 script를 통해서 BLAST 와 interacting 하는 위의 code를 따라 왔다면, 당신은 이미 blast_results를 가지고 있을 겁니다. 그것은 parser 에 건네주기에 적합합니다. 만약 result를 file 로 저장했다면, 아래와 같은 code를 써서 file을 열어줌으로써 원하는 것을 얻을 수 있을 것입니다.
blast_results = open('my_blast_results', 'r')
- 역자주 : 여기서 my_blast_results 는 아까 기록한 파일명으로 위에서는 "my_blast.out" 으로 나왔던 것임)
이제 우리는 handle을 가졌습니다. 우리는 결과물을 paring 할 준비가 되었습니다. Parsing 할 code 자체는 아주 작습니다.
from Bio.Blast import NCBIWWW b_parser = NCBIWWW.BlastParser() b_record = b_parser.parse(b_results)
- 역자주 : 위의 blast_results 가 여기서는 b_results 로 다시 바뀌었습니다. 여기 HTML 문서인 blast_results 가 들어가면 안됩니다. 핸들형식인 b_results 가 들어가야 합니다. 만약 위에서 blast_results = b_results.read() 를 해 주셨다면 역시 에러납니다. 왜냐하면, b_results 가 이미 비어 버렸기 때문입니다.
여기서 에러 나는 경우 보충 사항
한동안 이부분을 뚫어져라 봤지요.. 몰랐던 사실을 알아낸건 file like object(정확하게는 handle. file 과 완전히 동일하게 다룰 수 있으나 메모리에 있음)를 read()하고나면, 더이상 불러올 수 없다는것... 아시겠지만 file like object는 open('filename', 'r')명령으로 추출되거나 String을 StringIO.StringIO('aaaa')해도 같은 성질의 object가 생성됩니다. 위의 b_parser.parser()는 괄호안 인수로 file like object를 받습니다.
문제가 되는 부분은 Bio 루트디렉토리의 ParserSupport.py의 safe_readline 함수입니다.
def safe_readline(handle):
"""safe_readline(handle) -> line
Read a line from an UndoHandle and return it. If there are no more
lines to read, I will raise a SyntaxError.
"""
line = handle.readline()
if not line:
raise SyntaxError, "Unexpected end of stream."
return line에서 file like object를 readline()하는데, 더이상 line이 없으면 SyntaxError를 발생합니다. 이게 이해가 안되는 것이, 파일의 끝에 도달하면 line이 없는 것은 당연할텐데 왜 이걸 에러처리했을까요... 파싱명령을 내릴때 한참시간이 걸린것으로보아, 텍스트내용(blastoup 결과)를 읽고 끝부분에서 에러출력되는 것같습니다. 무엇이 문제인지...
쩝.. 지금막 BioPython 버젼 a3로 업그래이드하고나서 해보니깐 에러안납니다. ;
이것이 수행하는 것은 BLAST를 parsing 해서 output 을 Blast Record class 로 전환하는 것입니다.
아래에 있는 Record class 에 대한 section 은 이 class 에 어떤 것이 포함되는지를 기술합니다. 하지만 일반적으로 당신이 결과물로부터 추출하기를 원할 수 있는 모든 것을 가지고 있습니다. 바로 지금, 우리는 BLAST report 로부터 어떻게 정보를 얻는지에 대한 예를 보여줄 것입니다. 하지만 당신이 뭔가 특별한 것을 원한다면, 그것은 여기서 서술되지 않을 것입니다. Record class 에 대한 정보를 자세히 보고, code 나 자동적으로 생성된 documentation을 힐끗 보십시오. Docstrings 은 정보들 각각이 무엇을 저장하고 있는지를 보여주는 많은 유용한 정보를 가지고 있습니다.
우리의 예를 계속하기 위해서, 어떤 threshold 보다 큰 blast report 의 모든 hit 에 대한 요약된 정보를 출력하도록 합시다.
아래의 code 가 그 일을 합니다.
E_VALUE_THRESH = 0.04
for alignment in b_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print '****Alignment****'
print 'sequence:', alignment.title
print 'length:', alignment.length
print 'e value:', hsp.expect
print hsp.query[0:75] + '...'
print hsp.match[0:75] + '...'
print hsp.sbjct[0:75] + '...'이것은 아래와 같은 summary report 를 출력할 것입니다.
****Alignment**** sequence: >gb|AF283004.1|AF283004 Arabidopsis thaliana cold acclimation protein WCOR413-like protein alpha form mRNA, complete cds length: 783 e value: 0.034 tacttgttgatattggatcgaacaaactggagaaccaacatgctcacgtcacttttagtcccttacatattcctc... ||||||||| | ||||||||||| || |||| || || |||||||| |||||| | | |||||||| ||| ||... tacttgttggtgttggatcgaaccaattggaagacgaatatgctcacatcacttctcattccttacatcttcttc...
기본적으로 당신이 그것을 parsing 한 다음에는, BLAST report 에 있는 정보를 가지고 원하는 것은 무엇이든 할 수 있습니다. 이것은 물론, 그것을 어디 쓸 것인가에 달려 있습니다, 하지만 다행스럽게도 이것이 당신이 해야 할 것을 하는 것을 시작할 수 있도록 도울 수 있을 것입니다.
3.1.3 BLAST record class
BLAST report 로부터 정보를 추출하는데 중요한 고려대상은 정보가 저장되어 있는 object 의 type 입니다. Biopython 에서는 parser 는 어떤 것을 parsing 하는지에 따라 Blast 혹은 PSIBlast 의 Record object를 돌려줍니다. 이 object 는 Bio.Blast.Record 로 정의되고 아주 완벽합니다.
여기에 Blast 와 PSIBlast record class를 위한 UML class diagram 에서의 나의 시도가 있습니다. 만약 당신이 UML 에 능하고 오류나 개선할 점을 발견한다면 알려 주세요. Blast class diagram 은 Figure 3.1.3 에 있습니다.
PSIBlast 의 record object 도 유사합니다. 하지만, PSIBlast 의 iteration step에서 사용되어지는 round 들을 위한 support를 가집니다. PSIBlast를 위한 class diagram 은 Figure 3.1.4에서 보여집니다.
3.1.4 BLAST를 내 컴퓨터에서 돌리기
만약 당신이 sequence를 검색할 당신 자신만의 database를 만들고자 한다면, local BLAST (StandAloneBlast) 는 절대적으로 당신이 가야 하는 길입니다.
Online BLAST와 같이, Biopython 은 당신의 script로 local BLAST를 실행시킬 수 있는 많은 멋진 code를 제공하고, 이러한 프로그램들이 제공하는 많은 command line option을 완전히 access 할 수 있도록 합니다.
많은 platform에 대해 precompile 된 local BLAST를 ftp://ncbi.nlm.nih.gov/blast/executables/ 에서 구할 수 있으며, NCBI toolbox, ftp://ncbi.nlm.nih.gov/toolbox/ 에서 당신 스스로 compile 할수도 있습니다.
특히 그 이름이 의미하듯 BLAST executable 과 대응되는 blassall 과 blastpgp 에 대한 local BLAST를 다루는 code 는 Bio.Blast.NCBIStandalone에서 찾을 수 있습니다.
Local database 에 대해서 blastall을 돌리기 위해서 이러한 기능들을 사용하여 결과를 얻어 봅니다. 첫째, 우리는 blast를 하기 위해서 필요한 모든 것에 대한 path를 설정하고자 합니다. 우리가 알아야 하는 것은 검색할 database (이것은 formatdb를 사용해서 준비되어 있어야 합니다) 에 대한 path, 찾을 file 에 대한 path와 blastall executable 에 대한 path입니다.
import os my_blast_db = os.path.join(os.getcwd(), 'at-est', 'a_cds-10-7.fasta') my_blast_file = os.path.join(os.getcwd(), 'at-est', 'test_blast','sorghum_est-test.fasta') my_blast_exe = os.path.join(os.getcwd(), 'blast', 'blastall')
이제, 모든 것을 설정했습니다. 우리는 blast를 돌리고 결과를 수집할 준비가 되었습니다. 우리는 두줄로서 이것을 할 수 있습니다.
from Bio.Blast import NCBIStandalone blast_out, error_info = NCBIStandalone.blastall(my_blast_exe, 'blastn', my_blast_db, my_blast_file)
Local blast program 에 대한 Biopython interface는 두가지 value를 돌려준다는 것을 주목하십시오. 첫 번째는 save 되거나 parser로 pass될 준비가 된 blast output 에 대한 handle 이고, 두 번째는 blast command 에 의해서 생성된 가능한 error output 입니다.
Error 정보는 처리하기 곤란할 수 있습니다. 왜냐하면 만약 error_info.read() 명령을 시도할 때, 반환되는 error info는 없습니다. 그 다음에 read() call 은 block 되고 반환되지 않지 않아서, 당신의 script 를 멈추게 할 것입니다.
내 의견으로 error를 처리하는 가장 좋은 방법은, 만약 당신이 parse될 blast_out result를 얻고 있을 때가 아니라면 단지 그것을 print 하는 것입니다. 그러나 그렇지 않다면 내버려 두어야 합니다. (역자주 : 해석상에 약간 애매함이 있음.)
만약 당신이 result를 parsing 하기 전에 file 로 저장하는데 관심이 있다면, 이를 하기 위해서 copy module을 사용하는 방법에 대한 위에 있는 WWW blast section의 정보를 보세요.
이제 당신이 output을 얻었다면, 이를 parsing 해야 합니다. 그러므로 local BLAST output의 parsing 에 대해서 더 읽으세요.
3.1.5 Local BLAST 로부터 BLAST output 의 parsing
Local BLAST 에 의해 생성된 output은 web-based BLAST 에 의해 생성된 것과는 다릅니다. 우리는 이 결과를 다루고 parsing 하기 위해서 Bio.Blast.NCBIStandalone 에 위치한 parser를 사용해야 합니다.
WWW blast 와 같이 (바로 위에 있는 정보를 보세요) 우리는 parser 에 넘겨줄 수 있는 handle object를 필요로 합니다. Handle 은 readline() method를 제대로 실행해야 합니다.
이러한 handle을 얻는 일반적인 방법은 local blast를 돌리기 위해서 제공된 blastall 또는 blastpgp 기능을 사용하는 것과, commandline 으로 local blast를 실행하여, 아래와 같이 하는 둘 중의 하나입니다.
blast_out = open('my_file_of_blast_output', 'r')그러므로, WWW blast와 같이 당신이 원하지 않는다면 Biopython 의 편리한 기능을 사용할 필요가 없습니다.
이제, handle을 얻었다면 (이것을 blast_out 이라 부를 것입니다.) 우리는 그것을 parsing 할 준비가 되었습니다. 이는 아래의 code 로 가능합니다.
from Bio.Blast import NCBIStandalone b_parser = NCBIStandalone.BlastParser() b_record = b_parser.parse(blast_out)
이는 BLAST report를 parsing 해서 Blast Record class 로 만들어서 (무엇을 parsing 하는가에 따라 Blast 또는 PSIBlast record) 그것으로부터 정보를 추출할 수 있습니다. 우리의 경우에, 어떤 threshold value 보다 큰 모든 alignment 의 quick summary를 print out 하도록 해 봅시다.
E_VALUE_THRESH = 0.04
for alignment in b_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print '****Alignment****'
print 'sequence:', alignment.title
print 'length:', alignment.length
print 'e value:', hsp.expect
print hsp.query[0:75] + '...'
print hsp.match[0:75] + '...'
print hsp.sbjct[0:75] + '...'만약 당신이 이미 WWW version 의 BLAST parsing 부분을 읽었다면, 위의 code 가 그 section에서 본 것과 같다는 것을 알 수 있을 것입니다. 한번 어떤 것을 record class 로 parsing 했다면, parsing 한 original BLAST info 의 format 과는 별개로 그것을 다룰 수 있습니다. 아주 매력적이죠!
물론, 하나의 record를 parsing 하는 것도 훌륭하지만, 많은 record를 가진 BLAST file을 가지고 있다면, 그것들을 모두 어떻게 pasring 할 수 있을까요? 떨지 마세요. 답은 바로 다음 section 에 있으니까요.
3.1.6 BLAST runs 이 가득한 file 을 parsing 하기
물론, 많은 sequence 들을 한 database 에 대해서 돌려서 그 모든 것에 대한 report를 얻을 수 있기 때문에 local blast 는 멋집니다.
그래서, Biopython 은 memory problem 없이 humungous file 들을 절대적으로 쉽게 parsing 할 수 있도록 하는 기능을 가지고 있습니다.
우리는 이것을 blast iterator를 사용하여 할 수 있습니다. Iterator를 설정하기 위해서 우리는 먼저 parser를 설정해야 하고, 이는 우리의 Blast Record object 에 있는 blast report를 parsing 하기 위해서입니다.
from Bio.Blast import NCBIStandalone b_parser = NCBIStandalone.BlastParser()
그 다음에 우리는 blast record 들을 지닌 handle 이 있다고 가정하고, 이를 blast_out 이라 부릅니다. Handle을 얻는 것은 위의 blast parsing section에서 자세히 서술되었습니다.
이제 우리는 parser 와 handle을 가지고 있습니다. 우리는 다음과 같은 command로 iterator를 설정할 준비가 되었습니다.
b_iterator = NCBIStandalone.Iterator(blast_out, b_parser)
두 번째 옵션인 paser 는 optional 합니다. 만약 우리가 parser를 제공하지 않으면, iterator 는 단지 raw BLAST report를 매번 되돌려 줄 것입니다.
이제, 우리는 iterator를 가지고 있고, 우리는 next()를 이용해서 우리의 parser 에 의해 생성된 blast records retrieving을 시작합니다.
b_record = b_iterator.next()
next 에 의한 각각의 call 은 우리가 처리할 새로운 record를 돌려줄 것입니다.
이제 우리는 이 record를 통해서 iteration을 할 수 있고, 우리가 익히 보던 멋진 Blast report를 생성할 수 있습니다.
while 1:
b_record = b_iterator.next()
if b_record is None:
break
E_VALUE_THRESH = 0.04
for alignment in b_record.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print '****Alignment****'
print 'sequence:', alignment.title
print 'length:', alignment.length
print 'e value:', hsp.expect
if len(hsp.query) > 75:
dots = '...'
else:
dots = ''
print hsp.query[0:75] + dots
print hsp.match[0:75] + dots
print hsp.sbjct[0:75] + dotsParsing 할 record 가 떨어졌을 때 b_iterator.next() 가 None을 돌려준다는 것을 주의해서 보십시오. 따라서 record 의 존재를 check 하는 while loop를 가지고 전체 file을 통한 iteration 이 일어나는 것은 쉬운 일입니다.
한번에 하나씩 읽기 때문에 iterator 는 memory problem 없이 huge blast record를 다룰 수 있도록 해 줍니다. 나는 어마어마한 file을 이것을 이용해서 아무 문제 없이 parsing 했습니다.
3.1.7 Huge file 어딘가에 있는 bad record를 찾기
내게 일어난 아주 더러븐 문제는 커다란 blast file을 paring 하는 동안 parser 가 SyntaxError를 일으킬 수 있다는 것입니다. 이것은 심각한 문제입니다. 왜냐하면 당신은 SyntaxError 가 parser 의 문제인지, BLAST 의 문제인지 알 수 없기 때문입니다. 일을 더 꼬이게 하는 것은 어디서 parsing 이 실패했는지 알 수 없다는 것입니다. 그래서 이것이 important data point 일 수 있기 때문에 error를 무시해 버릴 수가 없습니다.
우리는 이 문제를 해결하는 작은 script를 만들어야 했습니다. 그러나 Bio.Blast module 은 이제 이 문제를 더 쉽게 만드는 BlastErrorParser를 포함하고 있습니다. BlastErrorParser 는 regular BlastParser 와 매우 유사하게 작동합니다. 그러나 parser 에 의해 생성되는 SyntaxError를 catch 하는 extra layer를 가짐으로써, error를 진단할 수 있습니다.
이 parser 에 대해서 살펴보기로 합시다. 먼저, 우리는 parsing 할 file 과 problem report를 저장할 file을 정의합니다.
import os b_file = os.path.join(os.getcwd(), 'blast_out', 'big_blast.out') error_file = os.path.join(os.getcwd(), 'blast_out', 'big_blast.problems')
이제, 우리는 BlastErrorParser를 얻고자 합니다.
from Bio.Blast import NCBIStandalone error_handle = open(error_file, 'w') b_error_parser = NCBIStandalone.BlastErrorParser(error_handle)
Parser 가 handle 의 optional argument를 취한다는 것을 주목하십시오. 만약 handle 이 pass 되면, parser 는 SyntaxError를 일으키는 어떤 blast record 도 이 handle 에 기록할 것입니다. 그렇지 않으면, 이러한 record 들을 기록되지 않을 것입니다.
이제 우리는 BlastErrorParser를 regular blast parser 와 같이 사용할 수 있습니다. 특별히, 우리는 우리의 blast record를 한번에 하나씩 error parser 로 parsing 하는 iterator를 만들기를 원할 수도 있습니다.
blast_out = open(b_file) iterator = NCBIStandalone.Iterator(blast_out, b_error_parser)
앞에서와 마찬가지로, 우리는 이러한 record 들을 한번에 하나씩 읽을 수 있습니다. 그러나 지금은 우리는 Blast 에 의한 (parser 그 자체에 의한 것이 아닌) problem 으로 인한 error를 catch 하고 처리할 수 있습니다.
try:
next_record = iterator.next()
except NCBIStandalone.LowQualityBlastError, info:
print "LowQualityBlastError detected in id %s" % info[1]바로 이제, BlastErrorParser 는 아래의 errors 들을 생성할 수 있습니다.
SyntaxError -- 이것은 regular BlastParser 에 의해 생성된 것과 같은 error이고, parser 가 특정한 file을 parsing 할 수 없기 때문입니다. 이것은 보통 parser 의 bug 이거나, 사용하는 BLAST 와 parser 가 다룰 수 있는 version 의 불일치 때문입니다.
LowQualityBlastError -- 정말 질나쁜 sequence (예를 들면, 기본적으로 한 nucleotide가 반복되는 짧은 sequence)를 BLASTing 할 때, Blast 가 전체 sequence를 masking out 해서 parsing 할 아무것도 없이 끝납니다. 이러한 경우에 그것은 불완전한 report를 생성하고, parser 가 SyntaxError를 일으키도록 합니다. LowQualityBlastError 가 이러한 경우에 report 됩니다. 이 error 는 아래의 정보를 돌려줍니다.
item[0] -- error message item[1] -- error를 일으킨 Input record 의 id
이것은 만약 당신이 problem을 일으키는 모든 record를 기록하기를 원한다면 정말 유용합니다. 언급된 바와 같이, 각각의 error 가 생성될 때, BlastErrorParser 는 특정한 error_handle 에 문제를 일으킨 record를 기록합니다.
당신은 더 나아가서 이것들을 보고 다룰 수 있습니다. (역자주 : 해석상의 어려움이 있음.)
당신은 single blast report를 가지고 parser를 debug 할 수도 있고, blast run에서 problem을 발견할 수도 있습니다. 어느쪽이건 절대적으로 유용한 경험이 될 것입니다. 다행히, BlastErrorParser 는 large Blast file 들을 다루고 debug 하는 것을 더 쉽게 해 줄 것입니다.
3.1.8 PSIBlast 다루기
PsiBlast를 script 로 직접 다루는 것을 더 쉽게 하기 위해서 some stuff를 써야만 합니다. (예. alignment 로부터 적절한 format 으로 align file 의 output) 나는 PsiBlast를 더 보고, 이를 할 수 있는 몇가지 좋은 방법들을 생각해야 합니다.