블라스트 XML에서의 검색

서열정렬 프로그램중 많이 알려진 프로그램이 BLAST가 있으며, HMMER 정도가 있습니다. 

HMMER은 은닉마르코프 모델 알고리즘을 사용하고 있고, 리눅스 기반이나 윈도우에서 콘솔에서 동작이 되기 때문에

처음 바이오인포매틱스로 서열정렬을 하는 입장에서는 접근이 편하지는 않습니다. 

 

여기서는 바이오 파이썬에서 Bio.SearchIO 를 활용한 간단한 서열정렬에 대한 명령어를 살펴 보는 내용을 진행해 봅니다. 

XML 형식으로 만들어진 파일로 분석하겠습니다. 

 

총 50개 데이터이며 id, def, accession 으로 데이터를 간단하게 나타냅니다. 

HSP 는 hig scoring pair의 약자로  a local alignment with no gaps that achieves one of the highest alignment scores in a given search. 으로 정의됩니다. 

 

 

from Bio import SearchIO
blast_qresult = SearchIO.read("my_blast.xml", "blast-xml")
print(blast_qresult)

 

my_blast.xml

0.13MB

결과: 

 

Program: blastn (2.12.0+)

Query: gi|8332116|gb|BE037100.1|BE037100 (1111)

MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold acc...

Target: nt

Hits: ---- ----- ---------------------------------------------------------- #

       # HSP ID + description

        ---- ----- ----------------------------------------------------------

       0 1 gi|1219041180|ref|XM_021875076.1| PREDICTED: Chenopodi...

       1 1 gi|1226796956|ref|XM_021992092.1| PREDICTED: Spinacia ...

       2 1 gi|731339628|ref|XM_010682658.1| PREDICTED: Beta vulga...

       ....

      49 1 gi|1216289774|ref|XM_021774063.1| PREDICTED: Manihot e... 

 

결과를 보면 50개의 서열결과가 나옵니다. 

이것을 psl 파일로 읽어서 출력해 보면 다음과 같습니다. 

 

blat_qresult = SearchIO.read("my_blat.psl", "blat-psl")
print(blat_qresult)

 

 

my_blat.psl

0.00MB

결과: 

Program: blat (<unknown version>)

Query: mystery_seq (61)

<unknown description>

Target: <unknown target>

       Hits: ---- ----- ----------------------------------------------------------

                  # # HSP ID + description ---- ----- ----------------------------------------------------------

                   0 17 chr19 <unknown description>

 

 

블라스트 쿼리 버전

print("%s %s" % (blast_qresult.program, blast_qresult.version))

 

결과:

blastn 2.12.0+

 

 

print("%s %s" % (blat_qresult.program, blat_qresult.version))

결과:

blat <unknown version>

 

 

blast_qresult.param_evalue_threshold

결과:

10.0

 

전체 파일 갯수 

len(blast_qresult)

결과:

50

 

len(blat_qresult)

결과:

1

 

blast_qresult[0]

결과:

Hit(id='gi|1219041180|ref|XM_021875076.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)

 

blast_slice = blast_qresult[:3] 

print(blast_slice)

 

결과:

Program: blastn (2.12.0+)

   Query: gi|8332116|gb|BE037100.1|BE037100 (1111)

 

 

블라스트 검색

blast_qresult["gi|1219041180|ref|XM_021875076.1|"]

결과: 

Hit(id='gi|1219041180|ref|XM_021875076.1|', query_id='gi|8332116|gb|BE037100.1|BE037100', 1 hsps)

 

 

 

blast_qresult.hits

 

결과: 

 

blast_qresult.hit_keys

 

결과:

 

블라스트 탐색 

"gi|1219041180|ref|XM_021875076.1|" in blast_qresult

결과
True

 

탐색

blast_qresult.index("gi|1389679838|ref|XM_016034586.2|")

결과

4     # 5번째 위치에 해당 파일이 있습니다. 

 

 

 

for문을 이용해서 총 5개의 파일을 찾아서 출력해 봅니다. 

 

for hit in blast_qresult[:5]: 
    print("%s %i" % (hit.id, hit.seq_len))

 

결과:

gi|1219041180|ref|XM_021875076.1| 1173

gi|1226796956|ref|XM_021992092.1| 672

gi|731339628|ref|XM_010682658.1| 847

gi|2031543140|ref|XM_041168865.1| 1020

gi|1389679838|ref|XM_016034586.2| 946

 

 

 

sort_key = lambda hit: hit.seq_len
sorted_qresult = blast_qresult.sort(key=sort_key, reverse=True, in_place=False)
for hit in sorted_qresult[:5]:
    print("%s %i" % (hit.id, hit.seq_len))

 

결과: 

gi|1585658262|ref|XM_028225077.1| 1176

gi|1219041180|ref|XM_021875076.1| 1173

gi|743838297|ref|XM_011027373.1| 1132

gi|1375871125|ref|XM_024605027.1| 1130

gi|1860377399|ref|XM_035077205.1| 1109

 

 

hsp값이 1보다 큰 서열은 총 몇개인지 찾아보는 명령어입니다. 

filter_func = lambda hit: len(hit.hsps) >= 1 
len(blast_qresult) 

 

결과:

50

 

 

 

filtered_qresult = blast_qresult.hit_filter(filter_func)
len(filtered_qresult)

 

결과:

0

 

 

 

 

 

for hit in filtered_qresult[:10]:
    print("%s %i" % (hit.id, len(hit.hsps)))

 

from Bio import SearchIO
blast_qresult = SearchIO.read("my_blast.xml", "blast-xml")
blast_hsp = blast_qresult[0][0] # first hit, first hsp
print(blast_hsp)

 

 

 

 

blast_hsp.query_range

결과:

(58, 678)

 

 

질의에 대한 서열

blast_hsp.query

결과:

SeqRecord(seq=Seq('ACAGAAAATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTGATCAATTGGC...GTA'), id='gi|8332116|gb|BE037100.1|BE037100', name='aligned query sequence', description="MP14H09 MP Mesembryanthemum crystallinum cDNA 5' similar to cold acclimation protein, mRNA sequence", dbxrefs=[])

 

응답에 대한 서열

blast_hsp.hit

결과

SeqRecord(seq=Seq('ACCGAAAATGGGCAGAGGAGTGAATTATATGGCAATGACACCTGAGCAACTAGC...TTA'), id='gi|1219041180|ref|XM_021875076.1|', name='aligned hit sequence', description='PREDICTED: Chenopodium quinoa cold-regulated 413 plasma membrane protein 2-like (LOC110697660), mRNA', dbxrefs=[])

 

 

두 서열 비교

print(blast_hsp.aln)

결과:

Alignment with 2 rows and 624 columns

ACAGAAAATGGGGAGAGAAATGAAGTACTTGGCCATGAAAACTG...GTA gi|8332116|gb|BE037100.1|BE037100

ACCGAAAATGGGCAGAGGAGTGAATTATATGGCAATGACACCTG...TTA gi|1219041180|ref|XM_021875076.1|

 

 

 

바이오파이썬에서 블라스트 검색어를 통한 간단한 서열 검색과 비교하는 명령어에 대해 살펴봤습니다. 

xml 파일안에서 검색과 서열비교를 보면서 서열비교에 대한 내용을 이해 할 수 있습니다.

 

코딩이 번거로울것 같아 아래 파이썬 소스파일도 첨부했습니다.

 

chapter_8.ipynb

0.02MB