[BioPython] BLAST

Bioinformatics

[BioPython] BLAST

rudgh99_algo 2024. 9. 18. 16:45

0. BLAST란?

BLAST란 Basic Local Alignment Search Tool의 줄임말이다. 쉽게 말해서, 내가 알고 있는 Seq와 유사한 Seq가 있는지 찾아주는 알고리즘이다. "바이오파이썬으로 만나는 생물정보학"에 좋은 예시가 있어, 여기서 설명한다. 어떤 질병이 있는 환자 A가 있다. 우리는 이 질병이 어떤 seq가 변이가 되어, 생긴 건지 알고 싶다. 따라서, 환자의 seq를 채취해 , sequencing 회사에 보냈다. 결과 파일을 받았고, Alignment 한 결과, 인간의 genome과 matching 되지 않는 부분이 있음을 발견했다. 그렇다면, 이 부분은 과연 어떤 종에서 유래된 seq란 말인가? 이때, BLAST를 이용한다. BLAST결과 streptococcus pneumoniae라는 것을 알았다. 이런 식으로 활용된다.

1. BLAST를 web에서 사용해 보기

web에 BLAST ncbi를 검색 후, blast를 실행해 볼 수 있다. 이때, blast program 이름이 여러 개가 있다.

1.1 blastn : DNA vs DNA 비교

1.2 blastp : Protein vs Protein seq 비교

1.3 blastx : 번역된 DNA 서열 vs Protein seq 비교 ( 설명: blastx는 뉴클레오타이드 서열을 6개의 가능한 아미노산 서열로 번역한 후, 이를 단백질 서열 데이터베이스와 비교합니다. 이는 DNA 서열이 단백질을 암호화하고 있을 가능성이 있을 때 사용됩니다.)

1.4 tblastn : Protein Seq vs 번역된 DNA Seq 비교 ( 설명: tblastn은 단백질 서열을 입력으로 받아, 데이터베이스 내의 DNA 서열을 6개의 가능한 아미노산 서열로 번역하여 비교합니다. 이는 단백질 서열을 알고 있을 때, 해당 단백질과 유사한 DNA 서열을 찾고자 할 때 유용합니다.)

1.5 tblastx : 번역된 DNA seq vs 번역된 DNA seq 비교 ( 설명: tblastx는 입력 DNA 서열과 데이터베이스의 DNA 서열을 모두 6개의 가능한 아미노산 서열로 번역한 후 비교합니다. 이는 두 DNA 서열 간의 유사성을 가장 포괄적으로 탐색할 때 사용됩니다.)

이중 필요한 프로그램을 사용하여, 시행할 수 있다.

2. BLAST를 BioPython으로 사용해 보기

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

record = SeqIO.read("file name", "file format") 
#handle = NCBIWWW.qblast("program name", db name, SeqRecord 객체)
## record.format("fasta"): 읽어들인 서열을 FASTA 형식의 문자열로 변환하여 전달합니다.
handle = NCBIWWW.qblast("blastn","nt",record.format("fasta")) 

blast_records = NCBIXML.parse(handle) 
E_value = 0.05 
for blast_record in blast_records: 
	for alignment in blast_record.alignments:
    	for hsp in alignment.hsps:
    		if hsp.expect < E_value:
           		print(alignment.title) 
                print(alignment.length)
                print(hsp.expect) 
                print(hsp.query[0:75]) 
                print(hsp.match[0:75]) 
                print(hsp.sbjct[0:75])

<부가설명 by chatGPT o1 mini :

alignment.title: 일치하는 데이터베이스 서열의 제목이나 설명입니다. 예를 들어, NCBI의 GenBank ID, 단백질 이름 등이 포함될 수 있습니다.
alignment.length: 일치 서열의 전체 길이입니다.
hsp.expect: HSP의 E-값으로, 이 값이 작을수록 일치가 통계적으로 유의미함을 나타냅니다.
hsp.query: 쿼리 서열에서 일치하는 부분의 서열입니다.
hsp.match: 쿼리와 데이터베이스 서열 간의 매칭 상태를 나타내는 문자열입니다. 일치하는 위치는 |로 표시됩니다.
hsp.sbjct: 데이터베이스 서열에서 일치하는 부분의 서열입니다. >

gi|123456|ref|NM_000123.4| Homo sapiens gene ABC, mRNA
1500
1e-10
ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA...
|||||||||||||||||||||||||||||||||||||||||||||||
ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA...

blast_record에 대한 좀 더 깊은 이해가 필요하다. 또한, blast가 어떤 algorithm에 의해 실행되는지 알 필요가 있다.

좀 더 파고들기

rom Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import SeqIO
records = SeqIO.read("../fasta/buccal_swab.unmapped2.fasta",format="fasta")
handle = NCBIWWW.qblast("blastn","nt",records.format("fasta"))
blast_records = NCBIXML.parse(handle)
i = 0
for blast_record in blast_records:
    i += 1
    print(f"{i}th blast_record : {blast_record.query}")
    j = 0
    for alignment in blast_record.alignments:
        j += 1
        print(f"{j}th alignment")
        k = 0
        for hsp in alignment.hsps:
            k += 1
            if hsp.expect < 0.05:
                print(f"{k}th hsp:")
                print(f"alignment_title : {alignment.title}")
				print(f"query_seq: {hsp.query[:75]}")
                print(f"matching_seq: {hsp.match[:75]}")
                print(f"sbjct_seq: {hsp.sbjct[:75]}\n")

1th blast_record : buccal_swab.unmapped2
1th alignment
1th hsp:
alignment_title : gi|2718571051|gb|PP680711.1| Mutant Human alphaherpesvirus 1 clone recombinant GFP-McKrae genomic sequence
query_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC
matching_seq: |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sbjct_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC

2th alignment
1th hsp:
alignment_title : gi|2627592510|gb|OQ077661.1| Mutant Human alphaherpesvirus 1 clone H129-Syn-G2, complete sequence
query_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC
matching_seq: |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sbjct_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC

3th alignment
1th hsp:
alignment_title : gi|2624849497|gb|OR833071.1| Human alphaherpesvirus 1 isolate HSV1/USA/WA-UW-2R10/2020
query_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC
matching_seq: |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sbjct_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC

4th alignment
1th hsp:
alignment_title : gi|2624849422|gb|OR833070.1| Human alphaherpesvirus 1 isolate HSV1/USA/WA-UW-2R6/2020
query_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC
matching_seq: |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sbjct_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC

5th alignment
1th hsp:
alignment_title : gi|2624849347|gb|OR833069.1| Human alphaherpesvirus 1 isolate HSV1/USA/WA-UW-2M3/2020
query_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC
matching_seq: |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
sbjct_seq: CCAGCCCCCCAGCCTCCCGATCACGGTTTACTACGCCGTGTTGGAGCGCGCCTGCCGCAGCGTGCTCCTAAACGC

위 결과로부터 추론할 수 있는 것이 하나 있다.

그렇다면 blast_record in blast_records에서 blast_record는 무엇인가?

내가 추측하기로는 이 홈페이지 전체다. query가 두 개 이상이라면 아마 유효한 것 같다. 따라서, 계층구조를 표현해 보면 다음과 같다.

blast_records (이터레이터)
└── blast_record (첫 번째 쿼리 서열 결과)
    ├── alignments (검색 결과 테이블의 히트 목록)
    │   ├── alignment1 (첫 번째 히트)
    │   │   ├── hsps (히트의 일치 조각 목록)
    │   │   │   └── hsp1 (첫 번째 HSP)
    │   └── alignment2 (두 번째 히트)
    └── blast_record2 (두 번째 쿼리 서열 결과)
        ├── alignments (히트 없음)
        └── ...

hsp에 대한 이해

쿼리 서열:   ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA
데이터베이스 서열: ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA

Alignment:
ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA
|||||||||||||||||||||||||||||||||||||||||||||||
ATGCGTACGTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA

HSP1:
쿼리: 1-25 (ATGCGTACGTTAGCTAGCTA)
Hit: 1000-1025

HSP2:
쿼리: 30-50 (CTAGCTAGCTAGCTA)
Hit: 1050-1070

blast_records (이터레이터)
└── blast_record (하나의 쿼리 서열에 대한 결과)
    ├── alignments (일치 서열 목록)
    │   └── alignment (하나의 데이터베이스 서열과의 일치 정보)
    │       └── hsps (일치 조각 목록)
    │           ├── hsp1 (첫 번째 일치 조각)
    │           └── hsp2 (두 번째 일치 조각)
    └── ...