BioPythonFaq about [BLAST]

<<TableOfContents>>

=== BLAST 결과의 파싱 ===

'''Q.''' WwwBlast를 돌리고 나서 얻어지는 결과를 파싱하면, GenBank ID가 안나옵니다. 결과 HTML화일을 들춰보니, 아래와 같이 GenBank부분이 href 링크부분 list_uids라는 곳에 숨어있습니다. 따라서 이내용이 파싱이 안됩니다. GenBank ID 항목을 따로 파싱하는 파서를 만들어야 할까요? 해결안을 아시는분?
{{{
<a href="http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?
cmd=Retrieve&db=Protein&list_uids=00166394&dopt=GenPept" >gb|AAA32660.1|</a>  
(M32884) leghemoglobin [Medicago sativa]          <a href = #166394> 81</a>  2e-15
}}}

'''A1.''' Biopython 으로는 해결방법이 없을듯 합니다. 처리를 안해주는 걸 어떻게 하겠습니까? 아래와 같은 별도의 code 로 list 를 만들어내는 수 밖에 없을 것 같습니다.
{{{
gi_list = []
b_results = open('xxx.html', 'r')
for line in b_results.readlines():
    if line[:92]=='<a href="http://www.ncbi.nlm.nih.gov:80/entrez/query.fcgi?cmd=Retrieve&db=Protein&list_uids=':
        gi_list.append(line[92:100])
}}}

'''A2. ''' 이문제로 고심하다가 StandAloneBlast로 nr돌려서 결과봤더니만, 여기도 그 id가 없더군요... 다시 WwwBlast로 돌아와서, 위 방법을 생각안해본건 아닌데, E value도 파싱을 해야하잖아요. 너무 복잡해진다 싶어서 제가 생각한 방법은 description.title을 쓰는 것이였습니다. 거기도 gi가 적혀있는 것이 있는 가하면, 안적혀있는 것도 있고... (걍 안적힌것 무시하고) 위의 쿼리내의 gi와는 다르지만 그래도 그걸로 NCBIDictionary하는데 이상없길래 그냥 그값을 썼었습니다. 참고로 [BioPythonTip/GenBank] 페이지에다 파싱되는 객체, 속성에 대해 정리하겠습니다.

=== BLASTP 돌린 결과로 Nucleotide record 를 얻을 수 있는가? ===
'''Q.''' BLASTP 로 protein sequence 에 대해서 돌릴 경우 위와 같이 GenPept DB 에서 받아오게 됩니다. 그런데 Nucleotide DB 에서 받아오는 방법이 있습니까? (Programming Assignment 2-2)

'''A.''' 방법이라면, GenPept DB에서 받은 내용을 파싱해서 핵산 gi를 뽑아가지고, 그걸로 NCBIDictionary를 돌리는 수밖에 없지 않겠습니까? 참고로 [[yong27]]은 위 PA2-2를 해결하는데 있어서, 엊그제 패치된 GenBank모듈을 써서 단백질 DB를 파싱했습니다. 단백질서열만 있으면 되기에 굳이 핵산 서열을 얻어오진 않았습니다. 
  
 ''형. 여기서 Programming Assignment 2-2가 뭐에요? list 가 있다면 링크 부탁합니다.''--[destine]
  예전 BTC 숙제였었슴. 문제를 찾아보려고, BioinformaticsInformation사이트에 가봤는데, 링크가 없어진듯... --[yong27], 2005-05-27


=== DB 가 틀리면 GI 도 틀린가? ===
'''Q.''' Peptide record parsing 할때, GI가 Protein DB 에서 쓰는 것과 Nucleotide DB 하고 틀린가요? 

'''A.''' NCBIDictionary(database='Protein')으로 해놓고 하니깐 그 gi로 단백질을 가져오던데요. 에러는 없었습니다.

=== Blastp 에러 문제 ===
'''Q.''' Blastp 를 돌릴때, "TypeError: popen3, argument 1: expected string, list found" 와 같은 에러가 납니다. BioPythonTip 에 나온대로 고쳐도 안됩니다.

'''A.''' formatdb 할때 protein 과 nucleotide 에 따라서 인수가 틀립니다. 이것이 제대로 지정이 안되어 위와 같은 경우가 생길 수 있습니다. {{{
formatdb -i xxx.fasta -p F -o T
--->
formatdb -i xxx.fasta -o T
}}}
이렇게 바꿔줘야 합니다. 

=== Pairwise Alignment 문제 ===
'''Q.''' Biopython 에서 Blast 를 access 하는 방법은 파일이름을 넘겨주는 식입니다. Fasta 포맷의 sequence 들의 집합이 하나의 파일로 저장되어 있을때, 가능한 모든 조합에 대해서 Pairwise Alignment 를 하려면, 따로 file 들을 만들어야 합니까? 만약 그렇다면 색인을 하나 하나의 file 별로 만들어 주어야 합니까? (Programming assignment 2. 3-2)

'''A.''' [[yong27]]이 사용한 방법입니다.
{{{
def exeBl2seq(combTup, giDict):
    sys.stdout = open('test1.fasta', 'w')
    print giDict[combTup[0]]
    sys.stdout.close()
    sys.stdout = open('test2.fasta', 'w')
    print giDict[combTup[1]]
    sys.stdout.close()

    parm = ' -i test1.fasta -j test2.fasta -p blastp'
    r, w = popen2.popen4('C:\ncbi\blast\bl2seq' + parm)
    w.close()
    return parseScore(r)
}}}
giDict는 사전형태로 해당 key에 gi, value에 fasta 객체로 저장한다음에, combTup에는 한쌍씩 모든 조합의 gi쌍이 하나씩 위 함수의 인수로 전달되고... 그래서 위 함수는 해당 쌍에 대해 각각의 파일을 만들고 bl2seq를 돌리게 됩니다. 그 결과는 parseScore라는 함수로 넘어가서 스코어를 파싱한 후 돌려주도록 했습니다. '''바뜨''' 위 코드가 실행되는데 문제가 있는지 윈도에서'''bl2seq 치명적에러'''하고 에러창이 뜬후, 프로세스가 멈춰있슴.. 아직 에러에서 헤메고 있습니다. 

'''A. ''' 두개의 for 문으로 감싸서 안쪽의 for 문이 하나 앞서가도록 하여, 전체 pair 를 만든뒤 i,j 를 파일 이름에 덧붙여서, alignment 각각을 하나씩의 파일로 만들수 있었습니다. 위의 에러는 정확히는 모르겠으나, 밑의 코드와 같이 각 file 의 full path 를 사용하지 않아서 생길 수도 있는 것 같습니다.
{{{
for i in range(len(record_list)):
    for j in range(i+1, len(record_list)):
        temp1 = open('d:\\blast\\temp1.fasta', 'w')
        temp2 = open('d:\\blast\\temp2.fasta', 'w')
        temp1.write(str(record_list[i]))
        temp2.write(str(record_list[j]))
        temp1.close()
        temp2.close()

        command = "d:\\blast\\bl2seq -i d:\\blast\\temp1.fasta -j d:\\blast\\temp2.fasta -p blastp -o d:\\blast\\pair\\pairalignment" + str(i) + str(j) + ".txt"
        os.system(command)
}}}

'''Q. ''' bl2seq 프로그램으로 alignment 한 결과를 biopython 으로 parsing 할 수 있습니까? blastall 이나 blastpgp 는 되지만, bl2seq 로 한 것은 안되는 것 같습니다. NCBIStandalone.py 에는 다음과 같이 나와 있습니다.
{{{
This module provides code to work with the standalone version of BLAST, either blastall or blastpgp, provided by the NCBI.
}}}