요즘 논문이나 기사를 보면서 Bioinformatics/Genomics 분야의 트렌드를 개인적으로 정리한 글로 본인의 무지에 의해 잘못된 정보를 전달할 수도 있음을 주의하시기 바란다. ㅋㅋㅋ
넘어가야할 허들 - 속도
지금까지는 속도 보다는 클라우드를 이용한 scale-up이 주된 테마였다면 이제는 기존의 파이프라인을 개선하는 방향의 speed-up이 주요 이슈가 되었다. 표준 파이프라인이라고 할 수 있는 BWA, GATK, SAMtools, Picard를 사용하는 경우 50X의 Human genome의 경우 variant call까지 16 코어 서버를 사용하는 경우 68-94시간이 소요된다. 물론 소프트웨어의 버전이나 병렬화를 어디까지 수행하는냐에 따라 이 시간은 달라질 수 있지만 말이다. 여기에 도전장을 내민것이 바로SpeedSeq으로 SAMBLASTER를 이용하여 메모리 상에서 바로 duplication을 마킹하고 Sambamba를 이용하여 멀티쓰레드를 이용한 sorting과 indexing을 수행함으로 기존의 Picard와 SAMtools 대비 비약적인 속도 향상을 노리고 있다.
다른 하나는 GATK를 병렬화를 통한 속도 향상을 꾀하는 ClusterK라는 스타트업이다. GATK SWE를 이용하여 2가지 부분에서의 GATK 파이프라인의 병렬화를 제공하고 있다. alignment 단계 (BWA-MEM)에서 500MB의 chunk들로 나누어 align을 수행(chunk당 15-20분 소요)한다. 각 chunk 단위의 align은 염색체별로 bam 파일로 combine을 수행한다.
염색체 단위의 alignment 파일은 30MB 단위로 split하여 각각 variant calling 작업을 수행한다. 이때 고려해야할 것이 바로 안전하게 30MB 단위로 split해야하는데 즉, indel등의 variant를 고려해서 안전하게 split을 수행해야 한다는 것이다. 잘못 잘랐다(아무런 옵션없이 even하게 자른다거나)가는 indel을 놓칠 수 있기 때문이다. Heng Li에 따르면 low complex regions (LCRs)를 고려하거나 pontial misassembly regions를 고려하여 safe split을 수행한다. 위에서 살펴본 2가지의 병렬화를 기반으로 GATK 파이프라인을 수행할 수 있도록 스케줄링을 하는 프레임워크를 제공한다면 GATK 파이프라인의 속도 향상을 노릴 수 있다.
일찍이 fastq 파일을 chunk로 분할하는 방법에 대해서는 많이들 사용하는데, 이게 BWA 알고리즘상 best hit 등에 문제가 발생할 수도 있으나 이따위쯤이야 무시해도 될만한 것이고 BAM 파일을 split하는 것 또한 기존에 염색체별로 나눈다 정도 였다면 safe split을 통해 안전하게 병렬화가 가능하다. 뭐 어쨌건간에 아직 human reference genome의 한계로 인해 이걸 우회할 수 위의 방법들이 유용하지만, 오늘 기사에 따르면 Seven Bridge Genomics에서 Graph Reference Genome과 Graph Aligner 등을 커머셜하게 만드는 펀드를 받았다고 하니 이것두 활용할 날이 얼마 남지 않았나 하는 생각도 들긴하다.
그림1. LCRs와 PMRs 사이트에 존재하는 soft clip reads
Graph Reference와 이를 지원하는 Aligner라고 하면 레퍼런스 게놈이 참 게놈스럽기 때문에 애초부터 염색체 1번부터 MT, X, Y까지... 1 bp 부터 주욱 일렬로 늘어선 레퍼런스는 우선 phasing 정보나 repeat 정보 등등을 표현하기는 태초부터 포기한지라 이를 그래프로 Ref를 만들고 이를 지원한다면야 clinical쪽으로 한발 앞서갈 수 있게 되는지라 이를 백만년 걸리지 않고 지금정도의 분석시간 수준으로만 해결한다면 굿 되겠다.
두 시퀀싱 센터의 걸림돌 제거 전략
위의 SpeedSeq이나 GATK SWE의 경우에는 기본적으로 BWA-MEM을 기반으로 하고 있다면 이번에는 아예 alignment 단계에서부터 속도 향상과 디스크 사용률을 낮추는 상용 소프웨어로 덴마크의 Genalice에서 만든 Genalice Map은 BAM 포맷대신 GAR 포맷을 사용하며, 처리속도나 용량이 GATK 파이프라인에 비해 월등히 우수하다. 이러한 속도와 중간/최종 결과물의 file size의 잇점으로 최근 마크로젠에서 해당 프로그램을 사용하기로 결정했다는 기사가 나오기도 했다.
그림 2. Genalice에서 제공하는 GATK 파이프라인과의 분석시간/File Size 비교
아마도 X10 데이터를 처리하기 위해서는 어쩔 수 없는 선택이었을 것 같기도 하다. 또 하나의 X10을 도입한 Garvan Institute는 AllSeq과 DNANexus 플랫폼을 도입해서 X10 데이터를 처리하는데, 공개된 HiSeq X Ten 데이터를 public하게 사용할 수 있다. HiSeq X Ten을 도입한 두 시퀀싱 센터의 경우 각각 클라우드 또는 상용 알고리즘으로 데이터 처리에 대한 걸림돌을 헤쳐나가는 모습을 볼 수 있다.
Genalice Map과 SpeedSeq 분석 좀 자세히 보기
코드설명: BWA 실행시 -R 옵션을 통해서 BAM 파일의 RG 헤드 필드와 태그를 추가하도록 한다. BWA를 수행하면 STDOUT으로 SAM 포맷이 출력되며 samblaster는 기본적으로 STDIN으로 SAM을 입력받아 STDOUT으로 Duplication이 마킹된(SAM FLAG 0x400) SAM 파일을 보낸다. 이제 SAMBAMBA를 이용하여 멀티쓰레드를 이용하여 sort를 수행한다.
자 그럼 split을 통한 병렬화와 결합한 완벽한 파이프라인은 좀만 기다리면 만들어 보여주겠다. 기다려라.