2019년 8월 7일 수요일

SAM flag에 대해서 좀 더 알아보자 - mapping 여부에 따른 read 추출하기

SAM/BAM 파일로부터 properly paired (mapped) read를 꺼내려면 -f 0x03 flag를 쓰면 충분하다. 더욱 간단하게는 두번째 비트가 1이면 된다. -f 0x02의 의미가 proper pair이기 때문이다. 그러나 paired read가 쓰인 경우 항상 0x01과 따라다니게 되므로 -f 0x03이 더욱 바람직하다.

이토록 간단한 flag이지만, 실제 매핑이 이루어진 read의 입장에서 4가지의 경우로 세분한다면 이해의 폭을 더욱 넗힐 수 있다. Paired fastq file(I & II)을 이용하여 매핑을 하였을 때 하나의 참조서열 내부에서 정상적인 pair를 이루는 모든 경우는 다음이 전부이다. 화살표는 참조서열 위에 매핑일 이루어진 read를, 별표(*)는 fastq 파일에 수록된 상태의 서열을 reverse complementary로 전환해야 함을 의미한다. 내가 알기로 first 및 second는 _1 및 _2 파일을 의미한다. pair를 이루는 시작(left)과 끝(right) read를 의미하는 것이 아니다.

First           Second        Second           First
-------▶......◀-------     -------▶......◀-------
from I         from II*       from II         from I*
  [A]            [B]            [C]            [D] 

(아, 도대체 html로는 부등호 기호를 써서 화살표를 표현하기가 어려우니...)

A-B의 네 가지 read는 단 한번씩만 mapping되어서 SAM 파일 내에서 하나의 alignment record로 존재한다고 가정하자. 그러면 각각의 alignment는 어떤 flag을 갖게 될 것인가? flag의 설명은 samtools flag 명령으로 표시되는 것과 SAM flag decorder에서 쓰이는 용어가 조금 다르지만 편의상 후자의 것을 쓰기로 하자.

기본적으로 여기에 보인 네 가지 alignment는 전부  웹사이트기본적으로 read paired(0x01)와 read mapped in proper pair(0x02)를 만족하므로 0x03을 기본값으로 깔고 있다. 다음으로 생각할 것은 이 read가 pair 중 첫번째 혹은 두번째 read인가, 그리고 이 read와 mate 중 어느 것이 참조서열에 대하여 reverse로 존재하는가의 여부이다. 이는 다음과 같이 4가지의 flag으로 표현된다.

  • 0x10: read reverse strand
  • 0x20: mate reverse strand
  • 0x40: first in pair
  • 0x80: second in par
A는 mate reverse strand(0x20), first in pair(0x40)이다. 따라서 0x03 + 0x20 + 0x40 = 0x63 = 99
B는 read reverse strand(0x10), second in pair(0x80)이다. 따라서 0x03 + 0x10 + 0x80 = 0x93 = 147
C는 mate reverse strand(0x20), second in pair(0x80)이다. 따라서 0x03 + 0x20 + 0x80 = 0xa3 = 163
D는 read reverse strand(0x10), first in pair(0x40)이다. 따라서 0x03 + 0x10 + 0x40 = 0x53 = 83

99/147, 그리고 163/83은 NGS read alignment를 다루는 사람이라면 잊어서는 안되는 숫자가 되겠다. samtools view -f 99로 추출한 read는 -f 147로 얻은 것과 mate 관계이고, -f 163과 -f 83 역시 같은 관계이다. 이렇게 얻은 모든 파일을 합치면 -f 0x03으로 추출한 read 세트와 동일한 것이 된다. 혹시 내가 전개한 이론이 맞지 않다고 생각하는 독자가 계시다면 언제든지 연락 주시길!

Unmapped read에 대해서는?

read unmapped는 0x04(0x4, 또는 십진수로 4라고만 써도 되는데 버릇이 되어서 계속 두 자리 16진수로 표현하고 있다), mate unmapped는 0x08이므로 이를 적절히 조합하면 unmapped pair와 어느 한쪽만 mapping이 된 pair를 찾을 수 있다. 내가 생각하는 가장 간단한 flag 값은 이러하다. 
  • both unmapped: -f 12 (0x04 + 0x08, 16진수로는 0x0c)
  • read mapped, mate unmapped: -f 0x08 -F 0x04
  • read unmapped, mate mapped: -f 0x04 -F 0x08
그런데 Novocraft의 웹사이트(Extracting unmapped reads from a BAM file produced by NovoAlign)에서는 좀 더 복잡한 flag 값을 제시한다.
  • both unmapped: -f 12 -F 256
  • read mapped, mate unmapped: -f 8 -F 260
  • read unmapped, mated mapped: -f 4 -F 264
숫자 단위가 큰 flag은 무엇을 의미하는지 알아보자.
  • 256(0x100): not primary alignment
  • 260: 256 + 0x04(read unmapped)
  • 264: 256 + 0x08(mate unmapped)
이것을 -F로 적용하니 직관적으로 이해를 하기가 좀 어렵다. 테스트를 해 보면 내가 제시한 flag로 하는 것과 결과에 차이가 없다. 이로써 이틀에 걸친 SAM flag 공부하기를 마치고자 한다. SAM and BAM filtering one-liner도 참고하기에 좋다.



댓글 없음: