Fork me on GitHub 단맛만좋아요 ::

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?

티스토리 툴바