바이오인포매틱스

Split Reads

hongiiv 2015. 6. 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)

반응형