samtools flagstat로 출력되는 숫자를 재현하기

samtools flagstat는 SAM/BAM 파일에 수록된 read alignment 정보를 출력해 주는 매우 유용한 프로그램이다. 그러나 출력되는 정보는 기본적으로 read가 아닌 'alignment'라서 직관적으로 이해하기 어려울 때가 많다. 특히 하나의 read가 reference 서열에 존재하는 반복 서열에 의해 여러 곳에 매핑될 경우 계산이 잘 맞지 않는 것처럼 느껴진다. 이것을 명확하게 이해하기 위하여 samtools view -f (또는 -F) INT 명령의 flag를 잘 조합하여 samtools flagstat로 표현되는 수치와 같은 값이 나오게 만들어 보았다.

솔직하게 고백하자면 나는 거의 항상 CLC Genomics Workbench에서 박테리아 게놈의 NGS data를 다루기 때문에(de novo assembly 위주의 업무) SAM/BAM 파일을 직접 다룰 일은 없다고 해도 과언이 아니다. 하지만 데이터 조작의 기본 원리를 이해하거나 사용 프로그램을 접하기 어려운 사람에게 설명을 하려면 CLC Genomics Workbench만을 고집할 수는 없다. 그래서 IGV(Integrated Genomics Viewer)도 가끔은 실행해 봐야 한다.

미생물 유전체 서열로부터 64780 x 2(129560 reads, 2x150 nt)개의 가짜 fastq file을 만들었다. 사용한 시뮬레이터는 ART MountainRainier이다. 레퍼런스 서열은 read를 만들어낸 원본 read와 가까운 유전체를 대략 같은 간격의 약 500개의 조각(크기는 1 kb)으로 나눈 것이었다. 매핑에는 bowtie2를 사용하였다.

$ samtools view -H mapping.sam # 헤더 정보만 보여줌
$ samtools view -b mapping.sam > mapping.bam
$ samtools view -b mapping.sam -o mapping.bam # 같은 결과
$ samtools flagstat mapping.bam 
[1] 131900 + 0 in total (QC-passed reads + QC-failed reads)
[2] 2340 + 0 secondary
[3] 0 + 0 supplementary
[4] 0 + 0 duplicates
[5] 21390 + 0 mapped (16.22% : N/A)
[6] 129560 + 0 paired in sequencing
[7] 64780 + 0 read1
[8] 64780 + 0 read2
[9] 15554 + 0 properly paired (12.01% : N/A)
[10] 15710 + 0 with itself and mate mapped
[11] 3340 + 0 singletons (2.58% : N/A)
[12] 122 + 0 with mate mapped to a different chr
[13] 105 + 0 with mate mapped to a different chr (mapQ >=5)

그러면 각 줄번호에 해당하는 숫자를 재현해 보자. 0만 표시된 것은 제외하였다. alignment의 수가 필요한지(samtools view -c), 혹은 이로부터 read name에 해당하는 것을 뽑아낼 것인지를 명확히 해야 된다.

[1]은 read의 수가 아니라 alignment의 수이다. 다음과 같이 하면 된다.

$ samtools view -c mapping.bam
131900

[2]는 primary와 secondary mapping을 이해해야 한다. 256을 16진수로 표현하면 0x100이다. 따라서 -f 0x100이라고 써도 된다.

$ samtools view -c -f 256 mapping.bam
2340

수치가 나오지 않은 [3,4]는 생략하자. [5]는 mapped read의 alignment 수를 의미하므로, SAM/BAM flag의 0x04(unmap)의 반대, 즉 -F 0X04를 쓰면 된다.

$ samtools view -c -F 0x04 mapping.bam
21390

만약 properly paired read에 대한 것만을 추출하고 싶다면 -f 0x03(1 for paired, 2 for proper_pair)으로 설정해야 한다.

[6]은 매핑 결과와는 관계없이 최초 read에 들어있던 read의 수를 의미한다. 이는 이미 알고 있는 숫자이다. 그러나 BAM 파일에서 굳이 구하고 싶다면 다음 계산에서 나온 수치에 2를 곱하면 된다. BAM file에서는 simulated read name 끝부분의 /1 또는 /2가 제거된 상태라서, uniq 명령을 거치고 나면 그 수가 절반으로 줄기 때문이다. sort는 꼭 하지 않아도 된다.

$ samtools view mapping.bam | cut -f1 | sort | uniq | wc -l
64780

[7]과 [8]은 fwd 및 rev read file에 들어 있는 read의 수이므로 서로 동일하며, 이미 알고 있는 값이므로 넘어간다. 전체 read 수의 절반임은 당연하다.

[9]는 조금 어렵다. Properly paired를 의미하는 flag은 0x02인데, 이는 반드시 0x01(paired)과 같이 쓰인다. 따라서 -f 0x03을 주어서 read의 수를 계산한 다음, pair를 만들기 위해 2를 곱하면 된다. samtools view -c -f 0x03 mapping.bam이라고 명령을 내리고 싶은 강렬한 욕구가 느껴지지만, 그렇게 하면 안된다. 왜냐하면 proper pair를 이룬 상태로 reference 내의 여러 곳에 결합할 수도 있기 때문이다. 실제로 이렇게 계산하면 15554가 아니라 16760이 나온다.

$ samtools view -f 0x03 mapping.bam | cut -f1 | sort | uniq | wc -l
7777

[10]은 read 쌍이 둘 다 매핑이 된 것을 의미한다. Read 쌍의 거리와 방향이 적정한지는 중요하지 않다. 이것은 어떻게 flag를 정해야 할까? read도 unmapped가 아니고(-F 0x04), 그의 짝도 unmapped가 아니여(-F 0x08)! 따라서 -F 0x12를 하여 계산되는 read의 수에 2를 곱한다.

$ samtools view -F 12 mapping.bam | cut -f1 | sort | uniq | wc -l
7855

다음으로 [11]은 꽤 까다롭다. sinlgeton, read pair의 한쪽만 mapping이 된 것을 산출하는 것이다. 여기에서는 -f 0x04 -F 0x08을 하거나 -f 0x08 -F 0x04를 하면 된다. 두 명령어를 따로 실행해서 합쳐야 할 것 같은 생각이 들지만, read name을 각각 추출하여 비교해 보면 똑같다. 이것은 read pair를 산출하는 것이 아니므로 곱하기 2를 하면 안된다.

$ samtools view -f 0x04 -F 0x08 mapping.bam | cut -f1 | sort | uniq | wc -l
3340

[12,13]은 아직 정확한 flag 조합을 찾지 못하였다. SAM/BAM flag의 의미를 파악하려면 매SAM/BAM speficiation 파일이나 samtools flags 명령어, 또는 SAM flag 디코더를 적절히 활용해 보라.

Paired read, 즉 [5]와 singleton[11]은 후속 분석에서 유용할 때가 많다. 특히 11번에 해당하는 read의 pair를 추출해 놓으면 read pair에 의한 연결 부위를 찾는데 활용할 수 있다. 이상에서 살펴본 바와 같이 SAM/BAM flag을 자유롭게 쓸 수 있게 되면, samtools bam2fq 명령을 이용하여 BAM 파일에서 조건에 맞는 read를 직접 추출할 수도 있다. pair가 깨진 상태라면 BBmap의 repair.sh 및 관련 명령어(Repair Guide)를 적절히 활용하면 매우 편리하다.