cross_match를 이용한 sequence alignment의 score 이해하기

cross_match는 아주 고전적인 염기서열 정렬(alignment) 프로그램이다. 짧은 한 쌍의 염기서열을 정밀하게 비교할 때 이보다 좋은 방법은 없다고 생각한다. 사용하는 알고리즘은 word-nucleated banded Smith-Waterman comparison이다. 수 megabase가 넘어가는 (박테리아)유전체 수준의 정렬, 또는 수많은 염기서열 데이터베이스에서 원하는 서열 단편을 찾을 때에 적합한 프로그램은 따로 있다.

cross_match에서 쓰이는 score matrix는 다음과 같다. N이 아닌 동일 염기의 매치에는 1씩의 점수가 부여된다.

Score matrix (set by value of penalty: -2)
    A   C   G   T   N   X
A   1  -2  -2  -2   0  -3
C  -2   1  -2  -2   0  -3
G  -2  -2   1  -2   0  -3
T  -2  -2  -2   1   0  -3
N   0   0   0   0   0   0
X  -3  -3  -3  -3   0  -3

프라이머 서열로 쓰기에 적합한 짧은 염기서열의 매치를 찾기 위해서 cross_match를 몇 차례 실행하고 있었다. mismatch가 전혀 없는 100 bp 길이의 alignment라면 당연히 score가 100이 될 것이라고 생각했는데 결과 파일을 보니 점수는 이보다 낮았다. 왜 그럴까? cross_match 문서를 찾아보았다.

점수가 100% identical match에 대해서 alignment length보다 짧은 이유는 기본 조건에서는 complexity-adjusted score가 쓰이기 때문이었다. 염기 조성에서 편이(bias)가 존재하는 경우 alignment score는 더 낮게 조정된다는 것이다. 매우 합당한 정책이다. 만약 조정을 하지 않은 score로 표현하고 싶다면 -raw 옵션을 주면 된다. 약 20년 전에 phred/phrap/consed의 라이센스를 받아 Sanger sequencing 결과물의 de novo assembly에 처음 쓰기 시작하면서 지금까지 별로 신경을 쓰지 않던 것을 NGS 시대에 다시 들춰보려니 감회가 새롭다. Sequencing chemistry는 두어 서대 진보하였지만, 염기서열을 커맨드 라인에서 조물락거리는 것은 그때나 지금이나 별로 다르지 않다.

Banded search의 기본 조건은 -minmatch 14이고, match의 필터 과정에서는 -minscore 30이 된다. 벡터 서열을 스크리닝하는 경우에는 좀 더 민감한 검색이 이루어지도록 -minmatch 10 -minscore 20의 조건을 준다.

그러니 primer 서열(보통 30 bp 미만)을 점검하기 위해 cross_match를 쓰는 경우에는 -minmatch 파라미터를 낮게 주도록 하자.