[1] Pairwise SNP difference와 [2] pairwise SNP distance 중 어느 것이 맞는 용어일까? R에서 어떤 데이터 매트릭스의 row 사이의 거리(distance)를 계산하는 것은 흔히 있는 일이다. 내가 알기로는 이때에 특별히 단위는 정의되지 않는다. 만약 pairwise SNP distance의 단위를 basepair로 나타낸다면, [1]과 [2]는 같은 의미가 될 것이다.
한 bacterial species에 속하는 균주 수십개의 유전체를 일루미나 플랫폼으로 시퀀싱했다고 하자. 이어서 나서 균주 사이의 유전적 거리 혹은 상관관계를 알기 위해서 SNP를 이용한 후속 분석을 계획한다고 가정하자. 가장 적당한 방법은 최적의 reference genome sequence를 선정하여 이를 NCBI 등에서 입수한 뒤, 여기에 raw sequencing read를 매핑한 다음 SNP을 발굴하는 것이다. 즉 모든 샘플에 대해서 i) read mapping과 ii) reference에 대한 SNP 발굴이라는 작업을 반복해야 한다. 이러한 일련의 작업을 수행하는 매우 편리한 도구로서 나는 요즘 Torsten Seemann의
snippy(GitHub
링크)를 즐겨 사용한다. Raw sequencing read는 없고 sample의 contig sequence만 존재한다면 매릴랜드 대학에서 개발한
Harvest suite(
링크)가 적절할 것이다. Snippy도 contig를 input로 받을 수 있지만, 250 bp 길이의 single end read를 contig로부터 만들어서 mapping을 하는 간접적인 방법을 사용한다.
오늘의 글에서는 snippy를 사용한 SNP 발굴 및
샘플 사이의 비교에 관한 사항을 정리하고자 한다. 먼저 reference genome sequence에 sequencing read를 매핑하는 단계부터 검토해 보자. 다음 그림은 이해를 돕기 위하여 그린 것이다. 검정으로 표시된 reference 위에 3 종의 샘플에서 유래한 sequencing read를 매핑하였다. Reference genome coordinate를 기준으로 하여 read가 붙어서 alignment가 생성된 곳을 파랑색 상자로 표시하였다. 회색 상자는 인접한 샘플을 1:1로 비교하여 얻어지는 공통적인 alignment block을 표시한 것이다. 각 샘플마다 차이가 존재하므로 mapping이 된 곳은 당연히 서로 이어지지 않는 여러 개의 상자로 나타날 것이다. block C는 모든 샘플에서 공통적으로 존재하는 영역을 의미한다. 단, 각 파랑색 블록이 실제로 해당 샘플의 유전체 내에서도 이러한 순서로 존재한다고 보장할 수는 없다는 점을 유의해야 한다. 왜냐하면 샘플의 유전체에서 얼마든지 inversion 등이 일어날 수 있기 때문이다. 게다가 recombination에 의해서 더욱 복잡하게 문제가 꼬이는 것은 오늘의 글에서는 다루지 않는다.
Snippy는 각 샘플(sequence reads)에 대해서 각각 실행해야 한다. 다음은 하나의 샘플에 대한 분석 사례이다. 결과 파일은 --outdir 로 지정한 디렉토리 안에 저장된다. 이는 어디까지는 reference와 sample 1번 사이의 1:1 비교에 해당한다.
snippy --cpus 16 --outdir sample1 --R1 sample1_1.fastq --R2 sample2_2.fastq
그 다음으로 할 일은 무엇일까? 모든 샘플에서 계산된 SNP table을 다 불러들여서 SNP alignment를 만드는 것이다. 출력 파일의 포맷은 nexus fasta phylip maf clustalw 등이므로 이를 이용하여 phylogenetic tree를 작성하면 된다. SNP alignment를 만드는 프로그램은 snippy 패키지를 구성하는 프로그램 중 하나인
snippy-core이다. 다시 위의 그림으로 돌아가보자. 각 샘플 read가 reference 위에 매핑하는 곳은 완벽하게 겹치지 않는다. snippy-core는 공통적인 alignment, 즉 core site로부터 염기서열을 뽑아내는 것이다.
snippy-core --prefix core sample1 sample2 sample3...
Core site 안에서도 모든 샘플이 동일한 염기를 갖는 monomorphic site가 있는가 하면 일부 샘플만 다른 polymophic or variant site도 있다. snippy-core로 만들어진 core.aln에서 각 샘플 유래 서열의 길이를 측정해 보면 당연히 그 reference genome sequence의 크기에 훨씬 못미친다. 그런데 결과물 중에 core.full.aln이는 것도 보인다. 이것은 무엇인가? Reference sequence의 크기와 동일한 alignment로서 align이 되지 않은 영역, 즉 zero coverage region or deleted region을 '-'으로 채운 것이다. 위 그림을 기준으로 설명하자면 빨강 블록을 연결하는 두꺼운 회색수평선을 '-'으로 채웠음을 의미한다. 실제로 full alignment 파일을 열어보면 ATGC와 '-' 이외의 문자도 보일 것이다. 이것 때문인지는 모르겠으나 core.full.aln을 FastTree로 처리하여 tree를 그리면 간혹 이상한 그림이 나오기도 한다. core.full.aln에 수록된 다양한 문자형태에 대한 설명은 snippy 웹사이트에서 가져온 설명으로 대체한다.
- 대문자: same as the reference
- 소문자: different from the reference
- -: zero coverage in this sample or a deletion relative to the reference
- N: low coverage in this sample (based on --mincov)
- X: masked region of reference (from --mask)
- n: heterozygous site in the sample (has GT=0/1/ in smps.raw.vcf)
다음으로 알아볼 것은
snp-sites(
GitHub 논문)라는 유틸리티이다. 이것은 multi-FASTA alignment에서 SNP를 추출하는 도구로서 snippy의 구동을 위해 꼭 필요하다. 논문을 보면 snippy의 개발자인 Torsten Seemann이 공저자로 포함되어 있기도 하다. snippy-core가 만들어낸 core.aln을 snp-sites로 처리하면 monomorphic site가 제거되는 것으로 이해하면 될 것이다.
마지막으로 오늘의 주제인
snp-dists(GitHub
링크)에 이르렀다. 이 유틸리티 역시 Torsten Seemann의 작품이다. 3 개의 sample에 대한 SNP 분석을 했다는 것은 reference에 대한 sample의 SNP 분포를 알아냈다는 뜻이다. 그러면 sample 1번과 sample 2번을 비교했을 때 SNP의 분포는 어떻게 되는가? 엄밀하게 말하자면 여기에는 두 가지 방법이 있다. 첫번째는 위 그림에서 공통적으로 존재하는 alignment(맨 아래의 빨강색 박스)를 이용하는 것으로서 snp-dists도 이에 해당한다. 다시말해서 core.aln 파일에서 샘플 사이의 모든 pairwise SNP를 재계산하는 것이다.
두번째 방법은 내가 실제로 어제 사용한 방법이었는데 결론적으로 말해서 바람직하지는 않은 것으로 여겨진다. 모든 샘플에 대해서 snippy를 다 실행해 둔 다음, 결과 디렉토리를 전부 두 개씩 짝을 지어서(Perl의 Algorithm::Combinatorics를 사용함) snippy-core를 돌린 것이다. 사용한 샘플의 수가 98개이니 여기에서 2 개씩을 모두 조합하면 C(98, 2) = 4753 개의 command line에 해당하는 snippy-core를 실행한 셈이다. 얼마든지 parallel하게 실행할 수 있는 일이었지만 아주 단순하게 순차적으로 command line을 실행하면서 분석이 끝나기를 기다렸다. '과연 이렇게 하는 게 맞는 것일까?'하는 의구심과 함께... 화면으로 출력되는 standard error를 파일로 리다이렉션한 뒤 다음 라인을 전부 찾아내어 parsing을 하였다.
[14:46:48] Found 3568 core SNPs from 6612 variant sites.
내가 원하는 것은 matrix를 만드는 것이 아니었다. 98 개의 샘플(strain)은 몇 개의 clade로 나눌 수 있는데, 동일 clade에 속하는 strain 간의 SNP distance가 다른 clade에 속하는 것과의 비교에서 나온 수치보다 현격하게 적다는 것을 보이는 것이 목적이었다. 지리한 계산이 다 끝나고 요즘 쓰고 있는 논문에 해당 수치를 기입하고 나서야 snps-dists라는 유틸리티를 발견하게 되었다. 이것에 core.aln(FASTA) 파일만 input으로 제공하니 순식간에 SNP distance matrix를 작성하였다. 98개의 snippy 결과 디렉토리를 두 개씩 전부 조합한 command line을 만들어서 훨씬 많은 시간이 걸려서 4700회가 넘는 snippy-core를 실행했던 내 자신이 한심하게 느껴졌다.
그러면 snp-dists가 만든 매트릭스는 실제 어떻게 생겼나? ooffice --calc로 파일을 열어보았다. 그런데 각 셀에 적힌 수치가 snippy-core 수작업으로 계산한 것보다 약간 작다. 왜 그럴까? 그 이유를 곰곰이 생각해 보았다. 다시 위의 그림으로 돌아가자. 스크롤하기가 불편하니 아래에 다시 불러다 놓고 글을 마저 적도록 하겠다.
snp-dists는 alignment file을 입력물로 사용한다. 이는 위 그림에서 맨 아래쪽의 빨강 박스를 기준으로 작동한다는 것을 의미한다. 즉 모든 샘플에서 공통적으로 mapping된 영역에 대해서만 pairwise SNP distance를 계산한다는 것이다. 그러면 내가 어제 했던 '비능률적인 방법'은 어떻게 돌아가는 것인가? sample_1과 sample_2의 SNP 분석물 1:1 비교를 생각해 보자. 위 그림에서 두 샘플의 read alignment 중 공통적인 영역, 즉 회색으로 표시된 부분에서 SNP를 찾아낸 것이다. 당연히 빨강 박스에서 찾아낸 것보다 더 많은 SNP가 발견되었다. 그 다음으로 sample_2와 sample_3의 snippy 결과를 snippy_core로 처리한 경우를 생각해 보자. 마찬가지로 두 샘플 사이를 연결하는 회색 영역에서 SNP를 찾은 것이다.
이상에서 열거한 두 가지 방법 중에서 어느 것이 옳은가? 나도 모른다. 사용하기에 편한 방법을 찾으라면 당연히 snp-dists이다. 모든 샘플에 공통적으로 존재하는 영역에 대해서만 SNP를 추출하여 그 숫자를 계산하므로 보다 객관적일 것이라는 생각도 든다. 하지만 샘플 대 샘플의 1:1 비교라는 측면에 집중한다면 내가 실행했던 다소 비능률적인 방법도 나름대로 타당하다고 생각한다. 특정 샘플의 쌍에서만 존재하는 공통 alignment 영역(위 그림에서 초록색으로 표시한 부분)에도 주의를 기울인다는 것이 이 방법의 중요한 특징이다. 이 영역은 core.aln에서는 당연히 취급되지 않는다.
