[samtools] BAM 파일에서 특정 chromosome 영역 추출하기

 

 

 

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

 

오늘은 samtools로 BAM파일에서 원하는 chromosome 또는 영역을 추출하는 방법

에 대하여 알아보겠습니다.

 

  samtools

 

samtools 로 BAM파일에서 원하는 chromosome 또는 영역을 추출하는 방법 - 준비물

준비물

1) samtools

samtools 설치가 되어 있지 않다면 다음 링크들을 참고해 주세요

  - samtools 공식 다운로드 페이지

http://www.htslib.org/download/

  - samtools 설치 방법

http://korbillgates.tistory.com/57

2) 실습 bam 파일

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam

Size: 92MB

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam.bai

Size: 4.9MB

 

 

samtools 로 BAM파일에서 원하는 chromosome 또는 영역을 추출하는 방법 - 실행

1) samtools 로 BAM파일에서 원하는 chromosome 추출하기

samtools 커맨드 중 view 커맨드를 다음과 같이 사용합니다.

samtools view sample1.bam chr21 > sample1.chr21.bam

단순히 chr 이름만 적어주면 딱 결과가 std out 으로 나오게 됩니다. 

이를 파일로 받으면 저장이 됩니다 ㅎㅎ;;

너무 간단하죠?

그러면 여러개의 chromosome 도 가능할까요?

네 됩니다 ㅎㅎ

samtools view sample1.bam chr21 chr22 > sample1.chr21_chr22.bam

정말 너무 간단하게 뚝딱 만들어 냈습니다 ㅎㅎ

2) samtools 로 BAM파일에서 원하는 영역 추출하기 - 직접 영역 지정

BAM파일에서 원하는 영역만을 추출하고 싶으시면 chr2:10100-10110 과 같이 영역을 써주시면 됩니다.

samtools view sample1.bam chr2:10100-10110 | head

결과)

범위에 걸쳐있는 read에 해당하는 영역이 나오게 됩니다 ㅎㅎ

SOLEXA-1GA-2_2_FC20EMB:5:127:459:202    0    chr2    10083    25    36M    *    0    0    CCACACCCACACACAATACCCACACCATACCCACAC    hhhhh[fh^hZhScHMc`USVGRHNZNIEDIHCNFL    NM:i:1    X1:i:1    MD:Z:14A1T10T8
SOLEXA-1GA-2_2_FC20EMB:5:25:447:574    0    chr2    10088    25    36M    *    0    0    CCCACACACGACACACACACGACACCCACACACCAC    RPZJMIOBLMGJFNHGFJAEBC>?G?D>A?=?=@>>    NM:i:1    X1:i:1    MD:Z:9G4A5G15
SOLEXA-1GA-2_2_FC20EMB:5:110:897:676    16    chr2    10088    25    36M    *    0    0    CCAACCCACCACACCCACACCACACACACACACCAC    @@??AFFBGFEACMKOHWMZh\b^UR^bh_h^hh_h    NM:i:1    X1:i:1    MD:Z:10A19C2A2

3) samtools 로 BAM파일에서 원하는 영역 추출하기 - BED 파일 사용

텍스트 에디터로 chf2    10100    10110을 bed 형식으로 저장합시다

$ cat range.bed
chr2    10100    10110

samtools view sample1.bam -L range.bed

SOLEXA-1GA-2_2_FC20EMB:5:127:459:202    0    chr2    10083    25    36M    *    0    0    CCACACCCACACACAATACCCACACCATACCCACAC    hhhhh[fh^hZhScHMc`USVGRHNZNIEDIHCNFL    NM:i:1    X1:i:1    MD:Z:14A1T10T8
SOLEXA-1GA-2_2_FC20EMB:5:25:447:574    0    chr2    10088    25    36M    *    0    0    CCCACACACGACACACACACGACACCCACACACCAC    RPZJMIOBLMGJFNHGFJAEBC>?G?D>A?=?=@>>    NM:i:1    X1:i:1    MD:Z:9G4A5G15
SOLEXA-1GA-2_2_FC20EMB:5:110:897:676    16    chr2    10088    25    36M    *    0    0    CCAACCCACCACACCCACACCACACACACACACCAC    @@??AFFBGFEACMKOHWMZh\b^UR^bh_h^hh_h    NM:i:1    X1:i:1    MD:Z:10A19C2A2

같은 결과를 얻을 수 있습니다 ㅎㅎ

 

Header 와 함께 bam 에서 bam 으로 영역 추출하기

samtools 에 -b 옵션은 bam으로 출력, -h 옵션은 header와 내용 포함이란 옵션입니다 ㅎㅎ..

이 두개를 같이 써서 -bh 로 쓰고 원하는 영역을 넣어 주고 > 기호로 standard out을 받습니다

$ samtools view -bh Sample.bam "chr1:123456-789012" > Sample.extracted.bam

생각보다 간단하지요? ㅋㅋ

오늘은 samtools를 사용하여 BAM 파일에서 부터 특정 chromosome 또는 영역을 추출하는

방법에 대하여 알아봤습니다.

여러분들께 도움이 되셨음 좋겠습니다 ㅎㅎ.

그럼 다음 시간에 만나요!

참고자료 : http://www.htslib.org/doc/samtools.html