Fork me on GitHub 단맛만좋아요 :: 단맛만좋아요 (6 Page)

Structural Variation
2015.03.19 07:28 | 바이오인포매틱스

NGS Short Reads를 이용하여 Strucural Variation을 찾는데에 있어서 depth of coverage (DOC), paried-end mapping (PEM, PE), split read (SR) 정보를 이용하게 된다. SV를 찾는데에 있어서 reference에 reads를 매핑하고 reference에 정확히 일치하는 않는 (not exact match to reference) read들을 SV를 찾는데 사용한다.


아래와 같이 60과 37이라는 부분은 reference에서 174 bases의 간격을 두고 있으며, 우리는 이 60과 37부분외에도 54 부분을 각각 서로 다른 말단에 가진 두개의 read를 가지고 있다고 하자. 

이 reads를 BWA를 이용하여 매핑한 경우 첫번째 read인 60, 37,54 에 대해서 reference에 60 부분만 매핑하고 나머지  37,54  부분은 soft clip으로 처리해 버린다. 두번째  read인 54,60,37 대해서도 중간의 60부분만을 매핑하고 54와 37부분에 대해서는 모두 soft clip으로 처리해 버린다. secondary alignments의 경우에는 앞서 soft clip을 처리하는 대신  hard clip으로 처리되는데 37(빨간색) 부분에 대해서 align되고 hard clip된 것을 확인할 수 있다.



같은 데이터에 대해서 GENALICE MAP은 좀 다른 전략을 취하는데 read에 deletion이 존재한다고 가정하고 read를 split하여 매핑한다면 60과 37이 아래와 같이 각각 reference에 매핑되게 된다. 본 예에서는 2개의 read만이 이러한 deletion을 지원하지만 이러한 read들이 충분히 많이 cluster를 이룰수 있을만큼 존재한다면 이부분은 deletion일 가능성이 크게 된다. 

위의 서로 다른 aligner를 보았지만, 이렇게 SV가 존재하는 부분에서는 clipping이 많이 관찰된다. 아래 그림과 같이 deletion 부분에 대해서 paired read 정보를 이용하는 경우 확연히 deletion을 확인할 수 있으며, 이러한 부분에는 soft clip이 다수 존재하는데 바로 soft clip이 이러한 deletion을 support하는 요인의 하나가 된다. 

현재까지의 SV 검출은 BWA를 이용한 후 그 뒤부터 SV 툴들이 진행하는데 aligner가 어느정도 받쳐준다면 훨씬 SV 검출에 있어서 수월하게 진행할 수 있게 된다. 또다른 aligner인 Spiral Genetics의 Anchored Assembly (이하 AA)를 보면 SV를 위한 노력의 흔적을 더 잘 볼 수 있다. AA는 4단계로 진행된다.


첫번째는 reference에 align전에 reference의 정보없이 read에 대해서 correction을 수행한다. 우선 read들에 대해서 k-merize를 수행하여 각 k-mer에 대한 read quality를 계산한다. 이때 low count k-mers는 에러로 간주하여 버려진다. 이것은 간단히 read trimming과도 유사한 것으로 대략 quality가 낮은 read의 뒷부분들의 일정 k-mer들이 N으로 표시되거나 하겠다. 다음의 두 단계는 reference와 match되는 read들은 제외하고, 나머지 read들만을 가지고 overlap graph를 생성한다. 이때에는 context보다는 kmer를 기반으로 de Brujin graph를 생성한다. 마지막으로 anchoring을 수행하는데 reference와 일치하는 그래프의 끝 단말을 각각 찾아내는데 이를 anchor라고 부르며 이는 SV의 breakpoint가 된다.



지금까지 aligner들이 SV를 찾는데에 어떠한 역활을 하는지에 대해서 간단히 알아보았다. 이제는 각 SV 툴들이 어떠한 전략으로 SV를 찾는지 다음에 알아보도록 하겠다.




저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2015.03.19 07:28
Currently 댓글이 없습니다. comments want to say something now?

요즘 논문이나 기사를 보면서 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 InstituteAllSeqDNANexus 플랫폼을 도입해서 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을 통한 병렬화와 결합한 완벽한 파이프라인은 좀만 기다리면 만들어 보여주겠다. 기다려라.

저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2015.03.16 13:39
Currently 댓글이 없습니다. comments want to say something now?

티스토리 툴바