유전체 염기 서열 교정하기

1) NGS에서 유래한 짧은 read를 참조 서열 위에 정렬하는 것,
2) 이를 바탕으로 변이를 찾는 것(variant calling),
3) 그리고 최종적으로 변이 정보를 사용하여 참조 서열을 교정하는 것.

이 세가지는 모두 다른 활동이다. 만약 인간 유전체를 대상으로 연구를 한다면 3)은 그다지 접할 일이 없을 것이다. 인간 유전체 정보를 교정한다면 거창한 일을 한다면 모를까? 반면 미생물 유전체 연구를 하는 사람이라면 꽤 많이 접하게 되는 일이다.

가끔 실험 진화를 거친 미생물 균주로부터 어떤 염기서열이 바뀌었는지를 착기 위해 resequencing analysis를 할 때가 많다. 1) + 2)가 이때 필요한 작업이다.

예전에 일루미나를 이용해여 대충(?) 유전체 서열을 얻어 놓고 최근에 PacBio로 완성 수준의 유전체를 새로 얻게 되는 일이 많다. 물론 미생물의 경우이다. 결과만 놓고 보면 분명히 개선이지만 일루미나 데이터가 잘못하면 사장될 수도 있다는 것이 나의 고민이다. 일루미나의 sequencing coverage가 훨씬 높다는 것을 감안하면 쓸데가 전혀 없지는 않다. PacBio assembly에 일루미나 read를 얹어 놓고 variant를 찾으면 많게는 수십개의 염기서열 차이가 검출된다. 이를 3)번 과정을 이용하여 교정하면 된다.

이것 말고도 일루미나 데이터를 쓸 곳은 좀 더 있다. 너무 주절주절 여기에 기록하면 최근에 투고를 앞두고 있는 간단한 논문에서 다루는 내용을 미리 누설(?)하는 셈이 되니 일단은 자제력을 발휘하겠다.

나는 거의 전적으로 모든 NGS 관련 작업을 CLC Genomics Workbench에서 한다. 아직도 workflow를 다루는 방법을 제대로 배우지 못해서 마우스 클릭만으로 대량의 데이터를 하나씩 처리하려면 조금 번거로운 면이 없지는 않다. CLC에서는 2)의 기능이 최근 1-2년 사이에 많이 향상되었다.

그런데 내가 아는 바로는 3)에 해당하는 기능이 없다. 만약 국내에서 CLC Genomics Workbench를 보급하고 있는 (주)인실리코젠에서 내가 잘못 알고 있는 것이 있다면 알려 주시기를! 그동안은 consensus sequence를 추출하는 기능을 대충 사용하여 왔다. 엄밀히 말하면 다양한 방법으로 variant calling을 할 수 있지만(CLC에서는 version 8.0.x를 기준으로 quality-based variant detection과 fixed ploidy variant detection 및 low frequency variant detection 등 세가지 방법을 제공한다) consensus extraction은 특정 위치에 align된 read에서 가장 많이 나타는 염기를 기준으로 최종 서열을 뽑아주는 간단한 방법을 쓰게 된다.

가장 바람직하게는 CLC에서 찾은 variant를 VCF로 전환한 뒤 다른 프로그램을 써서 reference를 교정하면 될 것이다. 그런데 아직도 CLC에서 VCF 파일을 제대로 뽑는 방법을 찾지 못했다. 분명히 메뉴에 기능이 있는데, 이를 클릭하여 따라가다 보면 'reference track'을 설정하는 단계가 나오면서 여기에서 막히고 만다. 인간 유전체를 다루는 사람에게 친숙할 GATK 패키지를 사용하면 되겠지만 호미로 막을 것을 가래로 막는다는 기분이 든다.

그래서 어제는 아주 '무식한' 방법을 동원하고 말았다. CLC의 variant table을 텍스트 파일로 출력한 뒤, 간단한 Perl 스크립트를 작성하여 reference를 고치게 하는 것이었다. 사실 이러한 일을 이미 몇년 전에 한 적이 있었고 당시에 만든 스크립트를 다시 꺼내 쓰면 되는 간단한 일이었지만, 스크립트에 항상 충분한 설명을 달지 않는 나쁜 버릇 때문에 이를 들여다보느니 새로 짜는게 더 나은 상황이 늘 벌어진다. 아마 자기 혼자 쓸 목적으로(주로 연구자들) 스크립트를 작성하는 사람들은 전부 비슷할 것이다.

CLC가 만들어 주는 variant table이 보여주는 reference 기준의 위치는 substitution에 대해서는 매우 명백하지만 insertion이면 고민이 필요하다. Sample은 이 위치 다음에 삽입이 되어 있다는 것일까, 혹은 그 뒤에 삽입이 되어 있다는 것일까? VCF의 specification 정보를 찾아보았다. 그 이벤트가 일어난 바로 앞의 위치를 기록한다고 되어 있다.

이를 참조하여 reference를 고친 뒤 다시 read mapping을 하여 variant를 찾아 보았다. 더 이상 variant는 나타나지 않는다. 실수는 없었다.