Fork me on GitHub 단맛만좋아요 :: '바이오인포매틱스' 카테고리의 글 목록

바이오인포매틱스

  1. BGZF (Blocked GNU Zip Format) - published: 2015.07.01
  2. Split Reads - published: 2015.06.26
  3. Simple consensus approach improves somatic mutation prediction accuracy - published: 2015.04.27
  4. Structural Variation - published: 2015.03.19
  5. Bioinformatics (genomics) 트렌드 - 지금 필요한건 스피두 - published: 2015.03.16
  6. Detecting Somatic Mutations - Ensemble Approach - published: 2015.01.14
  7. Samtools를 이용한 genotype likelihoods 구하기 - published: 2014.10.31
  8. 당신이 개발자라면... - published: 2014.10.10
  9. WGS에서의 확장을 고려한 유전변이 검출 파이프라인 - published: 2014.09.30
  10. Somatic mutation calling in Low-allelic-fraction - published: 2014.09.04
  11. RNA-Seq Applications - published: 2014.04.17
  12. somatic mutation 찾기 - published: 2014.04.08
  13. 베이즈 정리를 정리하고 넘어가자 - published: 2014.03.26
  14. NGS 데이터에서의 Genotype and SNP calling - published: 2014.03.22
  15. GATK의 incremental joint discovery를 위한 reference model pipeline - published: 2014.03.11
  16. comparison of variant detection methods - published: 2014.03.10
BGZF (Blocked GNU Zip Format)
2015.07.01 16:18 | 바이오인포매틱스

Random Access

BAM 파일의 경우에는 BGZF를 이용하기 때문에 원하는 곳으로 빠르게 access가 가능하다. 우리가 흔히 사용하는 GZIP (GNU ZIP) 보다는 압축효율 (압축했을때 용량)이 떨어지지만 random access가 가능하다는 잇점으로 인해 BAM 파일(BAM의 경우 재빠르게 자신이 원하는 position을 뷰잉하는데 많이 사용하기 때문)에서 사용하는 기술이다.


용량이 큰 텍스트 파일을 압축해 놓고 파일의 어느 부분이던지 랜덤하게 액세스 가능하기 때문에 그 활용도가 높은데 특히나 클러스터를 이용하는 경우 파일을 분할하는 등의 I/O 작업이 필요 없기 때문에 그 활용이 매우 높다고 할 수 있다. 여러 활용중 하나로 FASTQ 파일에 적용하여 사용하고 있다. 

FASTQ 응용

일루미나는 Base Call을 수행한 후 기본적으로 GZIP을 이용하여 압축하게 되는데 이를 다음과 같은 명령으로 BGZF로 변환한다. 뭐 언뜻 보기에는 BGZF가 GZIP과 호환되기 때문에 이게 BGZF인지 확인해 봐야 한다. grabix 툴의 check를 이용하면 간단하게 알 수 있다.


# zcat TN1506L0327_2.fq.gz |bgzip /dev/stdin -c > TN1506L0327_2.fastq.gz

# /BIO/app/grabix/grabix check TN1506L0327_2.fastq.gz

yes

# /BIO/app/grabix/grabix index TN1506L0327_2.fastq.gz

# /BIO/app/grabix/grabix grab TN1506L0327_2.fastq.gz 1 10

@HWI-ST840:481:HVKTWADXX:1:1101:1213:1943 2:N:0:ATTGGCTC

GGTACCAGTACCATGCTGTTTTGGTTACTGTAGCCTTGTAGCACAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

+

BB@DFFFFFHHHHJJJJJIJJJJJIIJJJJIIJIJJJJJJJJIJJI#######################################################

@HWI-ST840:481:HVKTWADXX:1:1101:1986:1934 2:N:0:ATTGGCTC

ATATCTACACCACGCAAAGTGATGTGTAAGTGTGGGTGTTGCTCNCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

+

@@CDDFFAHFFHHIIIIIIHHIHECDHIIEFHHHGI@DDGIIIG#########################################################

@HWI-ST840:481:HVKTWADXX:1:1101:1874:1947 2:N:0:ATTGGCTC

CCTAGGAATCCAACTTACAAGGGATGCGAAGGACCTCTTCAAGAAGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

BWA 활용

이렇게 BGZF 포맷의 FASTQ 파일은 굳이 압축을 풀지 않아도 원하는 라인을 바로 확인이 가능하게 된다. 이렇게 되면 우리는 클러스터 노드에 다음과 같이 직접 파일을 분할 할 필요 없이 "grabix grab 파일명 start_line end_line"을 통해 다수의 노드에 BWA 매핑 작업을 분할하여 수행할 수 있다.


# bwa mem -M -t 7  -R '@RG\tID:1\tPL:illumina\tPU:1_2014-08-13_dream-syn3\tSM:syn3-normal' -v 1 GRCh37.fa <(grabix grab normal-1.fq.gz 20000001 40000000) <(grabix grab normal-2.fq.gz 20000001 40000000)


많은 생물학 데이터 FASTQ, FASTA, SAM, BED, XML, Tab 구분 파일 등등 모두 적용이 가능하며 다양하게 활용이 가능하다. 물론 여러개의 BGZF 파일을 압축해제 과정이 필요없이 하나의 BGZF로 묶어주는 방법도 장혜식 박사의 블로그에 있으니 자신의 데이터에 적용해서 활용해 보기 바란다.

Python을 이용한 BGZF 헤더 다루기

RFC1952 문서에 정의된 GZIP에서 extra field를 활용하는 BGZF는 extra field의 두개의 subfield id에 각각 66, 67의 값을 가지며, EOF 마커로 "1f 8b 08 04 00 00 00 00 00 ff 06 00 42 43 02 00 1b 00 03 00 00 00 00 00 00 00 00 00"의 16진수를 가진다.

 


Python을 이용하여 파일을 읽고 BGZF 포맷이 맞는지 확인하는 코드를 작성해 보도록 하자. 대학원 시절 프로토콜 정의하고 구현하던게 갑자기 생각이 난다. 암튼 뭐 어쨌건간에 우선 파일의 첫번째 블록의 헤더 부분을 읽어오는데 기본 헤더와 함께 SI1,SI2의 extra subfield의 일부분을 읽는다. 읽어오는 부분은 총 14바이트 하나씩 출력해보면 위에 정의된 값들이 출력되는 것을 확인할 수 있으며, XLEN 부분은 총 6바이트가 extra subfield로 지정된 것을 확인 할 수 있다.


from struct import unpack


bgzf_file=open('TN1506L0327_1.fastq.gz')

header = bgzf_file.read(14)


id1, id2, cm, flg, time, xfl, os,extra_len, extra_id1, extra_id2 = unpack('<BBBBLBBHBB', header)


>>> from struct import unpack

>>>

>>> bgzf_file=open('TN1506L0327_1.fastq.gz')

>>> header = bgzf_file.read(14)

>>>

>>> id1, id2, cm, flg, time, xfl, os,extra_len, extra_id1, extra_id2 = unpack('<BBBBLBBHBB', header)

>>> id1

31

>>> id2

139

>>> cm

8

>>> flg

4

>>> extra_len

6

>>> extra_id1

66

>>> extra_id2

67


BGZF 파일의 EOF를 읽어들여 파일이 제대로 생성된 것인지를 확인해 본다. 이전에 언급했듯이 파일의 마지막 28바이트를 읽어들여 EOF 마커와 동일한지 확인한다. 왜 BGZF의 첫블록을 읽고 마지막 EOF를 확인하는지는 많은 데이터가 짤리는 경우가 많아서 확인 과정을 거친 후 진행하는게 건강에 좋아서랄까? ㅎ


>>> bgzf_file.seek(-28,2)

>>> eof=bgzf_file.read(28)

>> eof

'\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00BC\x02\x00\x1b\x00\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00'


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

Split Read (SR)

Split Read(SR)는 하나의 read가 분리 (split)된 것으로, 여기서 분리는 read의 일부분이 reference에 align되고 나머지 일부분은 또 reference의 다른 부분에 align된 것으로 Chimeric Alignment라고도 한다. SR은 deletion, insertion, inversion, tandem duplication과 같은 structural variation을 찾는데 유용한 지표로 사용된다. 



Identification of a deletion in an individual genome by split read analysis


SAM에서의 SR 흔적

SAM파일에서는 SR을 표시하는데 SA 태그를 사용한다. SA 태그는 Chimeric Alignment를 표시하는데 사용되며 앞서 말했듯이 하나의 read에서 일부는 "A" 위치에 align되고 일부는 "B" 위치에 align되는 경우에 사용된다. 이때 CIGAR String에는 H (hard clipping)나 S (soft clipping)과 같은 operation이 같이 나타나게 된다. Read ‘r003’이 각각 +/- strand로 reference genome 의 서로 다른 두 위치에 align된 경우 SAM 파일상에는 이를 각각 다음과 같이 표시한다.

r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0; 

r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1; 


첫번째 라인에서 CIGAR string 5S6M은 5개 ‘gccta’는 5S, 나머지 ‘ AGCTAA’는 6M 으로 표시된다. 이때 S는 soft clipping (reference에 정확히 align되지 않은 부분으로 일반적인 D를 이용하여 deletion으로 표시하지 않는다. 왜나면 read의 일부가 또 다른 위치에 붙는 즉 clipping에 의해 split된 read이기 때문이다)이라고 하며 soft clipping된 sequence는 SEQ 상에 표시된다. 두번째 라인을 보면 똑같은 read이지만 ref genome의 다른 부분에 6H5M으로 align 된것을 확인할 수 있다. 6H는 6개의 sequence가 hard clipping 되었다는 것을 의미하며, 이미 같은 read가 한번 soft clipping되었기 때문에 그다음의 read는 hard clipping으로 처리되고 이때 SEQ 상에 clipping된 부분을 표시하지 않는다. 

Non-overlapping base pair

Structural variation 분석에 있어서 이런 SR는 단순히 SA 태그를 지닌 read를 모두 사용하는 것이 아니다. 말그대로 read가 large deletion에 의해 둘로 split된 경우 split된 read는 서로 overlap이 되지 않고 정확히 sequence가 겹치지 않고 둘로 나뉘어져야 한다. 그러나 뭐 정확히 모든 read들이 그렇지 않을거고 암튼 그래서 SR만을 꺼내서 사용할 경우 SA 태그에 명시된 cigar와 서로 비교하여 두 chmeric read간에 overlap되 않는 non overlap이 일정 length 이상인 것만 추출하며 이때 기본적으로 20을 할당한다.


아래 그림은 36S54M (-)과 38M52S(-)로 서로 overlap되는 길이는 3bp이며, 첫번째 read 입장에서는 non-overlap length는 52bp, 두번째 read 입장에서는 36bp가 된다. 


Max Split alignment count

위에서 살펴본 non-overlap length에 이어 SR을 찾는 옵션 중 하나는 하나의 read에서 얼마나 많은 split alignment (split-read mapping)을 허용할 것인가에 대한 것이다. 하나의 read가 split되어 여기저기 refrence에 덕지덕지 붙는다는 것은 말이 안되고 기본적으로 2군데 align되어야 이걸 SV의 흔적으로 생각할 수 있다. 따라서 기본값은 2로 지정된다. 최대 split 갯수를 확인하는 방법은 SA 태그에 split된 위치가 ';'으로 분리되어 보여지는데 이것을 가지고 read가 몇개로 split되었는지 알 수 있다.

if "SA:" in samTages:

   if(len(samTages.split(";"))) <= 2:

      print "Yes, split read"

Strand normalize for non-overlap base pair

서로 split read간의 overlap되는 base를 계산하기 위해서는 strand normalized 과정이 필요한데 strand normalized를 이해하기 위해서는 우선 paired-end 시퀀싱에서의 mapping에 대해서 알고 넘어가야 한다. 아래 그림처럼 SAM에서는 Reference 방향으로 모두 통일해서 SEQ를 표시하고 이것이 원래 FASTQ에서 reverse된 것인지는 flag에 0x10을 넣어 주어 이게 reverse된 것이라는 것을 알려준다.


Spirt Read가 SAM 파일에서 Read는 F(orward) strand로 매핑되고 SR은 SA 태그에 R(everse) strand ('-'라고 SA 태그에 표시됨)라고 표시되었다면 서로 겹치지 않은 부분(non-overlapping base pair)을 계산하기 위해서는 strand를 맞추어 준 후에 CIGAR string을 가지고 strand에 맞추어서 계산해야 한다. 

Python을 이용한 구현

지금까지 우리는 Split Read라는 것의 정체를 알기 위해서 Soft/Hard Clipping, FASTQ 파일 포맷, F(orward)/R(everse) 등의 개념에 대해서 알아보았다. 이제는 우리는 이걸 코드로 표현해야 한다. 어려워말고 차근차근 따라와 주기 바란다. 여기서 사용할 코드는 python으로 진행하도록 한다.


우리의 구현할 프로그램의 목표는 SAM 파일을 입력받아 우리가 원하는 특정 Split Read만을 다시 SAM 포맷으로 출력하는 프로그램이다. 여기서 우리가 원하는것은 duplication reads는 뺀 나머지 read들 중 SA 태그가 존재하는 즉, split read에서 최대 split-read mapping이 2개이면서 서로 겹치지 않는 read가 20bp 이상인 것들만 출력한다.

Skeleton code

프로그램은 일반적으로 명령어를 실행하면 해당 명령의 help 페이지를 출력하고 이때 프로그램에서 사용할 수 있는 옵션을 보여주게 된다. 또한 명령어와 함께 사용되는 옵션에 대해서 처리할 수 있어야 한다.

#!/user/bin/env python

import ****


def 실제_수행할_메소드:

   #split read 추출


class 데이터를_저장할_구조정의:

   #실제_수행할_메소드에서 사용할 자료 구조 정의


def 부가적으로_필요한_메소드:

   #실제_수행할_메소드와 별도의 부가적인 메소드 정의


def main():

   #프로그램의 help 정보 표시

   #프로그램의 옵션들(인자들) 처리

   #실제 수행할 메소드 호출

    실제_수행할_메소드(인자들) 


if __name__ == "__main__":

   sys.exit(main())

Command line options

처음으로 우리는 프로그램의 옵션과 도움말을 처리할 것인데, 이건 main() 메소드에서 처리한다. optparse의 OptionParser를 이용한다. 우선 usage에 프로그램의 설명, 간단한 사용법을 넣고 paser를 이용하여 옵션을 처리한 후 꼭 필요한 옵션이 누락된 경우 도움말을 출력하여 옵션을 설정하도록 권유한다. 필수 옵션이 있는 경우에는 실제_수행할_메소드를 호출하는 것으로 마무리한다.


아무런 옵션 없이 프로그램을 실행하면 다음과 같이 위에서 지정한 도움말과 프로램의 옵션이 출력된다. 여기서는 -i 를 필수 옵션으로 지정해서 해당 옵션 없이 프로그램이 실행되면 도움말을 출력하며, 옵션이 설정되어 있다면 메인 메소드인 extractSplitRead를 실행한다.


SAM 포맷 저장을 위한 SAM 클래스 정의

메인 메소드인 extractSplitRead를 정의 하기 전에 SAM 파일을 한줄한줄 읽어들여와 각 read의 정보를 저장할 SAM 클래스를 정의해 보자. 물론 한줄을 읽어 들여 split()을 이용하여 리스트에 넣어 사용해도 되지만 말이다. SAM 클래스는 SAM 파일을 한 줄 읽어 이를 탭으로 분리한 리스트 형태를 입력으로 받게 되며 SAM 포맷에 맞게 query, flag, ref, pos, mapq 등을 설정한다.

메인 메소드 작성

우리가 입력으로 받은 SAM 파일을 담을 SAM 클래스도 정의했다면, 이제 extractSplitRead라는 이프로그램의 메인이 되는 메소드를 정의하도록 하자. 이 메소드는 입력 파일, split 갯수, duplication 허락여부, 최소 오버랩길이를 인자로 받아 제일 처음 입력 파일을 한줄 읽어 다시 탭으로 구분한 리스트에 담아 이를 SAM 클래스에 넘겨준다.


이렇게 읽은 한줄이 duplicate read이지를 확인하고, SA 태그도 있는지를 확인하여 SA 태그가 존재한다면, split = 1 이라고 지정하여 split read임을 표시한다. duplication을 허용하지 않고 (includeDups==0), flag 값을 통해 "read is PCR or optical duplicate (1024)"가 표시된 read라면 해당 read는 skip (continue)한다. split read는 다시 SA 태그 부분을 좀 더 파싱하여 split read의 CIGAR string을 mateCigar 변수에 저장하고 strand 정보를 mateFlag 에 저장('-' strand는 16)한다. 

Overlap 및 Non-overlap 길이 구하기

우리는 지금까지 SAM 파일을 읽어 해당 read가 split 된 read인지를 확인하고 해당 read를 SAM 클래스의 인스턴스로 만들고 split된 read의 CIGAR string과 strand 정보를 따로 보관하고 있다. 이제 마지막으로 하나의  read지만 split된 read가 서로 얼마나 overlap 되는지를 파악하는 것만 남았다. overlap되는 부분을 계산하기 위해 2개의 메소드 extractCigarOps, calcQueryPosFromCigar 메소드를 우선 정의한다.



extractCigarOps는 CIGAR string을 파싱하여 해당 read가 reverse strand (0x0010)인 경우 CIGAR string을 뒤짚어서 forward 형태로 만든 후 cigarOp라는 클래스의 인스턴스로 만든다. cigarOp는 CIGAR의 length와  operation을 가지게 된다. (ex, CIGAR string이 31S59M인 reverse strand의 경우 59M31S로 만들어 버린다) 이렇게 split된 read에 대해서 각각 readCigarOps와 mateCigarOps에 forward strand 형태의 CIGAR string의 정보가 담긴 cigarOp 인스턴스를 생성한다.


calcQueryPosFromCigar은 위에서 strand에 맞추어진 CIGAR를 입력 받아 각 read의 mapping된 start, end position과 전체 read의 length를 계산하여 queryPos의 인스턴스로 생성한다. 



이제 위의 queryPos의 인스턴스를 가지고 두 read가 서로 overlap되는 곳의 length와 non-overlap되는 곳의 length를 구하고 non-overlap이 기본값인 20 이하인 경우만 해당 read를 split read로 인식한다.

결론

지금까지 split read가 NGS에서 어떠한 의미이고 어떻게 활용될 수 있는지 그리고 실제 프로그래밍으로 어떻게 구현하는지까지 모두 알아보았다. 기나긴 글을 끝까지 읽어 주어서 고맙고 항상 생각했던 것을 코드로 직접 옮기면서 문제를 해결해본다면 더 많은 것들을 얻을 수 있을 것이다.

참고문헌

위의 코드는 LUMP 소스 코드를 기반으로 작성되었습니다.

LUMPY: a Probabilistic Framework for Structural Variant Discovery. 

Delly: structural variant discovery by integrated paired-end and split-read analysis

Genome structural variation discovery and genotyping (2011)

Computational methods for discovering structural variation with next-generation seuqncing (2009)

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

지난 동계유전체학회 워크샵의 암유전체 분석을 위한 Somcatic Mutation Calling에 관한 내용입니다. "A simple consensus approach improves somatic mutation prediction accuray"라는 논문을 바탕으로 SomaticSniper, VarScan2, MuTect을 이용하여 각각 somatic call을 수행한 후 각 툴에 대해서 filter를 적용하고 consensus 데이터셋을 만들어 이를 실제 validation하는 과정에 대한 내용입니다. 물론 데이터는 TCGA Benchmark 데이터셋을 이용했습니다. 마지막, validation 부분을 업데이트할 부분이 좀 있는데 우선 공유합니다.





저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2015.04.27 16:27
Currently 댓글이 없습니다. comments want to say something now?
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?
Detecting Somatic Mutations - Ensemble Approach
2015.01.14 13:50 | 바이오인포매틱스

예전에도 두어번 블로그글을 통해 variant calling에 있어서 Ensemble approach에 대해서 언급했었더랬습니다. comparison of variant detection methods, somatic caller는 뭘 사용해야 하나요? 를 참고하세요. 오늘은 그 끝판왕으로 준비했습니다. 왜냐구요? 잠시 광고 하나 하고 넘어가려고 합니다. 한국유전체학회 동계 심포지엄이 2월 4일부터 진행되는데 올해는 이틀간에 걸쳐 "Somatic Calling 알고리즘 소개 및 실습" 워크샵이 준비되어 있습니다. 현재 저조한 등록을 보이고 있다고 합니다. 오늘은 그래서 워크샵에서 진행할 내용을 미리 소개하는 시간을 가져 보려고 끝판왕으로 준비했습니다. ;-) 더 안오실지도 모르겠군요.


아! 그리고 저번 포스팅에서 구글애드센스 광고를 달았다고 말씀드렸는데요. 구글에서 광고수익을 지불 방법을 확인하는 우편도 도착했구요. 현재 무려 $15의 수익을 올리고 있습니다. 이자리를 빌어 이글을 보시는 모든분들께 감사하다는 말씀을 전합니다. ㅋㅋㅋ


Somatic Mutation Calling Tools

그럼 본론으로 들어가 보도록 하겠습니다. 뭐 요즘 핫한 분야는 바로 somatic mutation 찾기 아닐까 합니다. 암 샘플의 heterogeneous하고 normal cell과의 dilute 등등 그만큼 챌린저블한 분야이며, 무엇보다도 실제 clinical에 응용될 수 있는 부분이라 그런것 아닐까 생각이 듭니다.Ding 아주머니 (워싱턴 대학에서 SomaticSniper와 VarScan을 만드신)의 최근 리뷰 논문에 따르면 아래와 같이 수많은 cancer genome 분석 툴들이 존재합니다. (중간에 잘린 표입니다.)



이미 잘알려진 JointSNVMix, MuTect, SomaticSniper, Strelka, VasrScan 등이 우선 눈에 들어오고 최근 Dream Challenge의 Somatic Mutation Calling Challenge주인공들이 만든 툴들이 보이네요. 


 ICGC-TCGA DREAM Somatic Mutation Calling Challenge의 주인공들

역시 Broad-너네가 다해 먹어랏, Ding 아주머니도 보임


각각의 툴들에 대한 알고리즘이나 특징은 해당 논문을 찾아보시거나 여러툴들을 비교한 논문들이 꽤 존재합니다. 따라서 툴 소개는 이것으로 마무리하려고 했지만... Wang의 "Detecting somatic point mutation in cancer genome sequencing data: a comparison of mutation caller"를 잠깐 보면 다음과 같이 각 툴들을 소개하고 있습니다. 끝



하나로는 부족해 우리 합치자!

비단 bioinformatics 분야뿐만 아니라 다른 곳에서도 같은 목적을 위해 여러 툴들이 합쳐 각각의 단점은 상쇄시키고 장점을 부각하는 일명 ensemble 방법 (또는 consensus, combining)이 많이 사용되고 있다. 그런데 문제는 이 ensemble이라는 방법이 걍 사용한다고 해서 좋은 것 아니라는 것입니다. 좀 세련되게 각 툴의 결과와 feature들을 뽑아서 training을 시키거나 암튼 consensus를 구하는데에도 갖가지 어려운 수식이 오가면서 ensemble을 적용하는 것 또한 만만치 않다는 것입니다. 그래서 기름기를 쫙 뺀 바로 저 멀리 남반구 호주 멜버른의 Goode가 A simple consensus approach improves somatic mutation prediction accuracy"라는 제목부터 simple을 달고 나온 방법을 소개해보려고 합니다.


Data Sets

27 ovrian tumor와 그와 매칭하는 germline 샘플을 HiSeq 2000 으로 whole exome sequencing을 평균 102~225x로 100bp paired-end로 뽑아냅니다. 그리고 나서 JointSNVMix2, MuTect, SomaticSniper를 돌리죠. 왜 저 3개의 툴을 사용했냐? 메이저 tumor 논문에서 사용한 툴들이니까요 :-) 그리고 나서 총 9,226개의 somatic SNV를 prediction합니다. 평균 샘플당 321개이며, 툴별로 보면 S와 J는 샘플당 170개 정도 M의 경우 좀 더 보수적으로 115개의 SNV를 찾아내게 됩니다.



중요 포인트는 Non-reference allele frequency

그런데 somatic이라고 뽑은 mutation들의 해당 read가 reference와 같은 allele인지 아닌지를 한번 더 살펴보게 되며 germline에서 보면 J와 S에서 높게 나타나는 행패를 부립니다. 즉 somatic이라고 뽑은 놈이 germline mutation이게 되는 것입니다. 즉 이 놈들은 False Positive인 놈들인거죠. (그림의 빨간색 박스 참고) 반면 M은 somatic을 잘 뽑아낸것 처럼 보입니다. 여기서 그럼 somatic이라고 뽑은 놈들이 germline mutation 즉 germline에서 non-reference인 것들을 filterout 시킨다면 좀 더 True Positive에 가깝게 되겠지요. 이건 다음 filtering 부분에서 다시 보겠습니다. 




이제 tumor에서 non-refernece allele frequency를 보면 M,J,S와 JS, MJ의 경우 낮게 보이는데요. 요것 또한 somatic이라고 뽑은 놈들이라면 tumor에서 높은 non-reference  allle를 보여야 한다는 거죠.


즉, 3개의 툴이 공통으로 찾은 somatic mutation에 대해서 위에서 언급한 non-reference allele에 대한 필터를 거친다면 더욱 True Positive에 가까운 mutation을 찾을 수 있게 되겠죠. 그래서 이들은 찾은 mutation에 대해 sanger sequencing을 수행해서 validation set을 가지고 consensus call과 몇몇 filtering 과정을 거쳐 성능이 향상되는지를 확인하게 됩니다.


Additional filtering 왜? consensus의 specificity를 높이려구

앞서 언급했듯이 non-refence allele 부분과 관련한 부분을 해결코자 tumor/normal 샘플(.bam)에 대해서 GATK Unified Genotyper를 이용하여 tumor와 germline에 대한 SNV를 calling 합니다. 이걸가지고 consensus call에 대해서 우선 GATK UG의 germline SNV를 제거합니다. 그리고 나선 GATK UG가 tumor에서 찾은 SNV와 겹치는 것만 남기게 되죠. 이렇게 하면 아래 처럼 183개중 50개만 맞췄던 결과가 113개중 48를 맞추게 됩니다. 그에 더불어 mate-rescued된 read를 통해 필터를 더 추가시키게 되면 결과는 87개중에 48를 맞추게 됩니다. 결국 2개의 true positive를 잃게 되지만, 183개에서 87개로 확 줄어둔 즉 specificity를 확 증가시키게 되는 결과를 보입니다.



결론

결론적으로 간단한 consensus와 필터링만으로 validation rate를 높이고 sensitivity를 증대 시킬 수 있게 되었다는 것입니다. 하지만 그에 따라 다양한 툴들을 돌려야 한다는 즉 computation 리소스는 부담으로 남게 되며 multiple 툴을 병렬로 빠르게 돌리는 것에 대한 논의가 또 있어야 하겠죠.


뭐 어찌되었던간에 광고로 돌아가서 2월달 워크샵에서는 TCGA 데이터를 가지고 위의 consensus 메소드를 실제로 구현하고 검증하는 과정과 더불어 ranking을 통해 좀 더 연구자에게 보기 좋은 결과물을 얻는것이 실습으로 진행되니 관심있는 분들의 많은 참여를 바랍니다. 


안녕, ;-)

저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2015.01.14 13:50
Currently 댓글이 하나 달렸습니다 comments want to say something now?
Samtools를 이용한 genotype likelihoods 구하기
2014.10.31 13:32 | 바이오인포매틱스
samtools가 버전 0.1.19를 마지막으로 major 번호가 올라갔습니다. 바로 1.0 버전대가 탄생한것이죠. 홈페이지도 이제는 www.htslib.org를 사용합니다. 흔히 우리가 말하는 Samtools는 Samtools, BCFtools, HTSlib 3개로 구성되어 있습니다.  

Samtools는 SAM/BAM/CRAM 포맷의 파일을 읽고/쓰고/편집하고/인덱싱하고/볼 수 있는 툴입니다. BCFtools는 BCF2/VCF/gVCF 포맷의 파일을 읽고/쓸 수 있으며 SNP나 short indel의 sequence variants를 calling/filtering/summarising 할 수 있는 툴입니다. 그리고 그 기반은 high-throughput sequencing 데이터를 다루는 바로  C로 작성된 HTSlib이라는 라이브러리를 기반으로 하고 있습니다. Samtools를 이용한 Variant Calling에서도 다음과 같은 명령으로 변경되었습니다. 
samtools mpileup -ugf ref.fa sample1.bam sample2.bam sample3.bam | bcftools call -vmO z -o study.vcf.gz
잠깐 살펴보면 예전에는 bcftools의 view가 call이라는 명령어로 변경된 것을 확인 할 수 있을 겁니다. 그럼 간단하게 Samtools를 이용한 variants를 identifying하는 방법에 대해서 알아보도록 하겠습니다. mpipleup 명령은 aligned read 즉 BAM 포맷의 파일로 부터 genotype likelihood를 계산하는 명령입니다. 해당 샘플의 모든 site에 대해서 genotype likelihood를 계산하게 됩니다.

The Variant Call Format (VCF) Version 4.1 Specification을 보면 INFO 필드의 몇몇 키워드에 대한 설명이 나오는데요. GL 필드와 이와 관련한 PL, GP 필드에 대해서 우선 알아보겠습니다. GL 필드는 앞선 genotype likelihood에 대한 정보를 담고 있는 필드 (Floats)로 genotype likelihoods comprised of comma separated floating point log10-scaled likelihood for all possible genotypes given the set of alleles defined in the REF and ALT fields.

예를 들어 아래와 같이 reference allele가 T라고 할때 alternative allele가 C와 G 두 개가 있다고 합시다 (T=0, C=1, G=2) . 그리고 PL 필드를 보면 콤마(,)로 구분된 총 6개의 값이 보입니다. diploid인 경우 AA, AB, BB, AC, BC, CC 6개의 genotype에 대해서 likelihood 값을 계산하게 됩니다. REF는 0 , ALT는 1,2로 표현되는 경우 (0,0), (0,1), (1,1), (0,2), (1,2), (2,2)로 표시할 수 있겠습니다. 그럼 여기서 TG genotype은 몇 번째에 있는 값인가요?라고 물어보면 F(j,k) = (k*(k+1)/2)+j 즉, F(0, 2) = 3 이므로 'TG'는 genotype likelihood는 4번째에 존재하는 99가 genotype likelihood 값이 됩니다.
chr1 12887054 . T C,G 194 . DP=253;VDB=0.423501;SGB=-0.693147;RPB=2.7979e-06;MQB=2.50009e-15;MQSB=3.79278e-05;BQB=0.860152;MQ0F=0.41502;AC=1,1;AN=2;DP4=60,0,147,20;MQ=15 GT:PL 1/2: 255,100,99,213,0,231
 Genotype Phread Scaled Score   Likelihood p-value
AA, (0,0), TT 255
AB, (0,1), TC 100
BB, (1,1), CC 99
AC, (0,2), TG 213  
BC, (1,2), CG 0  
CC, (2,2), GG  231  




glf 포맷,,, glfsingle, samtools-hybrid,,, 등등등 추가해야함,,, 우선 skip 추후에 지속 업데이트하기로 한다. 

아이고 의미없다.  
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.10.31 13:32
Currently 댓글이 없습니다. comments want to say something now?
당신이 개발자라면...
2014.10.10 15:07 | 바이오인포매틱스
유전체데이터를 다루는 툴이나 스크립트를 만들고자 한다면, 다음의 프로그램, 라이브러리를 눈여겨 보고 응용할 것. 재미있는것은 이제 클라우드상의 데이터도 htsjdk 라이브러리에서 직접 핸들링이 가능해진다는 것이다.

gatk-tools-java

Tools for using picard and gatk with genomics API. getting reads from GA4GH genomics api and exposing them as SAMRecord "Iterable" resource. These will be used for subsequent work on enabling HTSJDK to use API data as input.
GA4GHPicardRunner wrapper around picard tools that allows for INPUTS inot Picard tools to be ga4gh:// usls.
https://github.com/googlegenomics/gatk-tools-java

htsjdk

java framework used by picard and GATK to access reads in bam files and variants in VCF files
https://github.com/samtools/htsjdk

picard

Java tool set using htsjdk to perform tasks on fastq, sam, bam and vcf files. Loop-based code, son no (easy) multiprocessor support.
https://github.com/broadinstitute/picard

GATK-framework

Java-framework future extending htsjdk with capabilities such as accessing pipleups. Defines a set of walkers for developers to extend to traverse bam or vcf files.
https://github.com/broadgsa/gatk

GATK-tools

A set of tools using the GATK-framework for analysis of bam and vcf files.
https://github.com/broadgsa/gatk-protected

htslib

C-framework for manipulation of BAM, CRAM, and VCF files.
https://github.com/broadinstitute/htslib

samtools

Tool set built using htslib for manipulation of sam files.
https://github.com/samtools/samtools.git

gamgee

C-framework further extending htslib. 
https://github.com/broadinstitute/gamgee


Source from: 
http://gatkforums.broadinstitute.org/discussion/4510/writing-custom-tools-htsjdk-vs-picard-vs-gatk
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.10.10 15:07
Currently 댓글이 없습니다. comments want to say something now?

WGS 분석에서의 bottleneck

alignment와 variant calling 단계에서는 cpu/mem 부분에서 bottleneck이지만, align post-process (base quality score recalibration, realignment around indels)과 variant post-process 단계에서는 disk의 io가 bottleneck으로 작용한다. 

대량 샘플에서의 병렬 네트워크 파일 시스템

단일 샘플 분석과 같은 경우 NFS가 유리하지만, 대량의 샘플을 분석하는 경우 Lustre나 GlusterFS와 같은 병렬 네트워크 파일 시스템이 유리하다. 당연히 io가 분산되기 때문에 적은 수의 샘플에서는 병렬 네트워크 파일 시스템이 불리하게 작용하지만, 다수의 대량 샘플 분석시에는 오히려 부하가 분산되면서 안정적이며 빠른 속도로 분석이 가능하다.

Post-alignment는 필요한가?

간단히 deduplication만 수행한 minimal BAM preparation과 앞서 이야기한 align post-process 단계를 거친것과 결과는 어떨까? 우리는 그동안 Broad의 GATK best-practice에서 이부분을 수행하는 것이 결과에 좋은 영향을 준다고 하여 무조건적으로 수행했다. 엄청난 IO와 분석 시간의 증가라는 cost적인 측면을 감내하면서 말이다. 일례로 30x의 whole genome sample의 경우 16 cores를 1샘플에 할당하여 분석할경우 약 12~16 시간이 이 부분에서 소요된다.

결과는 놀랍게도 deduplication만 수행한 BAM의 경우가 오히려 좋은 성능을 보인다는 것이다. recalibration과 realignment 과정없이 HaplotypeCaller와 FreeBayes에서 오히려 좋은 성능을 보인다는 것이다.

요즘 GATK 진영에 대해서 두가지로 요약된다. genome 데이터에 대해서 너무 많은 가공(?)을 한다는 것이다. 가령 위에서 언급한 alignment post-processing 등과 같이 말이다. 또다른 하나는 GATK의 핵심 멤버의 이탈(?)이다. 

모든 분야가 마찬가지겠지만, 영원한 1등도 없는 것이요. 1등은 그자리를 놓치지 않기 위해 비밀스런 고급 기술들을 쳐 넣기 때문에 오히려 역효과가 날 수도 있다는 것이다. 헨리가 말했지 않았던가 난 GATK 세부 알고리즘 몰라요. 그냥 samtool로 원하는것 하세요. 한번 곰곰히 생각해봐야할 부분인듯하다. 너무 GATK의 노예로 살아온건 아닌지 말이다. 

Variant calling pipeline의 성능 올리기

그럼, 이제 생각해 볼 부분은 low complexity regions, high depth filters, VQSR을 어떻게 variant calling에 적용할까?라는 부분이다. 최근 두개의 논문을 우선 살펴보면 어떻게 variant call을 평가하고 filtering 할지에 대한 도움을 받을 수 있다.

헨리의 "Towards Better Understanding of Artifacts in Variant Calling from High-Coverage Samples"와 Michael의 "Analytical validation of whole exome and whole genome sequencing for clinical applications"이다.


 
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.09.30 17:11
Currently 댓글이 하나 달렸습니다 comments want to say something now?
Somatic mutation calling in Low-allelic-fraction
2014.09.04 18:27 | 바이오인포매틱스
뭐 어쩌다보니 cancer까지 흘러 들어와 버렸다. cancer 분석에 대한 개념 정리는 뒤로 미루고 분석툴에 대한 내용으로 시작한다. 뭐눈에는 뭐만 보인다더니 어째 cancer 분석 논문보다 분석 툴(알고리즘)에 대한 논문만 넘쳐 난다는 생각이 든다. 

Somatic mutation detection

뭐 이것저것 많지만 우선 mutation dection 그것도 SNV만을 가지고 시작하자. 이 somatic mutation calling이라는게 germline mutation calling보다 복잡하다. 왜 그런지는 시간나면 지면을 할애해 설명하고 여기서는  cancer genomics의 somatic mutation의 SNV에 한정한다.

수많은 cancer genomics 관련 툴들이 나와 있지만, 필자 맘대로 selection한다. 뽑힌 툴들은 앞으로 나오기 바란다. Broad의 Mutect (Nature Biotechnology 출신), 워싱턴 대학의 VarScan2(Renome Research 출신), 연세대의 Virmid (Genome Biology 출신), 일루미나의 Strelka (Bioinformatics 출신) 4명 앞으로 나와! 

Low-allelic-fraction에 강한 엄친아 Mutect

베이지안 기반의 low-allelic-fraction (낮은 빈도의 체세포 변이를 발견)에서 강력한 힘을 보여주며, Broad 출신답게 java와 베이지안을 위해 dbSNP나 cosmic 파일을 입력으로 받는다. 오죽하면 논문 제목도 "Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples" 이다. 

그럼 진짜 sensitive한지 중립적인 상황에서 Mutect를 평가한 논문을 보도록하자. BMC Genomics에서 나온 "Comparison of somatic mutation calling methods in amplicon and whole exome sequence data" 를 보면 NIST-GIAB의 NA12878 데이터에 1,000 Genomes Project의 NA19129 샘플을 섞어 0%, 8%, 16%, 36%, 100%의 암샘플과 같은 역활을 하도록 데이터를 각각 만들어 평가를 진행했다. 아래 그림과 같이 C04 (8%) 즉, low-allelic-fraction에서  좋은 sensitivity를 보인다. (물론 Strelka가 amplicon의 경우 높지만)  - 원래 위의 논문은 추후 clinical에서의 응용을 고려해서 Amplicon과  Exome에 대해서 각각의 툴들이 어떠한지를 보여주려는 논문인데.. 암튼..


MuTect는 specificity에서도 가장 낮은 FPR(False Positive Rate)를 보인다. 


깔끔하게 jar 파일 하나로 실행이 가능하다는 장점과 친절하게도 상용 목적을 위해 별도의 라이센스와 지원까지 갖춘 놈이다. 역시 엄친아다. 별도의 매니저도 있고... 아래 그림처럼 결과  VCF 파일도 Normal, Tumor 두개의 정보가 보여지며 실행시 각각의 이름도 지정해 줄 수 있어 visualization할때 편리하다. 다만 thread 지원이 안된다는 것과 베이지안 계산 등등으로 Exome 1 pair의 경우 약 10시간정도의 running time을 가진다.


분석 속도 개선 여지는 있지만 MuTect과 견줄만한 Virmid

앞서 언급한 평가 논문에는 등장하지 않지만 Virmid 또한 Java로 구현된 caller로 아래와 같이 80%의 Concordnat Rate를 보이며 MuTect과 비슷한 성능을 보인다. 물론 분석 속도와 결과 VCF에 대한 약간의 개선이 필요하다. 아래 데이터는 위암 엑솜 샘플을 각각 MuTect와 Virmid를 이용하여 기본 옵션으로 분석한 결과이다.

MuTect과 Virmid

암샘플 분석에 있어서 Ti/Tv ratio가 큰 의미는 없겠지만, MuTect과 Virmid 각각 1.93과 1.95를 보이며, Known Site의 Novel site의 경우 1.84와 1.83의 비슷한 ratio를 보인다. 뭐 지금 상황에서 두개의 툴을 비교하기에는 어려움이 따르기 -.-;; 때문에 서로간의 비교는 그냥 calling된 숫자에 기반한다. 아울러 Virmid의 논문을 봐도 두개의 툴이 서로 비슷한 성능을 나타내는 것을 확인할 수 있다. 따라서 이번에 사용한 암샘플(위암) 분석이 비교적 잘?되었다고 하자.



위 그림은 Virmid 논문에 나온 그림으로 Virmid와 MuTect를 보면 Sensitivity나 call된 개수 및 false call이 비슷한 경향을 보이는 것을 확인할 수 있다. 시간만 된다면 각각의 툴들만이 call한 variants에 대해서 어떠한 특성을 가지는지 확인하고 싶지만, 필자는 더이상의 노력을 할 수 없음을 널리 이해해주기 바란다.

Strelka 일루미나에서 제대로 만든 S/W

타 툴들에 비해서 엄청나게 빠른 속도(multi-thread지원)를 자랑한다. 앞선 S/W들과는 달리 Java가 아니다. 논문은 "Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs"로 앞선 논문들이 언급한것과 같이 MuTect와 비슷한 성능을 보여주며, Virmid와 MuTect와 concordance를 비교한 Venn Diagram은 각각 다음과 같다. 아무래도 MuTect가 가장 많은 SNV를 찾았고 Virmid와 Strelka와 겹치는 부분이 상당하다. 

Virmid와 Strelka

MuTect와 Strelka

3개의 툴 비교

맺음말

지금까지 살펴본 결과를 보면 "Detecting somatic point mutations in cancer genome sequencing data: a comparison of mutation callers"를 보더라도 MuTect의 경우 low allelic-fraction에 강하며, VarScan2의 경우 high-quality sSNVs에 강한면을 보이고 있다.


아울러, 논문들을 리뷰한 결과 Strelka 역시 low-allelic-fraction에 강한 성능을 보이는 Strelka를 포함하여 Virmid, MuTect, Strelka 3개의 툴의 교집합과 각각의 툴들만 고유하게 calling한 variants에 대해서 features를 선별하여 훈련시킨다면 아마 가장 강력한 low-allelic-fraction에서의 Somatic Mutation 툴이 되지 않을까 한다.


어째 글이 쓰다만거 같지요. 끝

저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.09.04 18:27
Currently 3 comments want to say something now?
RNA-Seq Applications
2014.04.17 14:55 | 바이오인포매틱스
RNA-Seq을 수행하면 다양한 정보를 얻을 수 있습니다. 그중 첫번째 Annotation은 크게 Alternative Splicing Events와 Identify Known and Novel Transcripts입니다.

1. Alternative Splicing Events

1,2,3,4,5,6의 총 6개의 exon이 존재하는 gene이 있는 경우 genomic DNA에 read들을 매핑한 결과가 다음과 같다고 하자. paired-end read는 read간에 '---' 대시로 서로의 연결을 보여주고 있다. 맨 하단의 read 2개는 대시외에도 붉은색 원으로 보이는 부분은 read가 서로 끊겨 있다. 즉, 1,2,3이 연결되어 있고 5,6이 연결되어 있음을 알 수 있다. 그리고 4,5번에 걸친 read를 통해 현재 read의 정보를 종합하면 1,2,3,4,5,6이 모두 연결된 하나의 isoform이 존재한다는 것을 알 수 있다. 즉, a,b,c,d 4개의 read 정보를 통해 1,2,3,4,5,6으로 구성된  Isoform #1이 존재하는 것을 확인 할 수 있다.


또 다음과 같이 genomic DNA에 매핑된 read들이 있다면, 어떨까? 2번과 4번의 exon에는 read 정보가 존재하지 않으며, 붉은색 원의 정보를 통해 1,3번 exon과 3,5번 5,6번이 서로 연결된 것을 확인 할 수 있다. 최종적으로는 1,3,4,5번으로만 구성된 새로운 Isoform #2가 존재한다는 것을 알 수 있다. 

Ozsolak, F. and Milos, P. RNA sequencing: advances, challenges and opportunities Nature Review Genetics (2011)

그렇다면 어떻게 붉은색 원의 read들을 reference에 매핑할 것이냐는 문제가 남는다. 일반적인 mapping 툴은 붉은색 원과 같은 splice juction을 고려하지 않는다. 따라서 TopHat (A spliced read mapper for RNA-Seq)과 같은 툴을 이용해야 한다. 이제 일반생물학 시간에 배운 내용을 되짚어 보면 DNA는 intron 부분이 제거되고 exon부분만 서로 연결되며 이때 capping과 poly-A tail이 붙은 mRNA로 변하게 되며 intron 영역은 Donor와 Acceptor splice site, GT와 AG부분이 splicing 된다.


 

TopHat의 알고리즘은 gt와 ag를 보고 가능한 주변 exone간의 가능한 splice를 만들고 read들을 splice에 붙이게 된다. 

Trapnell, C., et al TopHat: discovering splice junctions with RNA-Seq Bioinformatics (2009)

2. Identify Known and Novel Transcripts

기존의 알려진 exon/gene에 매핑된 것이외의  mapping된 read들은 새로운 exon이나 gene일 수 있으며, mapping 되지 않은 read들은 새로운 splice junction일 수도 있다. 

일반적으로 reference에 기반한 transcriptome을 발굴하는 과정은 첫번째 그림과 같으며, 실제로 blat, TopHat, GSNAP, SpliceMap, MapSplice등이 splice를 고려한 align을 제공하며, TopHat에 의해 매핑된 read는 cufflinks를 통해 기존의 알려진 또는 새로운 transcripts들을 알아내게 된다. 아래 그림의 a,b,c를 참고하기 바란다. 


Jeffery A. Martine and Zhong Wang Next-generation transcriptome assembly Nature Reviews Genetics (2011)

Trapnell, C., et al Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isofrom switching during cell differetiation Nature Biotechnology (2010)

3. Quantification: Expression Profiling

이번에는 RNA-Seq의 3번째 응용으로 발현량을 profiling하는 방법입니다. NGS의 특성상 read를 통해 expression을 profiling해야 하고 이때 여러가지 고려할 사항들이 있다. 우선 Normalization에 관한 부분으로 많은 read가 존재한다는 것은 1) transcript의 사이즈가 길거나 (long) 2) 높은 depth of coverage로 시퀀싱 되어서 일 수 있다. 따라서 이러한 서로 다른 transcript의 length와 비교하고자 하는 각기 다른 condition으로 부터의 총 sequnce를 normalization해야 정확한 expression을 profiling할 수 있다.

이래서 나온 것은 바로 RPKM: Read Per Kilobase per Millon mapped reads 이다. RPKM=C/LN, C: Number of mappable reads on a feature, L: Length of feature (kb), N: total number of mappable reads (millions)

아래와 같이 kidney와 liver 두 샘플에서 ENSG00000212679 유전자의 read 수가 각각 620, 746이며, Kideny와 Liver의 총 sequencing depth가 9,293,530과 8,361,601이며, 해당 gene의 length가 1500 bases인 경우 각각의 RPKM 값은 44.48과 59.48이 된다. liver 너 이김. 이러한 normalization 방법으로는 RPKM외에도 RPKM이 read를 사용하는 대신 fragment (paired-end인 경우를 고려해서)를 사용하는 FPKM, Upeer-quartile, TMM 등이 있다. 이러한 quantify expression은 cuffcompare나 cuffdiff, R 패키지의 edgeR, DESeq 등을 이용하여 분석이 가능하다. (위 그림의 a,d,e 참고)



추가) 아래와 a와 같은 구조의  유전자가 있을 경우 TSS(alternative transcription start site)에 따라서는 A,B와 C 총 2가지 경우, CDS에 따라서는 A와 B,C 2가지 경우가 있다. 

isoform에 따라서는 A,B,C 모두를 비교가 가능하며, (a)
isoform내에서 유전자가 서로 다른 regulated를 보이는지를 비교하는 것은 A와 B만 가능 (b)
splicing의 differential을 보기 위해서는 condition A,B에 대해서 각각 A와 B를 비교 가능 (c)
TSS에 따른 promoter의 성능을 보려면 A+B와 C를 비교하는 것이 가능 (d)
promoter의 서로 다름을 비교하는 경우 condition A,B에 대해서 각각 A+B와 C를 비교 가능 (e)
유전자의 CDS의 결과물 비교는 A와 B+C에 대해서 가능(f)
protein 산물의 차이 비교는 condition A,B에 대해서 각각 B+C와 A를 비교 가능 (g)

Roberts et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks Nature Protocols (2012)

Qi Liu et. al RNA-Seq data analysis at the gene and CDS levles provides a comprehensive view of transcriptome response induced by 4-hydroxynonenal Mol. BioSyst (2013)

4. Gene Fusion

이제 마지막은 gene fusion이다. 가장 잘 알려진 fusion gene은 ABL1과 BCR gene으로 다음과 같이 두개의 유전자가 서로 fustion되어 질병이 발생할 수 있다. 현재는 Tophat2에 흡수된 tophat-fusion으로 분석이 가능하다. 

Ozsolak, F. and Milos, P. RNA sequencing: advances, challenges and opportunities Nature Review Genetics (2011)


아래는 TopHat-Fusion pipeline과 break point를 찾아 내는 그림이다.



Daehwan Kim, Steven L Salzberg TopHat-Fustion: an algorithm for discovery of novel fusion transcripts GenomeBiologoy (2011)

지금까지 초간단 RNA-Seq으로 할 수 있는 것들을 알아 보았다. 자세한 내용은 언급된 논문들을 살펴보면 되고 실제 분석에 관한 내용은 TopHat을 이용한 턱시도 프로토콜을 기반으로 나중에 하는걸로 하겠다.
 
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.04.17 14:55
Currently 댓글이 없습니다. comments want to say something now?
somatic mutation 찾기
2014.04.08 08:58 | 바이오인포매틱스
저번 포스팅에서는 일반적인 snp/genotype calling 메소드에 대해서 알아보았다. 이번에는 cancer분석에서의 somatic mutation 분석에 대해서 살펴보도록 하자. 이번 포스팅에서는 "Virmid: accurate detection of somatic mutations with sample impurity inference"라는  논문을 사용?할 것이다.

일반적으로 암 분석을 한다는 것 즉 somatic mutation을 찾는것은 variant calling의 하나로 NGS가 clinical로 가기 위한 기본적인 단계라고도 할 수 있다. somatic mutation을 찾는 전통적인 방법은 샘플 (normal/disease 또는 normal/cancer 또는 control/mixed disease sample이라고 각각 부른다.)에 대해서 각각 variant를 calling하고 두 샘플에서의 variants의 차이를 비교하는 것이다. 이는 단순한 비교를 통한 분석으로 최근에는 disease-control의 두 샘플을간의 joint probabilities를 계산을 통해서 true somatics mutation을 분리한다. 이와 더불어 de novo mutations을 찾기 위한 다양한 연구들이 진행되었다. 위에서 언급한 단순 비교를 지원하는 Traditional approaches 툴로는 VarScan이나 SomaticSniper가 있으며, joint probabilities를 지원하는 툴로는 VarScan2나 JointSNVMix가 있다. 
뭐 어찌되었든 간에 cancer 분석에 있어서 다양한 시도들이 이루어지고는 있는 상황에서 somatic muntation 분석에서의 가장 큰 걸림돌은 disease sample이 지니고 있는 impurity와 heterogeneity인한 분석 결과의 정확성 문제이다. 이러한 문제를 해결하기 위한 방법으로 나온것 중 하나가 바로 샘플의 contamination level을 계산하고 이를 분석에 적용하는 것이다. 바로 이러한 개념을 구현한 것이 오늘 이야기할 논문이다. 

샘플의 contamination level 계산하기

heterozygous mutation이 있는 경우 이 mutation에서 B allele frequency (BAF)는 50%에 가깝게 된다. BAF는 Alt와 Ref의 커버리지의 합을 Alt의 커버리지로 나눈 값이되며 BAF=coverage Alt / (coverage Alt + coverage Ref)가 된다. 암샘플이 정상샘플과 contamination이 심하게 되었다면 높은 BAF값을 보이게 되고 이를 통해 암 샘플이 정상 샘플이 contamination 되었다고 볼 수 있게 된다. 이렇게 mixed disease 샘플에서의 control 샘플의 비율 (α, 0 ≤ α ≤ 1)을 고려한 것이 VarScan2나 Strelka가 contamination level (α)를 고려하거나 이와 유사한 컨셉을 가지고SNV calling을 하고 있다.

Virmid (virtual microdissection for variant calling)

앞서 봤듯이 contamination 문제를 해결하기 위한 방법으로 sample contamination level (α)을 통해 somatic mutations에서의 disease genotype을 추정한다. 이를 통해 분석에 있어서 시퀀싱 에러, 매핑 에러, read mapping bias와 CNV의 영향 등의 bias를 줄일 수 있다. virmid를 이용한 분석 workflow는 input으로 short read sequence를 받는데, 하나는 순수한 control 샘플에 대한것과 다른 하나는 mixed되었다고 추정되는 disease 샘플이다. reference mapping의 preprocessing 단계를 거쳐 GATK를 이용하여 indel realign 등의 일반적인 과정을 수행한다. 이 작업이 끝나고 나면 모든 포지션에 대해서 BAF (B allele frequency)를 계산한다. BAF후에는 초기 필터를 적용하는데 이는 quality control 뿐 아니라 샘플 크기를 감소하기 위해 적용된다. (control sample은 B allele가 0인 것만 사용하고 disease sample은 1개 이상의 B allele가 있는 것을 사용하여 이 샘플들은 true somatic mutation을 찾기 위해 최소 BAF를 이용하여 필터링한다.) 마지막으로 두개의 control과 disease의 필터된 pileup 포맷의 alignment file이 input 파일로 준비된다. 이렇게 input file이 준비되면 이제 virmid는 contamination level (α)을 추정하는 과정을 수행한다. 이 값은 전역 파라미터로 모든 포지션에 대해서 동일하게 적용된다. 이 α는 작은 포지션이도 충분히 추정이 가능하게 된다.

Overall Virmid workflow

disease sample에서 11개의 각기 다른 비율로 두개의 genome이 섞인 시뮬레이션 데이터를 이용하여 (α = 1%, 5%, 10%, 20% ~ 90%) Virmid를 사용한 것과 그렇지 않은 (call-based method)의 contamination의 level의 추정한 결과 실제(검은색)와 virmid(빨간색)에 비해 call-based의 경우 60 ≤ α에서 성능이 떨어지는 것을 확인 할 수 있다.


실제 TCGA(The Cancer Genome Atlas)의 15개의 exome breast cancer 데이터를 이용하여 Strelka, MuTect, JointSnvMix2(JSM)를 이용하여 비교한 경우는 다음과 같다.


결론?,,,벌써?

짧게 나마 아주 짧게나마 cancer 분석을 위한 순수 국내기술?의 Virmid에 대해서 알아보았다. 어차피 정답은 없는법 한번 써보는 수밖에는 ㅋ Virmid는http://sourceforge.net/projects/virmid/
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.04.08 08:58
Currently 3 comments want to say something now?
베이즈 정리를 정리하고 넘어가자
2014.03.26 19:51 | 바이오인포매틱스
바로전에 포스팅한 variant calling에 대한 것에 후속으로 somatic mutation에 대한 내용을 정리하고 있다. 살짝 귀뜸해주면 VarScan이라는 툴에 대한 논문을 보려다가 "Virmid: accurate detection of somatic mutations with sample impurity inference"라는 논문으로 급선회했고 virmid는 간단히 cancer 분석에서 contol sample이 mixed된 disease sample에 대해서 control sample이 어느정도의 proportion을 차지하는지를 estimate하여 이것을 somatuc mutation을 calling하는데 사용하는 논문이다.

뭐 그거 그렇고,  바로 snp/genotype/somtic mutation calling에서 가장 중요한 역할을 하는 것이 바로 베이즈의 정리이고 이것에 대한 이해가 선행되어야 한다. 그래서 했다. 국내 번역본으로는 "집단지성프로그래밍"이라는 책에 바로 베이즈의 정리에 대한 내용이 나오는데 오늘은 genome data가 아닌 스팸데이터를 통해 베이즈의 정리가 무엇인지 한번 살펴보려고 한다. 벌써 번역서가 나온지도 08년도이니 6년이나 지난 책이 되어 버렸고 이책이 나에게 좀 더 친근하게 다가온것은 원저자가 생명공학회사에서 약제 발현 원리를 이해하기 위한 알고리즘 개발과 데이터 마이닝을 했던 사람이라 더욱 특별했다.

집단지성프로그래밍
카테고리 컴퓨터/IT > 프로그래밍/언어
지은이 토비 세가란 (한빛미디어, 2008년)
상세보기

우선 우리는 분류기에 대해서 알아볼 것이다. 이 이메일이 스팸인지 아닌지? 특정부분의 genotype이 AA, AB, BB인지?와 같은 분류 말이다. 분류를 하기 위해서는 우선 훈련을 시켜야 한다. 이미 스팸 또는 정상의 이메일을 가지고 훈련을 시킬 것이고 이 훈련은 고되지만 훈련량이 많을 수록 즉, 이미 분류된 문서가 많을 수록 추후 이 분류기가 올바르게 판단할 가능성이 높아진다.

우선 feature들이 각 분류에서 얼마나 나타나는지를 확인할 것이다. 예를 들자면, {'quick':{'bad':0, 'good':2}, 'the': {'bad':3, 'good':3}}과 같이 각 단어(feature)가 두개의 good, bad라는 분류에 얼마나 나타나는지를 feature count(fc)를 한다. 

그럼 이제 특정 분류에 특정 단어가 나타날 확률을 계산할 수 있다. Pr(단어 | 분류)라는 조건부 확률, 예를 들어 "Good으로 분류된 문서가 단어 quick을 가질 기회"인 Pr(quick | good)은 good이라는 문서에 quick이라는 단어가 나타난 횟수를 good이라는 분류에 있는 전체 문서 개수로 나눈 값이 된다. 이 값이 0.666이라면 2/3만큼의 기회가 있다는 것을 의미한다.

지금까지는 Pr(Word | Category), 즉 특정 카테고리에 대해  특정 Word가 나타날 확률만을 알고 있는데, 우리는 특정 단어를 분류하려는 것이 아니라 하나의 문서를 분류하려고 하는 것이다. 즉, Pr(Document | Categroy)를 알고 싶은 것이다. 그렇다면 우리는 간단히 문서내의 모든 단어들의 확률을 곱하여 알 수 있다.

특정 분류에 대해 문서가 속할 확률 (Pr(Document | Categroy)) 을 알아내었다. 하지만 우리가 하려고 하는 것은 특정 문서가 있을때 이것을 각각 어떻게 분류할지에 대한 확률을 알고 싶은거다. 즉 Pr(Category | Document)를 알아야 하는거고 그래야 스팸 필터기를 만들 수 있는 것이다.

젠장 뭔지는 모르겠지만 feature를 뽑아서 이리저리해서 "good으로 분류된 문서가 특정 단어들로 구성된 문서를 가질 기회"를 구했는데, 정작 필요한 것은 "특정 단어들로 구성된 문서가  good에 속할 확률"이라니... 자 이렇듯 조건부 확률을 부침개 뒤집듯 확 뒤집어야 하는 상황이 된 것이다. 그럼 뒤집자! 어차피 이 뒤집기는 250년전 영국의 수학자 토마스 베이즈가 해결했다.

Pr(A | B) =  Pr (B | A) * Pr(A) / Pr(B)로 일반적으로 표현하며 이는 Pr(Category | Document) = Pr(Document | Categroy) * Pr(Category) / Pr(Document)로 표현할 수 있다. 우리는 Pr(Category) 즉, 무작위로 선택한 문서가 이 분류에 속할 확률을 알아야 하는데 이건 분류내 문서 건수를 전체 건수로 나누면 되고 Pr(Document)는 생략,,,

마찬가지로 genotype/snp call에서도 문서(Documnet) 대신 "특정 site의 read정보들"을 사용하고 category (good/bad) 대신 "AA, AB, BB"냐의 genotype을 사용하는 것이 다를 뿐 기본적인 베이즈의 정리를 따른다.

우리가 스팸필터기의 성능을 높이기 위해서는 feature로 뽑힌 단어를 단순히 단어들만이 아닌 대소문자를 인식한 단어, 단어와 단어의 조합인 구나 절을 사용하여 성능을 높일 수 있듯이 genoytpe/snp call에서도 Pr(B | A) 즉, 흔히 말하는 genotype likelihood를 정교하고 세련되게 만듦으로서 genotyping의 성능?을 높일 수가 있다. 바로 이부분이 앞으로 챌린징이 가능한 부분이다.

다소 생략된 느낌이 있다면 책사서 보길 바람요. 일반적인 데이터 대신  genome 데이터를 가지고 이런 입문서 만들면 재미있을듯한데... 팔리려나??

저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.03.26 19:51
Currently 2 comments want to say something now?
NGS 데이터에서의 Genotype and SNP calling
2014.03.22 18:41 | 바이오인포매틱스
지금까지 수천 샘플에 대한 genotype/snp calling을 수행했음에도 2011년도에 발표된 "Genotype and SNP calling from next-generation sequencing data"라는 리뷰 논문을 이제서야 꺼내어 읽어본다. 이 논문에 대한 내용은 이미 "ideas should be in papers" 블로그에 소개되었으나 나름 다시 정리하는 차원에서 여기저기 살을 붙여서 작성했다. 간혹 이해가 가지 않는 부분은 참고하여 작성했다.

전체적인  genotype/snp calling 분석

Base calling

genotype/SNP calling에 있어서 가장 기초가 되는 것은 per-base quality score로 이는 일반적으로 NGS 장비의 기본 base calling software를 통해 제공된다. Illumina의 경우 RTA (Real Time Analysis software)를 사용하여 base calling을 수행하며, miss call error rate가 1%정도로 합성되는 시점에서 sync가 맞지 않아 signal이 맞지 않아서 발생하게 된다. 금방 언급했듯이 가장 기초적인 데이터로 Illumina의 base caller인 Bustard 대신 Ibis, BayesCall을 사용하게 되면 error rate가 5~30% 정도 향상된다고 한다. 

Summary of the available applications used for base-calling on the Illumina platofrm
Brief Bioinformatics(2011), Base-calling for next-generation sequencing platforms

 
Comparison of Base-Callers for the Illumina Platform 

Recalibration of per-base quality score

앞서서 base quality score는 genotype/snp call에서 가장 기본이라 했다. 여기서는 base calling 과정에서 생산된 raw Phread-scale quality score는 정확하지 않기 때문에 여러가지 covariate를 이용하여 재조정되어야 한다. 이러한 역활이 GATK의  BQSR (Base Quality Score Recalibration)이 담당하고 있다. 
  • Reproted quality score
  • Position within the read (machine cycle)
  • Preceding and current nucleotide (sequencing chemistry effect, dinucleotide context)
이러한 base quality를 조정하는 또 다른 방법은 Heng Li의 BAQ (Base alignment quality) 간단하게 INDEL 주변의 SNP들을 재조정하는 것이다. 여기서 잠깐 뒤에서 알아보겠지만, GATK나 Samtools에서의 genotype/snp calling은 기본적으로 bayesian model을 기반으로 하고 있지만, 그 결과는 slightly different하다는 것이다. 이러한 차이는 직접  Heng Li (아이디 lh3)가 언급한 내용을 정리하는 다음과 같다.

1. Alignemnt preprocess 부분

  • GATK는 low mapping quality를 가지는 read를 사용하지 않는 반면 samtools는 기본적으로 모든 read들을 사용한다.
  • Samtools는 base quality에 alignment quality를 덧씌워 사용하며, GATK가 현재 그렇게 하는지는 모르겠다.
  • GATK는 read pair에서 끝(ends)이 중복되는 경우 축소(collapses)하게 된다. 이게 맞는 것 같다. 하지만, 분석 정확성이 크게 향상되지 않을 것 같아 이것을 samtools에 구현하지 않았다.
  • samtools는 BAQ를 기본으로 사용하는 반면 GATK는 그렇지 않다. (앞서 언급했듯이 GATK에서도 구현되었음)

2. SNP genotype likelihood model 부분

  • GATK는 sequencing error가 독립적이라고 가정하는 반면 samtools는 두번째 에러 (second error)는 더 높은 기회를 가진다고 본다. 이러한 차이는 depth가 수백x가 되지 않는 이상 거의 영향을 미치지 않는다. 

3. Indel genotype likelihood model  부분

  • GATK와 samtool의 가장 큰 차이로 samtool은 BAQ를 기반으로 heuristic한 방법을 allele를 구분하는데 사용한다. GATK는 Dindel 모델을 기반으로 시간이 지나면서 많은 변화가 있으며 나(헨리)는 정확히 어떻게 동작하는지 모르겠다. 

4. Variant calling model 부분

  • samtools과 gatk 모델이 거의 동일하며, 최초 Rasmum Nielse과 Richard Durbin이 제안한 것으로 samtools 논문에 언급했으며 gatk 동일한 공식을 사용하여 구현했다. 추후 gatk는 multi-allelic 모델을 추가했으나 samtools는 포함하지 않았다.

5. 필터링 부분

  • 또다른 큰 차이로 samtools는 사람이 손으로 (hand-tuned)  필터를 사용하는 반면 gatk는 데이터로부터 filter를 학습한다.
  • 물론 model을 훈련하기 우히나 충분한 데이터와 non-polymorphic sites가 필요하지 않은 human variant calling에 있어서 gatk가 좀 더 편리하고 강력하다.
  • 게다가 효과적인 필터링을 위해 gatk는 samtools 보다  haplotype score와 같은 좀 더 많은 통계정보를 수집한다.
  • gatk가 진단의 목적 (diagnostic purposes)에 매우 유용하다.

6. GATK만의 특징

  • haplotype caller나 indel realignment 등
  • haplotype caller는 gatk v2의 가장 중요한 특징으로 longer indels에 대한 높은 민감도를 보여준다.

7 samtools만의 특징

  • genotype-free analysis와 physical phasing 등

8. technical 차이

  • gatk가 더 정교하고 더 복잡한 것은 사실이며 그럼에도 불구하고 samtools는 훨씬 간단한 프레임워크로 거대 데이터 셋 (huge data sets)에 동일하게 확장이 가능하다.

BAQ

초기 genotype calling 방법

초창기 시절의 genotype calling 방법을 알아 보기전에 우선 SNP calling과 genotype calling은 서로 다르다는 것을 염두에 두어야 할 것이다. snp calling은 polymorphisms이 있는 base를 찾는 것이며, genotype calling은 앞서 snp 부분에 대해서 개인의 genotype을 결정하는 것이다. 초창기 방법은 누가 들어도 충분히 구현이 가능한 방법으로 high-confidence bases (Q20) 남기는 fitering을 한 후 관찰된 allele를  고정된  cutoff (ex, heterozyous가 20-80%의 proportion을 갖는지?)를 사용하여 판단하는 것이다. 이러한 방법은 sequencing depth가 높은 (>20x) 경우 비교적 잘 작동하며 이를 software들로는 Roche's GSMapper, CLC Genomic WOrkbench, DNASTAR lasergene software들이 있다. 중요한 것이 바로 cutoff값인데 이는 경험적인 cutoffs를 사용하여 성능을 향상할 수가 있다.

확률론적 방법

앞서 살펴본 초장기 방법은 moderate 또는 low sequencing depth인 경우 heterozyouse gentotype에 대해 잘 작동하지 않으며 게다가 quality score에 기반한 filtering은 정보 손실로 이루어질 수 밖에 없다. 추가적으로 genotype을  추론하는데에 있어서 불확실성에 대한 어떠한 측정값을 제공하지 않는다는 것이다.
이러한 이유로 각각의 genotype의 quality score를 사용하여 사후확률 (posterior probability)를 제공하는 확률론적 방법이 개발되어 널리 사용되고 있다.

베이지안

경험을 통해 얻은 확률과 특정 선입견 (class)하에 들어온 관측치가 그 선입견 (class)에 들어맞을 확률을 통해 들어온 관측치가 특정 class일 확률을 계산하는 것이다. 예를 들어 아이가 침대에서 떨어질 수 있다는 특정 사건 (class)가 있다고 한다면, 침대 끝에 있다는 현재 상황에서 떨어질 수 있다는 class는 꽤 높은 사전확률 p(뒤로 떨어짐)을 갖는다. 이 class의 likelihood는 p(목이 뒤로 넘어감 | 뒤로 떨어짐)을 통해 목이 뒤로 넘어가는게 맞다라고 판단할 수 있다.  

바로 genotyp calling에 이 베이지안 방법을 적용하여 우리는 해당 position에 대해서 기존의 경험을 통해 해당 부분의 genotype에 대해 확률론적인 의미를 부여할 수 있게 된다. 

사전 확률 (Prior probabilty of genotype)

사건(class)에 대한 어떠한 정보도 없는 상태에서 뒤로 떨어지는 확률을 경험을 통해 알고 있다. 마찬가지로 genopype에서도 우리는 다양한 경험을 통해 정보를 알고 있다. 
  • genome wide SNP rate
  • SNP substitution type not equally likely
  • allele frequency
예를 들어, 다음과 같은 정보를 미리 알고 있다면 각 genotype에 대한 prior값을 계산할 수 있다.
heterozygous SNP rate: 0.001
homozygous SNP rate: 0.005
Reference allele: G
Transition/transverstion ratio: 2
라고 한다면 다음과 같은 prior probability을 갖게 된다. G/G가 나올 확률이 가장 높다는 것을 사전에 우리는 알게 되었다.


Use dbSNP prior probability

이외에도 다른 정보(reference sequence, SNP data-bases, available population sample) priors에 사용할 수 있다. 아래 그림은 dbSNP의 G/T site에 대한 prior probability이며(각 genotype에 대해 GG와 TT는 0.454, GT는 0.0909, 그외에는 10-4 이하)  G/G, T/T의 확률이 같은 것, 즉 G/T hetero site라는 것을 알 수 있다. 이러한 방법은 SOAPsnp가 사용하는 방법으로 MAQ또한 유사한 방법을 사용한다.

Use differnet polymorphism rate for different genomic regions
Consider different Ti/Tv rate for exonic regions


위의 경우에는 하나의 샘플 (single sample)에 대한 prior값을 구했지만, multiple sample을 고려하는 경우에는 prior값을 좀 더 improve할 수 있다. (추가해야함)이제 genotype G에 대한  사전확률 p(G)를 알았으니, genotype에 대한 likelihood p(X|G)값을 계산하기만 하면 우리는 Bayes' formula를 이용하여 p(G|X)를 계산할 수 있게 된다.

Genotype Likelihood

앞서 이야기 했듯이 read의 quality score를 이용해서 genotype likelihood 값을 계산한다. 한 개체에서 나온 특정 genomic site의 read i의 genotype G에 대한 확률 p(Xi|G)는 Xi의 quality score의 rescaling으로 계산이 가능하다. 그리고 genotype likelihood p(X|G)는 이 p(Xi|G)들의 product로 계산되어진다.

genotype likelihood를 계산하고 genotype prior(P(G))을 이용하여 Bayes' fomula로 각 genotype G에 대한 posterior probability (P(G|X))를 계산해서 최대의 확률을 갖는 genotype,G1을 그 개체의 genotype으로 할당한다. 이때 P(G1|X) 혹은 P(G1|X)와 P(G2|X) 사이의 비율(G2는 posterior probability가 두번째로 높은 genotype을 genotype confidence로 사용한다.

LD 정보를 포함한 calling

population genetics 관점에서 imputation은 SNP data set에서 missing data를 찾는 것으로 간단히 어떤 population에서 ATA, CGC라는 두개의 haplotype만 관찰된다고 하자. 어떤 개체의 샘플에서 첫번째 site에서 A나 C가 나오고 세번째 site에서 A와 C가 나온다면 이 샘플의 두번째 site를 모른다고 한다면 손쉽게 우리는 중간의 빈곳이 T또는 G라는 것을 알 수 있다. 바로 이것을 NGS 데이터에 접목하면 되는 것이다. 아래 그림은 다음의 3가지 방법에 대한 genotype call accuracy를 보여준다.
  • 각 개체에 대해 따로 genotype calling을 수행 (파란선)
  • LD기반의 분석을 사용하지 않고 모든 개체들을 jointly analysis (붉은선)
  • LD기반의 분석을 이용하여 모든 개체들을 jointly analysis (검정색)
multiple individual을 사용하는 것이 single sample에 비교해서 각각 80%에서 87%로 향상되는 것을 볼 수 있다. 약 40% 정도의 non-calls이 있다고 할때 LD 패턴을 고려하지 않는 분석과 LD를 고려한 분석이 비슷한 accuracy를 보이는 것으로 보아 이는 약 40% 정도의 missing data가 존재한다고 볼 수 있다. 이렇게 LD 패턴을 고려하는 경우 genotype calling에 있어서 많은 향상으로 보여주며 이는 hapmap이나 SeattleSNPs처럼 고품질의 LD정보가 있는 레퍼런스 데이터를 사용하는 경우 많은 도움을 받을 수 있다. (23andMe 데이터 또한 이러한 imputation을 통해 비어 있는 정보를 채워 넣을 수 있다) 하지만, reference panel에 없는 rare mutations에 대한 calling에서는 이러한 LD 정보는 성능 향상에 도움을 주지 못한다.

Filtering

각 site마다 사후확률이 정확히 계산되어졌다면야 필터링 등의 추가적이 데이터 조작은 필요없겠지만, 어디 그게 현실에 가능하단 말인가? 실제 데이터에서는 필터링을 통해 genotype/snp call이 향상된다. 1000 genomes project의 경우 이미 genotype이 알려진 HapMap 데이터와의 높은 차이로 인해 모든 sequencing batch를 제거해야 했다. 이러한 필터링은 genome-wide SNP genotyping이 이미 수행된 프로젝트에서나 가능하다. 인간 이외의 organisms이나 genotyping 데이터가 없는 경우라면 1000 genomes project 보다 높은 에러율이 나타날것으로 예상해야 할 것이다.

다른 형태의 필터링으로는 HWE로 인한 devaiations (일반적으로 low-quality scores, major/minor allele에 대한 sysmatic differnece in quality scores, LD aberrant, 너무 높은  read depth, strand bias 등)가 genotype과 snp calling의 정확도를 높이는데 활용될 수 있다.  적절한 필터는 시퀀싱 프로토콜과 upstream 분석에 의존한다. 예를 들어 strand bias (plus와 minus strand가 불균형적인 부분)가 존재하는 부분에서는 에러가 나기 쉬운 부분으로 이 부분은 필터 아웃 되어야 한다. exome 캡처를 사용한는 captured 시퀀스인 경우 이러한 bias는 capture array의 artefact로 인할 가능성이 크다.

Incorporating genotype uncertainty in analyses of NGS data

NGS 데이터 분석에서의 genotype에 대한 불확실성을 아예 안고 살아라, 즉 분석에 이러한 불확실성을 포함해서 분석을 하라는 것입니다. 다른 application은 서로 다른 genotype calling 전략을 사용하라는 것이다. 여기서는 association-mapping studies와 population-genetic studies에서 이러한 불확실성에 대한 이야기를 하고 있다. 

Recommendations

NGS 데이터 분석은 빠르게 변화하는 분야로 분석을 위한 새로운 통계 방법들이 계속 개발되고 있다. 더구나 1000 genomes project로 부터 많은 툴들이 개발되었지만 아직 publish되지 않은 것들도 많다. 이러한 시점에서 몇가지 사항에 대해 권장한다.

1. Base calling과 quality scores 계산은 테스트/벤치마크된 방법을 사용해라.
  • per-base quality scores의 recalibration은  GATK나 SOAPsnp를 추천
  • reference genome/short-read를 align하는 경우 Novoalign이나 Stampy와 같은 sensitive aligner를 사용해라.
  • 후에 bwa와 같이 효율적인 aligner와의 hybrid 사용 추천
2. SNP calling은 모든 샘플을 동시에 할 수 있는 것을 사용한다.
  • snp calling에 있어서 likelihood ratio나 베이지안을 사용하기 때문에 데이터로 부터 allele freqeuncy에 사전 분포를 사용하기에 모든 샘플을 동시에 분석하는 것이 좋다. 
3. genotype calling 또한 다수의 individual을 합쳐서 한다.

4. 가능하다면 LD-기반의 방법을 통해  genotype과 snp calling의 accuracy를 향상시킬 수 있다.

몇가지 추가적인 고정을 통해 성능을 향상 시킬 수 있는데, 가령 local realignments, 다수의 snp/genotype calling 알고리즘과 quality score기반의 post hoc 필터링을 사용할 수 있다. 추가로 아래는 1000 genomes project에서 variant/genotypes call의 성능 향상을 위해 사용한 방법으로 이미 위에서 설명한 방법들이 다수 포함되어 있다.
  • INDEL realignment
  • Per-base alignment quality (BAQ) adjustment
  • Robust consensus SNP selection strategy
    • variant quality score recalibration (VSQR)
    • support vector machine (SVM)
  • imporved genotype likelihood calculation
    • BAM-specific Binomial Mixture Model (BBMM)
    • leveraging off-target exome reads 
NGS 데이터에서 genotype/snp calling은 불확실성에 대한 확률적인 방법을 제공하는 것으로 allele의 count에 따라 간단한 방법에서 보다 정교한 많은 방법들이 존재한다. 확률적인 방법은 genotype likelihoods의 계산에 의존하므로 더 많은 연구가 필요하다. 또한 LD 기반의 방법의 개선을 포함한 다운 스트림 분석에서 genotype 확률을 통합하는 방법이 필요하다. NGS는 어찌되었든 향후 몇년동안 medical genetics 연구의 중심이 될 것이며 미래를 위해 충분히 가치가 있다.

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

좀 더 많은 genomes이 필요한 시대

$1000 게놈 시대가 진짜 도래했다. 이제까지 NGS 연구의 대부분이 하나의 genome 데이터를 가지고 연구(rare variant를 찾던)하던 것이 GWAS처럼 대규모의 cohort의 샘플을 수용하기 시작하면서 "Common Variant Association Study (CVAS)"에 눈을 돌리기 시작했다. 이는 가격뿐만 아니라 대량의 NGS 데이터를 다루기 위한 툴 또한 발전하면서 가능케 되었다.

이러한 CVAS 데이터는 cohort의 샘플들(individucal callsets)을 개별적으로 variant call을 하는 것이 아니라, joint callset을 만들어 joint variant discovery를 수행하여 흔히 말하는 power를 부여할 수 있다. 각각의 샘플 Sample #1을 분석해서 사용한다면 "A" site에 대해서 진짜 변이인지에 신뢰할 수 없다 하지만, Sample #1~Sample #N까지 모든 샘플을 보면 cohort내에서의 "A" site에 대한 real variation에 대한 좀 더 신뢰할 수 있게 된다.


또한 join discovry를 통해 bias 이슈 또한 어느정도 해결이 가능하다. 하나의 샘플만을 본다면, strand에 따라 FWD reads는 모두 REF와 같은 A(4개)를 REV reads는 REF allele A가 5개  VAR allele인 C가 3개 즉, allele ratio = 9: 3 = 3:1이 된다. 


자, 이제 joint analysis에서는 어떻게 systematic biases를 줄이는지를 보면, 여러샘플 데이터를 통해 allele bias와 strand bias를 filter out 할 수 있게 된다. 이것은 위와 같이 single sample에서는 해결 할 수 없는 부분이었다. 


근데 문제는 있다.

첫째 이렇게 joint analysis를 하기 위해서는 모든 site에 정보가 있어야 한다. 하나의 샘플만을 분석한다면 단지 variant에 대한 정보만 존재하면 되지만 joint calling을 위해서는 모든 샘플의 모든 site에 대한 정보를 필요로 하게 된다. 둘째, 현재까지의 접근으로 multisample을 분석할때 매우 comutationally intensive한 작업으로 computing cost와 "N+1 문제"에 봉착(모든 샘플을 동시에 분석하기 때문에 샘플이 하나 더 증가하면 +된 1개의 샘플을 포함해서 다시 N+1개의 샘플을 다시 분석해야하는...)하게 된다. 



기존의 여러개의 BAM(Analysis-Reday Reads)을 joint variant calling을 하는 경우 N+1 문제가 발생합니다. 물론 Reduced Reads를 이용해서 어느정도 computing resource에 대한 해결을 보긴 했지만,,, 새로운 방법은 각 sample별로 variant calling을 합니다. 이때 HC(Haplotyper Caller)의 gVCF mode를 통해 각 sample별 variant calling을 수행한 후 개별 gVCF를 통해 Joint Genotyping을 수행하면 기존의 Joint Varaint Calling을 동일한 결과물을 얻을 수 있게 된다.

gVCF

gVCF는 GATK 2.7버전부터 HC에 들어있는 기능으로, --emitRefConfidence (-ERC) 옵션을 통해 사용이 가능하며 이 옵션을 통해 이전에 언급한 모든 site에 대한 정보가 있어야 joint analysis가 가능하다라는 부분을 해결할 수 있다.
  • NONE: don't emit anythong (default) 아무것도 넣지 않는것으로 기본
  • BP_RESOLUTION: emit detailed information for each BP 각 Base Position의 모든 정보를 포함 (ref와 같던 같지 않던 모든 bp에 대한 정보를 포함하며 마치 UG(Unified Genotyper의 --output_mode의 EMIT_ALL_SITES와 유사한)
  • GVCF:  emit a block summarized version of the BP_RESOLUTION data BP_RESOLUTION의 모든 bp에 대한 정보를 넣는 반면 block 단위로 데이터를 압축(summarize)하는 것으로 마치 Reduced Reads가 불필요한(정보성이 없는) read에 대해서 필요한 정보만 남겨 데이터를 압축하는 것처럼, 일정 block이 reference와 같으면 해당 region 정보를 bp레벨로 다 정보를 입력하지 않고 합쳐서 보여주는 형태
마지막 GVCF를 다시 보면, 첫라인의 염색체 20번의 10000435 bp는 ref가 'T'이고 <NON_REF> 즉, 변이가 없고 해당 block은 길이가 4 즉, 10000435~10000438까지가 모두 Reference와 같다는 의미이다. (BLOCK_SIZE=4;END=10000438) 이 block은 대충 GQ가 비슷한 놈들끼리 묶이게 되며 MIN_DP와 MIN_GQ값으로 해당 block의 정보를 표현하게 된다.

20 10000435 . T . . BLOCK_SIZE=4;END=10000438 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:53:99:53:159
20 10000439 . T G 1553.77 . AC=2;AF=1.00;AN=2;DP=56;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=221.40;MQ0=0;QD=27.75 GT:AD:DP:GQ:PL 1/1:0,55:55:99:1582,164,0
20 10000440 . T . . BLOCK_SIZE=6;END=10000445 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:56:99:55:165

그럼, 왜 BP_RESOLTION이 있는데 왜 굳이  gVCF모드를 따로 두었느냐? 바로 데이터의 크기 때문이다. 모든 bp의 정보를 표현하다 보면 데이터가 커지게 되며 이는 곧 cohort내의 모든 샘플을 분석하는데에 비용이나 시간적인 문제로 직결되기 때문이다. 실제로 exome 데이터를 NONE, BP_RESOLUTION, GVF 모드로 각각 수행했을때 결과 파일의 용량을 비교해 보면 gvcf에 비해 bp resolution이 60배 이상 나는것을 확인 할 수 있다. 즉 몇천샘플을 가지고 분석을 한다면 gvcf가 훨씬 데이터 reduction에 의한 효과를 볼 수 있게된다.


gVCF 파일의 joint genotyping

각 샘플별로 HC를 통해 calling된 gVCF (genomic VCF, record every single locus in the genome, 각 샘플은 variant가 아닌 sites에 대해서도 reference에 대한 confidence estimation 정보를 갖게 됨(reference model)) 파일은 마지막으로  GenotypeGVCFs를 통해 하나의 joint된 VCF  파일을 만들게 된다. 이후 recalibration을 수행하면 끝

정리하자면, joint discovery(gVCF와 reference model을 통해)는 추후 지속적으로 늘어나는 샘플에 대해서 incremental하게 함으로써 cost를 줄이고, 기존에는 가능하지 못했던 large cohorts에 대한 study가 가능하다.

좀 혼동되어 사용된면이 있는데, 우선 나혼자라도 용어정리를 하고 난 주욱 이걸로 커뮤니케이션 하겠다.
  • individual callsets: 샘플별로 call된  variant들 셋트 (복수 결과들)
  • joint callset: 다수의 샘플이 한번에 call된 세트 (단수 결과)
  • joint variant discovery (join discovery): joint callset을 얻기 위한 일련의 과정(gvcf,,,reference model등을 활용)
  • join variant calling: 각 개별 샘플을 한번에 묶어서 calling하는 방법(sample #1 + ... + sample #N = one bam)
  • joint genotyping:  각 개별 샘플별로 calling 한 후 이를 묶는것 (sample #1, sample #N---->joint genotyping)
  • multi-sample calling: 다수의 샘플을 joint variant calling이나 joint genotyping등의 기법을 통해 분석하는 것
  • reference model: reference와 같은 site에 대해서 reference에 대한 confidence est 정보
     
참고
GATK HaplotypeCaller document
GATK blog - The GATK reference model pipeline for incremental joint discovery, in full detail
GATK Using the reference confidence calculation model & generating a GVCF document
GATK GenotypeFVCFs document
GATK presentations in DropBox

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

Clinical Sequencing을 위한 준비 - 표준 variants

미국 NIST (National Institute of Strandards and Technology)의 Div. Biosystems and Biomaterials 에서는 추후 임상으로서의 NGS 데이터 적용을 대비하기 위한 작업을 하고 있다. 이와 유사하게 국내에는 한국표준과학연구원 국가참조표준센터 (NCSRD)이 생명과학 관련 참조 표준 제정하고 있는데요. 아직 NGS 데이터와 관련한 표준은 없는 상태이다. (하단의 생명과학 관련 참조 표준 목록 참고) 

NIST는 NA12878에 대해서 자세한, 표준의 variants call set을 만들어 calling 알고리즘의 벤치마크나 기타 여러 분야에 활용할 수 있도록 하고 있다. 이 작업은 NIST, "Genome in a Bottle Consortium (GIAB) ", Harvard School (Dept. of Biostatistics) 공동으로 진행하였으며, 해당 결과는 CGAT (Genome Comparision & Analytic Testing)을 통해 활용이 가능하도록 하였다.

이와 관련되어 진행된 내용은 Blue Collar Bioinformaticsbcbio-nextgen에 설명되어 있으니 참고해보면 된다. 또한 이 작업은 Nature biotechnology에 "Integrating human sequence data sets provides a resource of benchmark SNP and indel genotype calls"에 논문으로 출판된 상태이다.

앞서 언급했듯이 NIST는 잘 알려진 CEU female NA12878에 대한 11 whole genome, 3 exome data set을 5개의 서로 다른 sequencing platform들을 가지고 Highly confident set of reference call를 생성했다.  (각 데이터셋은 GIAB의 ftp site를 통해 다운로드 가능) 서로 다른 sequencing platform과 variant caller로 인해 서로 다른 결과값들을 보여주고 이들간의 bias를 줄이기 위한 여러가지 방법이 동원되었다. 아래 표는 이러한 표준을 만들기 위해 사용한 플랫폼과 mapping algorithm들을 보여주고 있다. 최종적으로 다양한 sequencing technologies를 통합, bias를 최소화하고 machine learning을 적용하여 NA12878에 대해 2,741,014 SNP와 174,718 indels을 정의했다.

Referece call을 활용한 Calling method 평가

이렇게 만들어진 reference call을 가지고 FreeBayes, GATK UnifiedGenotyper, GATK HaplotypeCaller, 그리고 Ensembl calling method의 총 4가지 calling method들을 평가를 수행한 결과에 대해서 살펴 보도록 하겠다. 우선 각 calling method에 대한 간략하게 정의하면 다음과 같다.
  • Bayesian calling approach that handles sumultaneous SNPs and indel calling via assessment of regional haplotypes.
  • GATK UnifiedGenotyper: Bayesian approache to call SNPs and indels, treating each position independently.
  • GATK HaplotyperCaller: Performs local de-novo assembly to call SNPs and indels on individual haplotypes.
  • Ensembl calling method: Consensus approaches (하단의 "Machine Learining(SVM Ensemble)을 이용한 variants calling 성능 향상" 참고)
Blue Collar Bioinformatics의 글에 따르면, 각각의 평가 결과는 reference와의 concordant, dicordant에 따라 다음의 4개의 category로 나누어 평가를 진행하였다고 한다.
  • Concordant: reference material과 일치
  • Discordant (missing): NA12878의 reference에는 존재한 평가 데이터 셋에는 없는 변이 (these are potential false negatives)
  • Discordant (extra): 평가 데이터 셋에는 존재하나 reference에는 없는 변이 (these are potential false positives or missing calls from the refrerence materials)
  • Discordant (shared): 평가 데이터 셋과 reference에 모두 존재하나 서로 다른 변이 (allele differences, such as heteo vs. homo calls, or variant identification differences, such as indel start and end coordinates)
각 calling method에 대한 평가 결과는 다음과 같다.
  • FreeBayes: SNP과 indel 모두에서 GATK보다 우수한 성능을 보여줌, 최근 버전은 sensitivity와 specificity가 향상되어 GATK HaplotpyeCaller와 동등한 위치를 가짐, het/hom call의 경우 더 정확하여 이는 discordant shared variants의 lower number가 말해줌
  • GATK HaplotypeCaller: UnifiedGenotyper에 비해 전반적으로 나은 성능을 보여줌, 이전의 평가에서는 UnifiedGenotyper가 SNP에서 좋은 성능을  보였고, HaplotypeCaller가 indel에서 좋은 성능을 보였으나, GATK 2.7에서 SNP calling 부분을 해결함. GATK pipeline을 사용하는 경우 UnifiedGenotyepr의 indel 부분의 realigning에 성능이 떨어지므로 HaplotypeCaller를 사용할 것을 권장
  • Ensemble Calling approache: SNP와 INDEL에 있어서 가장 좋은 방법. homeozygote/heterozygote call에 대한 약간의 뒤쳐짐의 영향으로 HaplotypeCaller와 FreeBayes 모두 해당하는 것으로 이것을 반영하게 되어 discordant(share)가 높게 나오게 된다.

이미지 출처: Blue Collar Bioinformatics

calling sensitivity와 specificity와 함께 중요한 요인은 바로 분석에 걸리는 시간이다. 대략적으로 whole genome sequencing 데이터의 경우 HaplotypeCaller가 UnifiedGenotyper에 비해 7배 정도 느리며, FreeBayes보대 2배 느리다. Ensembl calling은 3가지 caller를 모두 사용하애 하기 때문에 시간이 많이 소요 된다. 그러나 이 시간은 compute infrastructure, coverage 등에 매우 변동적이다. 

Machine Learining(SVM Ensemble)을 이용한 variants calling 성능 향상

앞서서 여러가지 varinat calling method들이 서로 concordnat가 떨어지는 것을 확인했다. 그렇다면 다수의 caller를 통해 발굴된 variants를 모두 사용한다면 false positive가 많게되고, 그렇다고 concordant한 것들만 사용한다면 진짜 real variants를 놓칠 수 있게된다. 그렇다면 이시점에서 우리는 ML(Machine Learining)을 이용하는 방법을 고려해 볼 수 있다. 

잠깐 아래의 표와 벤다이어그램은 SevenBridge Genomics에서 발표한 내용으로 2가지 서로 다른 pipeline을 이용하여 Ion PGM  exome 데이터를 비교한 내용이다. 결과적으로 BWA-SW+GATK 조합의 pipeline의 경우 novel snp에 대한 Ti/Tv 비율이 1.05로 TMAP+IVH 대비 많은 true positive를 call하였다는 것을 알 수 있지만, TAMP이 놓친 ture variant를 찾았다고도 볼 수 있다. 따라서 무조건적으로 하나의 pipeline보다는 여러 pipeline(caller)를 이용하여 최대한의 성능을 낼 수 있는 방안이 필요한 것이다.


그렇다면 무조건적인 둘 혹은 그이상의 caller의 결과를 무조건 합한 결과를 사용하는 것이 아니라, 다음의 순서에 따라 ML을 적용한다.
  • training을 위해 모든 caller들이 공통적으로 찾은 variants를 true positive set(pass)으로 지정하고 하나의 caller만이 찾은 variants를 false positive(fail)로 지정한다. 
  • 위의 traing을 통해 만들어진 모델은 variant type, zygosity, reiongal sequence complexity(?)의 classifier(분류기)를 이용하여 분류를 수행한다. (ensemble) 
  • 최종적으로 분류를 통과한 variants를 사용한다.
그럼 이렇게 생성한 모델은 잘 만들어진것인지?에 대한 질문의 답은 위에서 언급했듯이 ensemble을 적용한 것을 referece call(정답을 알고있는)에 대해서 검증한 것이 있으므로 pass.

지금까지 표준 variant와 이를 통한 variant caller 비교에 대해서 알아보았다. 그리고 더 나아가 이러한 caller들을 ensemble하는 방법에 대해서도 알아 보았다. NGS 기술 발전과 보급으로 지속적으로 새로운 caller들이 저마다의 알고리즘으로 무장하고 계속 생겨나고 있다. 어느 하나의 caller를 사용하는것 보다는 이들을 모두 활용하는 방법을 통해 즉, 여러 모델의 출력을 결합하는 ensemble 방법을 통해 극적으로 성능을 향상할 수 있다. NGS variant를 검출하는데에 있어서 이러한 consensus 방법은 낮은 false positve rate에 의존하는 rare variant, pathway 분석, de-novo filtering에 유용하다.

출처: Blue Collar Bioinformatics, An automated ensemble method for combining and evaluating genomic variants from multiple callers 

지지벡터머신(SVM, support-vector machine)

분류하고자하는 범주의 평균을 사용해서 계산된 선형분류기와는 달리 각 범주와 가능한 멀리 떨어진 선(maximum-margin hyperplane, 최대허용 초평면)을 이용한다. 가장자리 점들을 구분선의 위치를 결정하는데에 사용되기 때무네 다른 데이터들을 제거하더라도 선은 동일한 장소에 존재하게 되며 이 선의 근처에 위치한 점들을 support-vector(지지벡터)라고 한다. support-vector를 찾고 이를 이용해서 구분선을 찾는 알고리즘이 바로 SVM이다.

SVM은 고차원 데이터 셋트에서 잘 작동하기 떄문에 데이터 중심의 과학 문제나 복잡한 데이터 셋트를 다루는 문제에서 자주 적용된다.
  • 얼굴표현분류
  • 국방 데이터셋트를 사용한 침입탐지
  • 단백질 서열 내의 단백질 구조 예측
  • 필기체 인식
  • 질병 예측 및 CCDS(Clinical Decision System)

Soft-maring classifier

SVM은 기본적으로 binary classifier(분류기)로 일반적으로 SVM이라고 한다면 몇개의 에러는 용인하고 줄을 긋게 된다.

Kernel

직선만 그릴 수 있기 때문에 직선으로 나뉘지 않는다면 데이터를 변형해서 둘 사이를 떨어뜨려 놓은 후 선을 긋는것으로 kernel은 떨어 뜨려 놓는, 비슷한 것끼지 뭉치게 하는 방법을 지정하는 것으로 kernel중 RBF(radial basis function) 방사함수, linear kernel, polynomical kernel이 존재한다. 즉 선형 분류로 불가능한 데이터에 대한 처리를 가능하게 한다. 계산시간이 증가하지만, 선형 분류기의 정확도를 높일 수 있다.

Radial Basis Kernel 왼쪽과 같이 선형으로 분류가 어려운 데이터를 오른쪽과 같이 커널을 적용

Python with libsvm

svm_train을 이용하여 svm_problem을 훈련 시킨후 [1,1,1]이라는 데이터를 해당 모델(m)을 이용하여 분류한다. 예측 결과 [1.0]로 분류되었으며 100%의 accuracy를 보인다.  [-1,1,1]을 예측하면 [1]이 아니며(accuracy 0%) [-1.0]으로 예측한다.
>>> from svmutil import *
>>> prob = svm_problem([1,-1],[[1,0,1],[-1,0,-1]])
>>> m = svm_train(prob, '-t 0 -c 10')
*
optimization finished, #iter = 1
nu = 0.025000
obj = -0.250000, rho = 0.000000
nSV = 2, nBSV = 0
Total nSV = 2
>>> predicted_labels, _, _ = svm_predict([1],[[1,1,1]],m)
Accuracy = 100% (1/1) (classification)
>>> print "Predicted value: " + str(predicted_labels[0])
Predicted value: 1.0
>>> predicted_labels, _, _ = svm_predict([1],[[-1,1,1]],m)
Accuracy = 0% (0/1) (classification)
>>> print "Predicted value: " + str(predicted_labels[0])
Predicted value: -1.

참고1: Ensemble

Machine Learning의 분류를 통해 여러개의 분류기(classifier)를 생성하고 이것들을 결함하여 학습하는 방법으로 다양한 분류기의 예측 결과를 결합하여 단일 분류기보다 신뢰성이 높은 예측값을 얻는 것이 목표이다. Ensemble은  SVM(Support Vector Machine)과 같은 한개의 학습 알고리즘을 사용하는 방법에 대해 주로 이루어졌으며, 학습 알고리즘을 조작하여 다양한 분류기를 생성한 후 다수결(Majority Voting)이나 가중치 투표(Weighted Voting)에 의해 예측값을 결합한다.
위에서는 ML을 이용하는 방법이었다면, 이건 오버랩되는 것들을 선택해서 하는 방법 또한 생겨나고 있다.
  • Consensus Genotyper for Exome Sequencing(CGES) Improving the Quality of Exome Variant Genotypes
    132 samples seuqenced at Hudson Alpha Institute for Biotechnology using the Numblegen Exome Capture and Illumina sequencing technology.
    In 2013, 1,184 published papers were indexed in PubMed with the key words "exome sequencing" (more than twice the number published in 2012)

참고2: 국가참조표준센터의 생명과학 분야 참조 표준

  • SNP Genotyping에 의한 1촌~6촌 Genetic Distance 유전체 생명정보
  • 한국인 차씨가계도의 유전적 거리
  • 조직적합성 판별 영역별 변이체의 상동성 기반 유전적 거리 
 
저작자 표시 비영리 동일 조건 변경 허락
신고
Software enginner of GenomeCloud. Covers bioinformatics, computational biology, and life science informatics.
Posted in : 바이오인포매틱스 at 2014.03.10 18:44
Currently 댓글이 없습니다. comments want to say something now?

티스토리 툴바