[Samtools] BAM 파일에서 Duplicated Read 찾기

안녕하세요 한주현입니다.

오늘은 bam 파일로 부터 duplicate read만 가져오는 방법에 대해 알아보겠습니다

필요한 툴, 파일

samtools

실습할 bam 파일

samtools 설치는 다음 링크를 참고해주세요.

설치 방법은 버전이 달라도 대개 비슷비슷합니다 ㅎㅎ;

http://korbillgates.tistory.com/57

NGS 분석에서는 기술특성상 필연적으로 PCR duplicate(또는 optical duplicate) 가 생길 수 밖에 없습니다.

이 부분에 대해서는 다음 포스팅에서 만나도록 하겠습니다 ㅎㅎㅎ

duplicate의 존재는 samtools로 stat을 확인해보면 나오는데요,

각자 가지고 계신 bam파일을 가지고 아래와 같은 command를 실행해봅시다.

Command

$ samtools flagstat in.bam

Result

106085813 + 0 in total (QC-passed reads + QC-failed reads)
98499 + 0 secondary
0 + 0 supplementary
8745257 + 0 duplicates
105925730 + 0 mapped (99.85% : N/A)
105987314 + 0 paired in sequencing
52993657 + 0 read1
52993657 + 0 read2
104310118 + 0 properly paired (98.42% : N/A)
105690822 + 0 with itself and mate mapped
136409 + 0 singletons (0.13% : N/A)
991440 + 0 with mate mapped to a different chr
858383 + 0 with mate mapped to a different chr (mapQ>=5)

stat 결과의 네 번째 bold 된 줄을 보시면 duplicates가 있지요?

duplicate는 picard등의 tool이 duplicate로 marking을 해줘서 저렇게 flagstat 명령어로 알 수 있습니다.

marking은 bam 파일의 flag 필드에 1024로 마킹이 됩니다.

bam 파일을 samtools view 로 열어서 보면 아래처럼 나오는데 두번째 컬럼인 FLAG 를 보시면..

Query template NAME	bitwise FLAG	Reference sequence NAME	1-base leftmost mapping POSition	MAPping Quality
ABCD:123:12345AA:1:1234:1234:1234	145	chr1	51	60
ABCD:123:12345AA:1:1234:5234:2213	1105	chr1	51	60
ABCD:123:12345AA:1:1234:6234:4241	81	chr1	51	60

이런식으로 duplicate를 FLAG 형식으로 표시해줍니다.

우리는 FLAG 정보로 duplicate를 골라낼 수 있는 것 이지요..

flag의 자세한 내용은 아래를 참고 하세요 ㅎㅎ..

https://samtools.github.io/hts-specs/SAMv1.pdf

Bit	Description
1	template having multiple segments in sequencing
2	each segment properly aligned according to the aligner
4	segment unmapped
8	next segment in the template unmapped
16	SEQ being reverse complemented
32	SEQ of the next segment in the template being reverse complemented
64	the first segment in the template
128	the last segment in the template
256	secondary alignment
512	not passing filters, such as platform/vendor quality controls
1024	PCR or optical duplicate
2048	supplementary alignment

자.. 그럼 이제 duplicated read를 찾으러 가볼까요?

Command

$ samtools view -f 1024 in.bam

Result

ABCD:123:12345AA:1:1234:5234:2213  1105    chr1    51      60
ABCD:123:12345AA:1:2421:4334:2241  1187    chr1    68      60
ABCD:123:12345AA:1:3131:2134:5313  1123    chr1    98      60
ABCD:123:12345AA:1:4124:3214:1143  1107    chr1    169     60

flagstat 의 결과와도 비교해볼까요?

Command

$ wc -l in.onlydup.sam

Result

8745257 in.onlydup.sam

samtools를 가지고 duplicate read를 골라내는 방법에 대해 알아보았습니다.

그럼 다음 시간에 만나요~~