2024-05-24 ~ 2024-05-28
1. History of search sequence
1) Needleman-Wunsch의 global alignment algorithm (1970)
- 그러나 전체적인 유사도보다는 보다는 지역적인 유사도가 보다 생물학적인 의미를 지니고 있기 때문에 local alignment 도구 필요
2) Smith-Waterman의 local alignment algorithm (1981)
- dynamic programming 사용
- 최적의 비용을 가지는 정렬을 찾아주지만 계산 시간이 오래 걸림
3) FASTA(1988), BLAST(1990)
- word based and heuristic algorithhms(모든 서열을 다 조사하는 것이 아님)
2. Sequence alignment
서열 정렬은 DNA, RNA 또는 단백질의 서열을 배열하여 서열 간의 기능적, 구조적 또는 진화적 연관성을 모두 파악하여 유사성을 확인하는 것이다. 일반적인 정렬방식은 아래와 같이 나눌 수 있다.
1) Global alignment vs Local alignment
- Global alignment: 서열 비교 시 각 서열의 모든 길이를 비교하여 전체적인 유사도(similarity)를 파악하는 방법이다. 같은 길이의 서열을 비교할 때 좋음.
e.g. f Needleman and Wunsch
- Local alignment: 가장 유사한 서열이 많은 부위만 비교하는 방법으로, 검색하고자 하는 서열이 데이터베이스에 비해 query 서열이 짧을 때 좋음.
e.g. Smith and Waterman algorithm/ FASTA/ BLAST
추가적으로는 dynamic programming, heuristic algorithms 혹은 probabilistic methods도 존재한다.
2) Pairwise Alignment vs Multiple Sequence Alignment
- Multiple Sequence Alignment: 여러 서열 간의 비교
- Pairwise Alignment: 두 서열 간의 비교
- dot-matrix methods
- dynamic programming
- word methods
3. Substitution matrix
1) 치환 행렬(Substitution matrix)이란?
뉴클레오타이드의 정렬에서 일치, 불일치를 단순하게 고려하면 된다. 실제로 BLAST는 match에서 +2, mismatch에서 -3의 점수를 부과한다. 그러나 아미노산의 서열 정렬에서는 물리 화학적 특성을 고려한 치환행렬을 사용해야 한다.
즉 치환행렬은 아미노산의 진화적 관계(보존과 비보존)를 고려하기 위한 점수 표이다.
아래 예시를 통해 치환행렬 내의 값이 어떻게 반영되는지 알아보자.
두 서열을 정렬하였을 때, 일치하지 않는 부분을 확인할 수 있다.
그러나 위, 아래 서열 모두 같은 단백질을 암호화하고 있다.
즉, 단순한 서열의 유사성이 아닌 단백질로 번역된 후를 고려해야 한다.
또한 우리는 단백질의 구조를 고려해야 한다. 아래의 세 가지 아미노산을 관찰해 보자.
Phenylalanine, Aspartic acid는 Tyrosine과 오직 하나의 염기 차이를 보인다.
그러나 구조적으로 유사한 Tyrosine과 Phenylalanine과 달리 Aspartic acid와는 큰 차이가 존재한다.
이러한 아미노산 사이의 관계를 점수로 나타낸 행렬이 바로 Substitution matrix (혹은 scoring matrix)이다.
2) PAM matrices
최초로 substitution matrix를 제안한 사람은 Margaret Dayhoff이다. Dayhoff와 동료들은 1978년에 scoring matrix인 Accepted point mutation(PAM) matrices을 개발했다.
Dayhoff는 단백질에 따라 진화속도가 다르다는 결과를 발표했다(e.g. 히스톤은 1억 년 동안 조금밖에 변하지 않지만, 카파 카세인의 경우 매우 빠르게 변한다).
이 데이터를 기반 삼아 만들어진 아미노산 사이의 치환 및 보존 값을 PAM matrices로 나타냈다.
PAM 행렬은 identity가 85% 이상인 단백질의 70개 그룹에서 아미노산 치환 패턴을 분석한 결과이다.
- PAM1는 100개의 아미노산에서 한 개의 아미노산이 치환될 확률 (1% divergence)
- PAM100는 100개의 아미노산 당 100번 치환이 일어날 확률(모든 아미노산이 치환될 확률이 아님, 이때 여러 번 치환으로 원래 아미노산으로 되돌아갈 수도 있음) (PAM1기반)
- PAM250는 100개의 아미노산에서 250번의 치환이 일어날 확률(PAM1기반)
PAM250
가장 대표적인 matrix는 PAM250를 살펴보자.
PAM matrix의 기본 가정은 "한 단백질 내에서의 특정 위치의 아미노산의 변화는 다른 위치의 아미노산들에 의해 영향받지 않으며, 또한 특정 위치에서 아미노산이 교체될 확률은 그의 위치와는 무관하다"는 것이다. 그러나 이러한 가정은 생물학적으로 위배되었으며, 이에 대한 한계를 가졌다고 볼 수 있다.
3) BLOSUM matrices
BLOSUM(BLOcks SUbstitution Matrix)은 1991년에 Steven Henikoff Altschul과 Jorja Henikoff에 의해 제시되었다. 이때 사용한 데이터는 BLOCKs database로 아미노산 서열 중 conserved 한 지역의 데이터를 담고 있다.
- 어떠한 진화모델에 근거하지 않고, 직접적으로 관찰에 의해 획득된 값이다.
- 초기 blast는 PAM사용하였지만, 이제는 BLOSUM matrices로 대체되었다.
BLOSUM도 PAM처럼 여러 버전에 존재하며, 적힌 숫자는 계산에 사용된 서열들의 보존(유사성)을 나타냄
- BLOSUM62 서열이 적어도 62% 동일할 때 아미노산이 쌍을 이룰 확률 계산하여 구한 행렬값
- BLOSUM90 ~ 90%.
- BLOSUM45 ~ 45%.
PAM vs BLOSUM
| PAM | BLOSUM |
| global alignments 를 기반으로 한다 | local alignments 를 기반으로 한다 |
| PAM1은 차이가 1% 이하이지만 identity 99%에 해당하는 서열의 비교로부터 계산된 행렬값이다. | BLOSUM 62는 62% 이하의 pairwise identity 을 갖는 서열의 비교로부터 계산된 행렬값이다. |
| score가 높을 수록 진화적 거리가 더 크다는 뜻 | score가 클 수록 유사성이 높고 진화적 거리가 짧다는 뜻 |
(https://en.wikipedia.org/wiki/BLOSUM#The_relationship_between_PAM_and_BLOSUM)
4. BLAST
1990년 미국 NCBI의 연구자들에 의해 개발된 유사성 검색 프로그램으로 heuristic search algorithm와 Karlin-Altschuls의 통계방법을 채택하고 있다.
- BLAST의 혁신점은 neighborhood words(확장 단어)을 도입하며 검색 속도를 높인 점이다.
- 또한 FASTA에 비해 통계적으로 엄격하고, 더 sensitive 하다.
1) Word search
우리가 알고자 하는 query 서열에서 길이 w에 해당하는 작은 조각인 word를 만든다.
- 왜? 서열을 하나씩 비교하는 것보다 훨씬 빠르기 때문이다.
- 아미노산에서 기본 word의 길이는 3, 염기에서는 11이다.
PQG에 집중해 보자. PQG 단백질 서열이 얼마나 잘 보존되는지 알아보기 위해 BLOSUM 62 치환행렬을 기반으로 살펴보면, 총합은 18이다.
이때 우리는 치환행렬 값을 기반으로 비슷한 단백질 서열들을 나열할 수 있다. 총 합(T) 값이 13 이상을 PQG와 유사한 word라고 정의하자. 이때 T 값이상인 유사한 word를 Neighborhood words라고 한다.
2) Identification of exact word match method
이제 Neighborhood words를 모두 구했으면, 이 words list를 기준으로 database 서열과 일치하는지 알아본다. word와 database가 일치되는 부분을 seed라고 한다.
3) High segment pair (HSPs) alignment method
seed가 생성되었다면, 이를 기준으로 양 옆으로 더 일치하는지 확인한다.
- similarity score(raw score)는 BLOSUM을 기반으로 각 염기들의 관계에 따라 match 혹은 mismatch점수가 부과된다.
- 전체적인 score가 최대가 되는 지점(high scoeing pairs; HSP)까지 확장한다.
만약 PQG-PMG 쌍에서 한쪽 방향으로만 유사한 서열을 확인한다고 할 때, 어디까지 확장해야 할까?
BLOSUM62를 기반으로 한 total score는 아래와 같은 추세를 보일 것이다.
이때 우리는 시작점부터, 상승 추세가 꺾이는 부분까지를 High score라고 판단한다.
(실제로는 한 방향이 아니라 양 방향으로 뻗어 나가면서 계산).
Gap penalty
blast 1.0에서는 gap을 허용하지 않지만, 2.0에서는 일부 허용한다. 이때 Gap에 대한 페널티는 어떻게 주어야 할까?
모든 gap 한 칸 마다를 부과해야 하는가? 연속된 Gap 이 많을수록 더 많은 페널티를 부과해야 할까?
아래 3개의 Gap이 존재하는 정렬을 살펴보자.
왼쪽 오른쪽 모두 3개의 Gap이 존재한다. 왼쪽은 각각 3개의 사건로 인해 Gap이 생성되었으며, 오른쪽은 하나의 사건로 인해 생성될 수 있다. 즉, 왼쪽에 더 많은 페널티를 부과해야 한다.
이를 위해서는 Gap의 생성과 연장 모두에 페널티를 부과하는 것이 바람직하다( Gap open = -11, Gap extend = -1).
실제로 Gap penalty를 고려한 similarity score를 살펴보자.
즉, HSP에 대한 raw score(S) 은 match 시 점수 + mismatch 시 점수(감점) + gap시작에 대한 페널티 + gap 연장에 대한 페널티를 모두 더한 값이다.
5. BLAST statistical significance
1) E-value
길이 m과 길이 n을 가진 두 서열이 S이상의 score를 나타낼 기댓값(Expected value; E value)을 아래와 같다. 이는 Karlin-Altschul 방정식에서 유래하였다.
즉, E value 은 query마다 우연히 S값 이상의 정렬된 서열을 찾을 경우 hit 수에 대한 기댓값이다(확률이 아니다).
- m: length of query sequence
- n: length of database sequence
- K: search space에 따른 매개변수(m과 n의 영향받음)
- λ: scoring system에 따른 매개변수
이 식을 보면 m과 n의 길이가 길 수록 HSP를 찾을 기댓값이 증가될 것이다. 반대로 HSP의 raw score(S)가 증가할수록, 높은 S 값 이상의 HSP를 찾기 어려워지니 E value는 지수적으로 감수한다.
즉, E값은 0에 가까울수록 발견된 HSP보다 더 높은 값의 HSP 점수가 나타나지 않을 것을 뜻한다.
| E-value | description |
| > 10^-1 | insignificant |
| 10^-5~10^-1 | distantly homologous |
| 10^-50 10^-5 | homologous proteins |
| 10^-100~10^-50 | nearly identical |
| <10^-100 | identical proteins |
(https://www.youtube.com/watch?v=S3gr8gjKHhc&t=2s)
만약 0.01을 선택하면 100개의 query 마다 한 번의 오류가 발생한다는 것을 의미한다. 만약 데이터베이스가 두 개라면 두 번의 에러를 얻을 수 있다. 이러한 이유로 E-value 값에서 데이터베이스의 크기가 중요한 것이다.
(https://www.researchgate.net/post/What_is_E-value_exactly_means_and_what_is_it_1e-63_representing)
2) Bit score
K, λ 값은 HSP 마다 다르기 때문에 다른 서열과의 비교를 위해 S값을 K, λ로 정규화 해준다. 이를 bit score(S')라고 한다.
bit score를 사용한 E value추정은 아래와 같다.
이를 보면 E-value값에서 중요한 변수는 m과 n에 따른 검색공간(search space)이 임을 다시 한번 확인할 수 있다.
만약 λ = 1.37, S = 23, k = 0.771일 때, bitscore의 값은 약 46이다.
또한 m이 34, n이 5,854,611,841이라고 가정할 때, E value는 34*5854611841*2^(-46) = 0.003이다.
3) P-value
HSP 점수의 통계적 유의성 평가는 Gumbel extreme value distribution (EVD)를 사용하여 평가된다.
x보다 크거나 같은 점수 S를 관찰할 확률 p는 다음 방정식으로 제공된다.
이때 μ값은 아래와 같다.
혹은 아래와 같이 표현 가능하다.
이때 P value는 결과적으로 E value와 유사하다.
Example ➕
human p53를 UniplotKP/Swiss-Prot database(C.elegans)에서 정렬결과를 살펴보자.
- length of quert = 305
- length of database = 1,414,501
- K = 0.0410
- λ = 0.267
- similarity score: 77
Q. 이때의 bis score와 E-value는?
Q. 이때 p-value를 값은?
결과적으로 E value가 작으면, P value값과 유사하게 나타난다.
참고
BLAST Statistics
- https://www.ncbi.nlm.nih.gov/BLAST/tutorial/Altschul-1.html 👍 👍 👍
- 한국어 강의 (https://www.youtube.com/watch?v=c7xqId2Sto0)
BLAST algorithm
- blast용어: https://www.ncbi.nlm.nih.gov/books/NBK62051/
- https://www.youtube.com/watch?v=wfi_KimrNQM
- MIT open course, Foundations Of Computational And Systems Biology, 2014 (https://ocw.mit.edu/courses/7-91j-foundations-of-computational-and-systems-biology-spring-2014/resources/lecture-2-local-alignment-blast-and-statistics/)
- RobEdwards. BLAST 2 seed and extend (https://www.youtube.com/watch?v=kAAme1fBanc)
+) 요약 강의: https://www.youtube.com/watch?v=g0nSH17psDc 👍 👍
- http://ibi.zju.edu.cn/bioinplant/temp/BSA/BSA_Blast_BLAST.pdf
- 유재천 외. (2022). 유전자 발굴 알고리즘 개발, 한국과학기술정보연구원
- Mount DM. (2004). Bioinformatics: Sequence and Genome Analysis. Cold Spring Harbor Laboratory Press: Cold Spring Harbor, NY. ISBN 978-0-87969-608-5.
Alignment
-Altschul SF. (1991). Amino acid substitution matrices from an information theoretic perspective. Journal of Molecular Biology. doi:10.1016/0022-2836(91)90193-A
- Weik TC, Bio Informatics (https://www.fh-muenster.de/eti/labore_forschung/db/bio-informatics.php) 👍
PAM
- Dayhoff MO, Schwartz RM, Orcutt BC. (1978). A model of Evolutionary Change in Proteins. In: Atlas of protein sequence and structure. National Biomedical Research Foundation. Washington, DC. ISBN 978-0-912466-07-1.
