본문 바로가기
Bioinformatics/Algorithms

GenomeMapper 색인 생성

by 임은천 2012. 11. 29.

본 문서는 내 동료 Jorg의 diploma thesis의 내용 중 일부를 번역한 것이다.


색인 구조 생성(mkindex)


이 프로그램은 각 지놈 혹은 조합된 지놈에 대해서 오직 한 번만 수행되어야 한다. 이것은 해시 색인을 생성하여 씨드의 발생과 위치의 빠른 룩업을 가능하게 하고, 참조 서열에 대해 계통([생물]strains)의 다른점들과 관련이 있는 모든 서열 정보를 포함하는 시퀀스 그래프를 생성한다. 입력으로 프로그램은 단지 FASTA 포맷으로 된 참조 서열 하나 혹은 다수의 다른 계통 서열들을 받을 수 있다. 각 계통들은 하나의 분리된 입력 파일 안에 있어야 한다. 이 프로그램은 모든 변이에 대해서 염색체, 서열의 위치, 그리고 계통적 다형성(삽입이 삭제와 SNP보다 더 선호된다)에 따라 오름차순으로 표시한다. 정확한 명세는 부록에 있는 메뉴얼에 나와 있다. 계통에서 돌연변이는 제공되는 참조 서열에 존재하는 염색체와 동일한 염색체에 대해서만 처리될 수 있다. 이 프로그램의 출력 파일은 genomemapper 매핑 프로그램의 입력 파일로 사용될 수 있다.

GenomeMapper에 의해서 사용되는 시퀀스 그래프의 구조를 설명하기 전에, 먼저 어떻게 지놈 서열이 소위 말하는 "블록들"에 저장되고 어떻게 인코딩되는지 설명한다. 다음에 그래프 생성 알고리즘이 설명되고 최종적으로 그래프로부터 모든 씨드 발생 정보로 색인 구조를 채우기 위해 사용되는 방법이 설명될 것이다.


지놈을 블록들로 쪼개기


해시 테이블은 모든 씨드들을 그것이 지놈 상에서 발견되는 위치와 모든 씨드들을 연관시켜야한다. 씨드의 위치 정보는 더 복잡하다. 이 정보는 어떤 지놈에서 씨드가 나왔는지, 어떤 계통의 염색체에서 나왔는지, 참조 서열의 염색체 안의 어떤 위치에서 나왔는지, 계통 염색체의 상대적인 위치는 어디인지(단일 가닥의 방향 정보는 암시적으로 저장되는데, 이는 전향과 후향 색인에 씨드가 저장되기 때문이다). 각 씨드에 대해서 이 값들을 정수 값으로 저장하는 것은 너무 많은 씨드의 숫자 때문에 수 많은 공간을 요구할 것이다.

대신에 각 씨드 발생에 대해서 색인 안에 하나의 단일 4바이트 정수가 저장되고, 보조 해시 테이블이 제시되었다. 4 바이트는 3 바이트의 블록 번호와 1 바이트의 오프셋으로 나눠진다. 참조 서열과 계통은 256 길이의 블록으로 쪼개진다. 참조 서열과 계통 서열들은 256 길이의 블록으로 쪼개진다. 블록 번호는 블록 테이블이라고 불리는 해시 테이블의 키로 이용된다. 오직  블록 테이블 엔트리가 있을 수 있기 때문에(키의 길이가 24 비트 혹은 3 바이트 이므로), 수많은 데이터가 각 테이블 엔트리 안에 많은 메모리 혹은 디스크 공간을 사용하지 않고 저장될 수 있다. 예를 들어, 발견된 씨드가 있는 염색체를 표현하기 위한 블록 테이블 안의 정수 하나를 저장하기 위해서는 오직 byte = 64 MB만을 필요로 한다. 애기장대에 대해 생성된 1억 천 9백만 씨드들의 모든 염색체 번호를 저장하기 위해서는 최대 450 MB가 필요할 것이다.

1 바이트의 오프셋은 블록 안에서 씨드의 시작 위치를 나타낸다. 그렇기 때문에, 블록들은 지놈 서열을 위한 저장소로 사용될 수 있고, 각 블록은  염기 쌍(1 바이트)을 저장할 수 있다.

그래서, 씨드의 위치는 색인으로 부터 블록 번호와 오프셋을 통해서 얻어질 수 있다. 블록 번호 덕분에, 계통, 염색체, 블록 서열의 시작 위치를 상수 시간에 얻을 수 있다. 씨드 시작 위치는 블록 테이블에서 블록 시작 위치들과 1 바이트 오프셋을 더함으로써 계산할 수 있다.

씨드 위치를 표현하는 4 바이트 조합된 숫자가 지놈에서 발견된 모든 씨드들에 대해서 하드디스크 뿐만 아니라 genomemapper 프로그램을 사용할 때 RAM에 저장되기 때문에, 가능한 작게 유지시키는 것이 매우 중요하다. 3+1 바이트 조합으로 인간 지놈을 12.5 Mio(메비, 옥텟의 단위로 8비트, 즉 1 바이트의 개수이다)에 저장할 수 있다. 블록 숫자 공간은 더 많은 디스크 공간과 메모리 사용으로 확장할 수 있다.


블록 서열 표현하기


다양한 계통 서열 정보를 GenomeMapper의 입력으로 하는 것은 참조 서열에 대해 계통 서열 정보를 전사 지놈 정렬 개념으로 이해할 수 있다. 이 정보는 계통에 대해 리드를 정렬하려는 것뿐만 아니라, 부수적으로 참조 서열에 대해서 정정렬할 때 사용될 수 있다. 참조 또는 계통에 있는 미스매칭 위치는 다르게 표현된다. 이것은 가령, 다른 계통들에서 발견되지 않은 계통 다형성으로부터 한 게통에서 발견된 참조에 대한 변이를 발견하게 해줌으로써 차후에 있을 컨센서스 분석을 쉽게 만들어 준다.

참조에 모든 계통 서열이 정렬되도록 하는 작업은 참조 서열 정보를 계통에 있는 서열의 표현에 합침으로써 구현되었다. 그런 까닭에, 계통 서열은 다음 방법으로 인코딩 되어야 한다.

8 비트 문자의 3 최하위 비트(LSB, least significant bit)는 계통에서 발견되는 한 염기를 인코딩하고, 4~6의 비트는 참조의 염기를 인코딩한다. ASCII 코드의 3 최하위 비트는 한 계통의 정보를 인코딩하기 위해서 설계되었다. 이것은 게통과 참조 서열에서 동일한 서열에 대해서는 인코딩이 필요 없음을 암시한다. 오직 다른 염기들, 즉 에를 들어, SNP, indel의 서열은 인코딩 되어야 한다. 그들의 최상위 2 비트들은 0이고, 이것은 이들을 일반적인 ASCII 코드 문자와 구분짓게 한다. 왜냐하면, 일반 문자의 경우 2 번째 최하위 비트(LSB)는 설정되어 있기 때문이다. 예를 들어, ASCII 코드 중 '0'과 같은 숫자 값을 표현하는 문자는 내부적으로 언제나 64보다 큰 정수 값을 가진다. 문자 중에 64 보다 작은 값을 가지는 문자들은 서열 추출 과정 동안에 디코딩 되어야 한다. 그림 2.2는 결과로 나타난 인코딩 테이블을 보인다.


ASCII code A/a: dec: 65/97, bin 01000001/01100001 -> A/a: 011

ASCII code A/a: dec: 65/97, bin 01000001/01100001 -> A/a: 011

ASCII code A/a: dec: 65/97, bin 01000001/01100001 -> A/a: 011

ASCII code A/a: dec: 65/97, bin 01000001/01100001 -> A/a: 011

ASCII code A/a: dec: 65/97, bin 01000001/01100001 -> A/a: 011

ASCII code -      : dec: 0        , bin 00000000                      -> -     : 000

(A)


 

A

C

G

T

N

-

 A

01 000 001
65, A

00 001 110

11, (AC)

00 001 111
15, (AG)

00 001 110
12, (AT)

00 001 110
14, (AN)

00 001 000

8, (A-) 

 C

00 011 001

25, (CA)

01 000 011
67. C

 00 011 111
31, (CG)

00 011 110
28, (CT)

00 011 110

30, (CN)

00 011 000
24, (C-)

 G

00 111 001
57, (GA)

00 111 011
59, (GC)

01 000 111
71, G

00 111 110
60, (GT)

00 111 110
62, (GN)

00 111 000

56, (G-) 

T

00 100 001
33, (TA) 

00 100 011
35, (TC)

00 100 111
39, (TG)

01 000 100
84, T 

00 100 110
38, (TN) 

00 100 000
32, (T-) 

 N

00 110 001

49, (NA) 

00 110 011

51, (NC) 

00 110 111
55, (NG) 

00 110 100

52, (NT) 

01 000 110
78, N 

00 110 000
48, (N-) 

 -

00 000 001
1, (-A) 

00 000 011
3, (-C) 

00 000 111
7, (-G) 

00 000 100
4, (-T) 

00 000 110
6, (-N) 


(B)

그림 2.2. 시퀀스 인코딩 (A) 코드는 상응하는 문자의 ASCII 코드로부터 유도되었다. (B) 정렬된 서열 정보를 위한 인코딩 표. 참조 서열의 염기는 행에, 계통 서열의 염기는 열에 나타냈다. 표는 이진 표현을 보이고, 각 정렬된 문자 쌍에 대해 그에 상응하는 정수와 문자열 표현을 보인다. *는 갭에서 갭으로의 정렬이 허락되지 않기 때문에 --와 같은 조합은 사용되지 않는다.


디코딩은 비트 연산으로 간단히 수행될 수 있다.(게통 서열을 위해서는 비트 마스크 00000111이, 참조 서열을 위해서는 비트마스크 000111000이 AND 연산과 사용될 수 있다). 사람이 읽을 수 있는 형태로 보이기 위해서 표에서는 괄호 안에 그것을 표시하고 있다. 오직 3비트가 각 계통 정보를 표현하기 위해서 사용되기 때문에 오직 8(=)개의 다른 가능한 조합을 가지며 이는 4개의 염기에 대해서는 충분하나 모든 IUPAC 기호를 표기하기에는 역부족이다. 그것이 A, C, G, T를 제외한 문자는 N으로 대치되는 이유이다. 이 변화는 서열의 표기를 변화시키지만, 모든 유일하지 않은 염기는 매핑 단계에서 미스매치로 인지되기 때문에, 이 변환이 정렬의 품질에 영향을 미치지는 않는다.

참조 또는 계통 서열에 정렬된 리드의 매핑된 영역을 얻기 위해서는 다른 연산을 할 필요가 없이, 서열을 다른 방식으로 파싱하기만 하면된다.


댓글