안녕하세요 한헌종입니다!
오늘은 sequencing data 분석에 굉장히 많이 쓰이는 samtools 라는 툴을 사용하는 예제를 적어보고자 합니다.
samtools 는 BAM, SAM 형태의 파일을 읽고, 쓰고, 조작할 수 있게 해줍니다.
영어로 된 설명은 여기서 찾아볼 수 있습니다: Samtools 사이트 링크
내용이 매우 길어서, 아직 미완성이라 계속 보충해 나갈 예정입니다.

부연설명을 하자면, Samtools 는 BCFtools, HTSLib 와 관계가 있습니다.
HTSLib 은 high-throughput sequencing 데이터를 읽고 쓰게 해주는 라이브러리 이구요,
이를 Samtools, BCFtools 내부에서 사용하고 있습니다.
Samtools 는 BAM/SAM 을 다루는 툴이고, BCFtools는 BCF/VCF 파일을 다루는 툴입니다.
(BCF는 VCF의 binary 형태입니다. BCF-VCF 관계는 BAM-SAM 관계와 비슷하죠.)
Samtools, BCFtools 는 서로 다른 repository 로 각각 설치 가능하고, 이 때 HTSLib이 자동으로 같이 설치됩니다.
물론 HTSLib은 따로 설치할 수도 있습니다.

samtools 설치부터 알려주세요

samtools 를 설치하는 방법은 다양합니다.

  1. source code 를 사용해서 설치:
    다음 링크에서 samtools source code 를 받을 수 있습니다: 다운로드 링크
    다운로드 한 뒤 압축을 풀고, 다음 명령어로 설치 가능합니다. 
    $ cd samtools-1.x    
$ ./configure --prefix=/where/to/install
$ make
$ make install
  2. github 를 통해서 설치:
    다음 링크에 samtool github repository 가 있습니다. samtools github
    (samtools github에는 이것 말고도 bcftools, htslib 등 다양한 repository 가 있군요.)
    clone으로 저장소를 불러온 다음 설치하면 되겠습니다. 
    $ git clone https://github.com/samtools/samtools.git
$ cd samtools
$ ./configure
$ make
$ make install
  3. conda 를 사용해 설치
    만약 anaconda 를 사용하고 계시다면, 더 쉽게 설치가 가능합니다.
    다음 명령어를 사용해보세요 (링크 참조) 
    $ conda install -c bioconda samtools

다양한 samtools 명령어들

설치가 끝나셨다면 samtools 명령어를 실행해볼까요?
터미널에서 samtools라고 작성하면 다음과 같이 나옵니다:

$ samtools

Program: samtools (Tools for alignments in the SAM format)
Version: 1.11 (using htslib 1.11)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags
     markdup        mark duplicates
     ampliconclip   clip oligos from the end of reads

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region
     coverage       alignment depth and percent coverage
     depth          compute the depth
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)
     ampliconstats  generate amplicon specific stats

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM

명령어가 굉장히 많군요!
이 명령어들을 하나씩 살펴보기로 하죠.
그 전에! BAM/SAM 파일 형식을 잘 숙지하고 계셔야 쉽게 이해가 될 것입니다. (참고링크)

*참고로 데이터는 1000 genome project 에서 NA12878 샘플의 chromosome 20번 데이터를 골라봤습니다. (링크)

Indexing: fasta/fastq/bam 파일의 index 를 만들고 싶을 때

  1. dict
    fasta 파일에서 sequence dictionary 를 만들 때 쓰입니다.
    이를 통해 SAM 형식의 header 가 나오게 되는데, 여기에는 fasta 파일의 contig 수, 길이 등이 나오게 됩니다.
    다음 명령어를 통해 진행합니다: 
    $ samtools dict chr20.fa
@HD	VN:1.0	SO:unsorted
@SQ	SN:20	LN:63025520	M5:0dec9660ec1efaaf33281c0d5ea2560f	UR:file:/path/to/fasta/chr20.fa
  2. faidx
    faidx는 fasta 파일을 indexing 할 때 씁니다.
    이 indexing 과정은, read data (fastq) 를 reference 에 빠르게 mapping 할 때 쓰입니다.
    reference fasta 파일을 활용해서 무언가 작업을 하고 싶을 땐, fasta index를 먼저 만들어놓아야 합니다.
    다음 명령어를 실행하면 fasta index file 인 reference.fa.fai 파일이 만들어 집니다. 
    $ samtools faidx reference.fa
  3. fqidx
    fqidx 는 fastq 로 되어있는 reference sequence 의 index를 만들 때 쓰입니다.
    fastq 형태로 된 reference sequence는 전 아직 본적 없지만요…
    이 명령어의 경우, short reads 가 많은 fastq를 사용하면 거의 원본만큼의 크기를 가진 큰 idx 파일이 나오고, 이를 통해 연산할 경우 메모리를 매우 많이 잡아먹는다고 합니다.
    그래서 input으로는 적은 수의 sequence 를 가진 fastq 파일을 쓰는 것이 좋습니다. 
    $ samtools fqidx reference.fq
  4. index
    여기서의 indexing 은 bam file indexing 을 말합니다.
    이 과정은 bam 파일에서 특정 영역의 alignment 정보를 빨리 불러올 수 있게 만드는 bai 파일을 만들어 줍니다.
    bam 파일에서 어떤 작업을 하고싶을 땐 이 index 파일이 있어야 합니다.
    사용법은 다음과 같습니다: 
    $ samtools index example.bam

Editing: bam 파일을 수정하고 싶을 때

  1. calmd
    BAM 파일의 MD tag 와 NM tag 를 다시 계산하고자 할 때 사용하는 명령어입니다.
    MD tag 는 reference sequence 와 얼마나 일치하는지를 표현하는 방법입니다. 링크 참조
    NM tag 는 edit distance to reference 라는건데요, 쉽게 말해서 number of mismatches 입니다.
    사용법은 다음과 같습니다. 
    $ samtools calmd -b example.bam reference.fa > result.bam

    결과물이 standard output으로 나오므로, output file 을 지정해주는 것이 좋습니다.
    또한, default output 형식은 sam 파일이므로, -b 옵션을 설정해주면 bam 결과물을 얻을 수 있습니다.

  2. fixmate
    read alignment 를 수행하는 툴은 여러가지인데요, 어떤 툴은 가끔 FLAG (2nd column) 값을 이상하게 만드는 경우가 있습니다.
    paired-end reads 분석을 위해서는 fragment 를 이루는 read pair 정보를 정확히 가지고 있는것이 중요합니다.
    때문에 fixmate 기능으로 FLAG 값을 수정해주는 작업이 필요할 때가 있습니다.
    fixmate 를 사용하려면 먼저 read 들이 genomic coordinate 가 아닌, 이름순으로 나열되어 있어야 합니다.
    다음과 같이 실행하시면 됩니다. 
    $ samtools sort -n example.bam | samtools fixmate - example.fixmate.sam

    (‘-‘ 표시는 이전 명령어의 standard output을 현재 명령어의 input으로 받는다는 뜻입니다.)
    만약 이름으로 이미 나열되어 있다면 다음과 같이 실행하시면 됩니다.

    $ samtools fixmate -O bam example.bam example.fixmate.bam

    input 은 sam, bam 아무거나 상관이 없는데, output 형태가 bam 파일이 되려면 -O bam 이라는 옵션을 추가해줘야 합니다.

  3. reheader
    bam 파일에 header 가 없거나, header 가 이상해서 이를 수정하고 싶을 때가 있습니다.
    그런 경우에는 header sam 파일을 만들어 놓고, 다음과 같은 명령어를 수행하시면 됩니다. 
    $ samtools reheader new_header.sam example.bam > result.sam
  4. addreplacerg
    이 명령어는 read group 명을 더하거나 바꾸고 싶을 때 사용합니다.
    bam header 에 @RG 로 표현된 부분이 바로 read group 에 대한 정보이며, 각 read 들은 이런 read group ID 중 하나를 가지고 있습니다.
    read group 에 대한 자세한 설명은 다음 링크에서 보실 수 있습니다. 링크
    -r 옵션으로 @RG 전체 string을 지정해줄 수도 있으며, -R 옵션으로 이미 있는 RG ID 를 지정해줄 수도 있습니다.
    -m 옵션으로 덮어쓸지 추가할지를 결정할 수 있는데요, -m orphan_only 옵션으로는 추가, -m overwrite_all 으로는 덮어쓰기가 가능합니다. 
    $ samtools addreplacerg -r '@RG\tID:NG01G\tLB:NG01G\tPL:ILLUMINA\tSM:NG01G' -m orphan_only example.bam -o example.addrg.bam

    위 명령어를 통해 ‘@RG\tID:NG01G\tLB:NG01G\tPL:ILLUMINA\tSM:NG01G’ 라는 read group 이 bam file에 추가됩니다.

  5. markdup
    markdup 은 duplicated alignments 를 표시하기 위한 기능입니다.
    먼저 input bam file 은 coordinate 순서대로 나열되어 있어야 합니다.
    또한 위에서 언급한 fixmate 를 -m option 추가한 상태로 돌린 뒤에 가능합니다.
    그래서 1) name sort 진행 후 2) fixmate -m 진행하고 3) coordinate sort 를 진행한 뒤에 4) markdup 기능을 사용할 수 있습니다.
    다음 명령어로 진행합니다: (‘-‘ 표시는 이전 명령어의 standard output을 현재 명령어의 input으로 받는다는 뜻입니다.) 
    $ samtools sort -n example.bam | samtools fixmate -m - example.fixmate.bam
$ samtools sort example.fixmate.bam | samtools markdup - example.markdup.bam

    이 때, duplicated alignments를 없애고 싶다면 markdup 에 -r 옵션을 사용하면 됩니다.

File operations: bam 파일을 합치고 나누고, 나열하는 방법

  1. collate
    collate read 들을 read group 별로 나열해주는 작업입니다.
    이 작업은 ‘read groups’ 를 알고나면 이해하기 쉽습니다. (관련 링크)
    요약하면, read group 이란 각 single lane 마다 부여되는 이름입니다.
    최근 시퀀싱 기술은 여러 sample libraries 들을 pooling 해서 한번에 분석하는 경우가 많은데요,
    이런 경우 sample 들을 여러 lane 에 나누어서 분석하는 것을 바로 multiplexed sequencing 이라고 합니다. 더 빠르고, 효율적인 방법이죠.
    결과로 나온 reads 에는 어떤 lane에서 분석된 것인지를 이름으로 구분해주는데요, 우리는 이 정보를 활용해서 각 lane 에서 오는 영향(bias 등)을 관찰하고 활용할 수 있습니다.
    collate 는 read 들을 이런 read group 별로, 즉 sequencing lane 별로 묶어서 나열해주는 작업을 진행합니다.
    다음 명령어로 진행합니다. 
    $ samtools collate example.bam --output-fmt sam -o example.collate.sam

    output 형식은 –output-fmt 옵션으로 지정해주시고, output 파일은 -o 옵션으로 지정해주시면 됩니다.

  2. cat
    cat 은 concatenate 의 약자입니다. 이름대로 여러개의 bam 파일을 합칠 때 사용합니다.
    합치는 방식은 정말 간단하게 각 파일을 이전 파일의 맨 뒤에 붙이는 방식으로 진행합니다.
    단, 조건이 있습니다. 먼저 파일들의 형식이 모두 같아야 하고, 파일들은 모두 똑같은 sequence dictionary 를 가져야 합니다.
    사용법은 간단합니다. -o 옵션으로 결과를 저장합니다. 
    $ samtools cat example_1.bam example_2.bam -o example_cat.bam

    합칠 파일들을 특정 텍스트파일에 적어놓고 사용할 수도 있습니다. 한 줄에 파일이름 하나씩 적어놓은 file_list.txt 를 만들고, -b 옵션으로 다음과 같이 쓰시면 됩니다:

    $ samtools cat -b file_list.txt -o example_cat.bam

    header 는 default 로 첫번째 파일의 header 를 사용하나, 이를 -h 옵션으로 다른 파일의 헤더를 사용할 수도 있습니다.

  3. merge
    merge 역시 여러개의 bam 파일을 합칠 때 사용합니다.
    cat 명령어와 같다고 생각하실 수 있으나, 한 가지 차이점이 있습니다.
    바로 ‘sorted output’을 결과물로 준다는 것인데요.
    그러기 위해서는 input 파일들이 모두 같은 알고리즘으로 나열되어 있어야 합니다.
    즉, 한 파일을 samtools sort -n 옵션으로 나열했다면, merge 하려는 다른 파일들 역시 -n 옵션으로 나열되어 있어야 합니다.
    사용법은 다음과 같습니다: 
    $ samtools merge example_1.bam example_2.bam -o example_cat.bam

    여기서도 -h 옵션으로 헤더를 지정해줄 수 있는데요, 지정하지 않는 경우엔 input header 들이 모두 합쳐집니다.

  4. mpileup
    bam 파일을 사용해 pileup 형식의 파일을 만들어내는 명령어입니다.
    pileup 형식이란? SNP/indel 등의 variant 정보를 한눈에 확인할 수 있는 파일 형식입니다.
    자세한 설명은 다음 링크에서 확인하실 수 있어요! (링크)
    간단한 사용법은 다음과 같습니다. -o 로 output 파일을 만들어 줍니다. 
    $ samtools mpileup example.bam -o example.pileup

    옵션을 살펴보니 몇가지 알아두면 좋을 항목이 있었습니다:
    -f [fasta 파일]: indexing 되어있는 reference fasta 파일을 사용할 수 있습니다.
    -l [bed 파일]: 지정해준 bed 파일에 나타난 영역에 대해서만 mpileup 을 진행합니다.
    -b [bam 파일 목록이 적힌 텍스트 파일]: 여러 bam 파일에 대해 작업할 때, 파일 목록을 만들어서 줄 수 있습니다.

  5. sort
    bam 파일에 있는 read 들을 나열할 때 사용합니다.
    간단한 사용법은 다음과 같습니다. -o 옵션으로 ouput 파일을 지정해 줍니다. 
    $ samtools sort example.bam -o example_sorted.bam

    sort 에는 여러 유용한 옵션들이 있으니 한번 참고해볼까요?
    -n: read 를 나열할 때 read name 순으로 나열합니다.
    -t: read 를 나열할 때 TAG 숫자 순으로 나열합니다.
    -T [prefix]: sort 를 진행하면 temporary 파일들이 많이 사용되는데요, 이 때 temporary 파일의 이름을 지정해서 prefix.nnnn.bam 으로 만들도록 합니다.
    -@ [thread 수]: 나열할 때 사용할 thread 를 지정해줍니다.
    -m [memory]: 나열할 때 사용할 thread 별 최대 메모리를 지정해줍니다.

  6. split
    bam 파일에 있는 reads 들을 read group 별 독립된 파일들로 나눠줍니다.
    이 작업은 ‘read groups’ 를 알고나면 이해하기 쉽습니다. (관련 링크)
    사용법은 다음과 같습니다. 
    $ samtools split example.bam

    결과물로는 read group 마다 example_0.bam, example_1.bam … 으로 나옵니다.

  7. quickcheck
    이름 그대로 bam 파일이 온전한지를 빠르게 체크하는 명령어입니다.
    설명에 따르면 header 를 제대로 가지고 있는지, target sequnece를 최소한 하나 가지고 있는지, 그리고 파일 끝 부분이 EOF (end-of-file) 을 제대로 가지고 있는지를 확인한다고 하는군요.
    파일 하나씩 체크할 수도 있지만, 매뉴얼에는 여러 파일을 한번에 체크하는 방법을 추천하고 있습니다.
    다음 명령어 입니다. 한번 확인해볼까요? 
    $ samtools quickcheck example.bam
$ samtools quickcheck *.bam && echo 'all ok' || echo 'fail!'
$ samtools quickcheck -qv *.bam > bad_bams.fofn && echo 'all ok'  || echo 'some files failed check, see bad_bams.fofn'

    첫 번째 명령어는 파일 하나만 체크하는 방법입니다.
    두 번째 명령어는 모든 bam 파일 확인 후, 정상이라면 ‘all ok’를 출력하고, 그렇지 않으면 ‘fail!’을 출력합니다.
    세 번째 명령어는 모든 bam 파일 확인 후, 비정상 bam 파일 목록은 bad_bams.fofn 에 저장하며, 모두 정상이라면 ‘all ok’를 출력하고, 그렇지 않으면 ‘fail!’을 출력합니다.

  8. fastq
    bam 파일을 fastq 로 바꿔주는 명령어입니다.
    다음 명령어로 사용합니다: 
    $ samtools fastq example.bam -o example.fastq

    만약 paired-end sequencing 데이터라면 read1, read2 파일을 따로 지정해주면 됩니다.
    이 때, 두 개의 fastq 파일에 같은 read 순서로 만들고 싶다면 먼저 collate 혹은 sort -n 으로 나열되어 있어야 합니다.

    $ samtools fastq example.bam -1 example_1.fastq -2 example_2.fastq -0 example_0.fastq

    이 때, read1, read2 이외의 read 들은 standard out으로 출력이 될 겁니다.
    어마어마한 수의 reads 정보가 출력될 것이기 때문에, 만약 출력하고 싶지 않다면 위 예시처럼 숫자 ‘0’ 옵션을 사용해서 -0 [파일이름] 으로 지정해주셔야 합니다.

  9. fasta
    bam 파일을 fasta 로 바꿔주는 명령어입니다.
    매뉴얼에서는 samtools fastq 와 samtools fasta 를 같이 설명하고 있군요!
    fastq 에서와 마찬가지로 read1, read2, other reads 를 구분해서 쓰시면 되겠습니다.
    다음 명령어로 사용합니다: 
    $ samtools fasta example.bam -o example.fasta
$ samtools fasta example.bam -1 example_1.fasta -2 example_2.fasta -0 example_0.fasta

Statistics: bam 파일의 통계치를 확인하고 싶을 때

  1. bedcov
    특정 영역에 대해서 read base count 를 세어주는 명령어 입니다.
    read base count 는 the sum of per base read depths, 즉 mapping 된 총 base 수를 합친 값입니다.
    영역은 bed 파일로 지정해주면 되고, 각 영역별 read base count 를 세어줍니다.
    다음 명령어로 실행합니다. 
    $ cat example.bed # 이 형태의 bed 파일을 사용한다고 하면...
20	60000	70000
20	70000	80000
20	80000	100000
$ samtools bedcov example.bed example.bam # 이 명령어로 아래와 같은 결과가 나옵니다.
20	60000	70000	51006
20	70000	80000	163986
20	80000	100000	51060
  2. coverage
    특정 영역에 대해서 read coverage 를 표시한 histogram 을 만들어주는 명령어 입니다.
    아래 명령어로 실행합니다. 
    $ samtools coverage -r 20:100000-200000 example.bam
#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
20	100000	200000	7577	87538	87.5371	5.32686	33.7	54.3
$ samtools coverage -r 20:100000-200000 -m example.bam
20 (63.03Mbp)
>  88.57% │ ▆▃  ▄▂▄     ▄▃▁▄▆     ▃█ ▃▁  ▆▅▂▃     ▄ ▃▃  ▄▃  ▄  ▁  ▂▄  █▂│ Number of reads: 7577
>  78.73% │▅███▁███▅▆ ▇▂█████▂▂  ▃██▆██ █████▇▅▆▅▅████▆███▄▄█▂▃█  ██ ▂██│     (887 filtered)
>  68.89% │██████████ █████████▄▇███████████████████████████████▅▄██▄███│ Covered bases:   87.5Kbp
>  59.05% │█████████████████████████████████████████████████████████████│ Percent covered: 87.54%
>  49.21% │█████████████████████████████████████████████████████████████│ Mean coverage:   5.33x
>  39.37% │█████████████████████████████████████████████████████████████│ Mean baseQ:      33.7
>  29.52% │█████████████████████████████████████████████████████████████│ Mean mapQ:       54.3
>  19.68% │█████████████████████████████████████████████████████████████│
>   9.84% │█████████████████████████████████████████████████████████████│ Histo bin width: 1.6Kbp
>   0.00% │█████████████████████████████████████████████████████████████│ Histo max bin:   98.414%
     100.0K    116.4K    132.8K    149.2K    165.6K    181.9K     200.0K

    어떤 옵션이 있는지 살펴볼까요?
    -r [chromosome:start-end]: -r 옵션으로 영역을 지정해줄 수 있고, 영역을 지정해주지 않으면 모든 chromosome 에 대한 결과를 각각 만들어서 제공해줍니다.
    -m : table 형태의 결과 대신 histogram 을 보여줍니다.
    -w [bin 개수]: historgram 을 몇 개의 bin 으로 나눠서 보여줄지 결정합니다.
    -o [output 파일]: 결과를 output 파일에 저장합니다.
    -b [bam 목록 파일]: 여러 bam 파일을 사용할 때, bam 파일들의 이름이 적힌 목록 파일을 사용할 수 있습니다.

  3. depth
    지정해준 영역 혹은 position 에 대해 read depth 를 계산해줍니다.
    영역은 bed 파일로 제시해주면 되는데요, 만약 제공하지 않으면 bam 파일이 걸친 모든 영역의 position 마다 결과를 제공합니다.
    그럼 매우 큰 결과가 나오게 되니 주의해주세요!
    다음 명령어를 사용합니다: 
    $ samtools depth -b example.bed example.bam -o example_depth.txt

    옵션들을 한번 살펴보시죠.
    -a : 모든 position 에 대해 결과를 제공합니다. 이는 depth 가 0 인 position 도 포함합니다.
    -aa : 사용하지 않은 reference sequence에 있는 영역까지, 정말 전체 영역에 대한 결과를 모두 제공합니다.
    -b [bed 파일]: 결과를 만들 영역을 bed 파일로 제시해줍니다.
    -f [bam 목록파일]: 여러 bam 파일을 사용할 때, 파일 이름을 목록으로 만들어 제시합니다.
    -o [output 파일]: 결과를 output 파일에 저장합니다.
    -g [flags]: 특정 flag 를 가진 read 는 제외합니다.

  4. flagstat
    bam 파일의 기본적인 상태를 보여줍니다.
    주로 몇 개의 read 중 몇 개가 mapping 되었는지 등의 결과를 확인하고자 할 때 쓰입니다.
    flagstat 에 대한 더 자세한 설명은 이 포스트 를 참고해보세요!
    아래 결과를 보시면 더욱 이해가 쉬우실 겁니다. 
    $ samtools flagstat example.bam
4579959 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
371078 + 0 duplicates
4532991 + 0 mapped (98.97% : N/A)
4569463 + 0 paired in sequencing
2284756 + 0 read1
2284707 + 0 read2
4413160 + 0 properly paired (96.58% : N/A)
4475527 + 0 with itself and mate mapped
46968 + 0 singletons (1.03% : N/A)
21515 + 0 with mate mapped to a different chr
15292 + 0 with mate mapped to a different chr (mapQ>=5)

    결과를 살펴보면, 약 450만 개의 read 가 bam 파일에서 발견되었군요.
    그 중 371078 개는 duplicate 으로 밝혀졌구요, 대부분의 reads (98.97%) 가 reference genome 에 mapping 된 것을 알 수 있습니다.

  5. idxstats
    bam 파일의 통계 결과를 index 파일에서 찾아 출력해줍니다.
    즉, 이 명령어를 쓰기 위해서는 samtools index 가 선행되어야 합니다.
    결과는 다음 네 개의 컬럼으로 나옵니다.
    1) reference sequence name, 2)sequence length, 3) # mapped read-segments, 4)# unmapped read-segments
    다음 명령어로 실행합니다. 
    $ samtools idxstats example.bam
1	249250621	0	0
2	243199373	0	0
...
19	59128983	0	0
20	63025520	4532991	46968
21	48129895	0	0
22	51304566	0	0

    제 예시 파일에는 20번 염색체에 대한 정보밖에 없어서 결과가 하나밖에 안나오는군요. ^^

  6. phase
    이 기능에 대해서 알려면, 먼저 phasing 을 이해하고 있어야 합니다. Wikipedia 링크
    간단히 말하면, 우리가 얻은 genotype 이 두 상동염색체 중 어느쪽에서 왔는지 알아내는 과정을 phasing (혹은 haplotype estimation) 이라고 합니다.
    일반적인 시퀀싱 정보로는, 각 reads 들이 두 상동염색체 중 어느쪽에서 온건지 알 수 없습니다.
    그런데 특정 시퀀싱 기술을 통해서, 각 염색체에서 온 reads 들을 서로 다른 barcode 를 달아주면 reads 들을 두 그룹으로 나눌 수 있습니다.
    이 정보를 haplotype이라고 하는데요, 이런 실험기법을 통해 reads 들의 haplotype 정보를 알 수 있게 되는거죠. (이 과정을 phasing 이라고 합니다.)
    samtools phase 명령어는, 만약 우리가 가진 bam 파일이 phasing 되어 있는 파일이라면, bam 파일에서 변이를 찾아주고 이 중 heterozygous SNPs 에 대해서 haplotype 을 표시해주는 기능입니다.
    다음 명령어로 사용합니다: 
    $ samtools phase example.bam > example.phase

    만약 우리가 가지고 있는 bam 파일의 2,3,8번째에 heterozygous SNPs 들이 존재하고, 한쪽 allele (0번) 은 각각 A,C,T 이며 다른쪽 allele (1번) 은 각각 T,A,C 라면 결과물은 다음과 같이 나옵니다.

    CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
CC
CC
PS      frame-1   2    9
M1      frame-1  60    2      A       T       1       3       0       23      0
M1      frame-1  60    3      C       A       2       35      1       24      0
M1      frame-1  60    8      T       C       3       46      0       30      0
  7. stats
    stats 명령어는 bam 파일의 통계치를 텍스트 파일로 제공해줍니다.
    위의 flagstat 가 alignment 관련 정보를 제공했다면, stats 명령어는 checksum, read quality, ATGC contents 등의 정보를 보여줍니다.
    다음 명령어로 사용합니다: 
    $ samtools stats example.bam > example.stats

    output 은 텍스트 파일로 나오므로, 이를 “>” 로 저장해주는 것이 좋습니다.
    결과로 나오는 항목은 매우 많습니다. 링크를 참조해보세요.
    몇가지를 확인해보면 다음과 같습니다:
    CHK: 파일 checksum
    SN: Summary numbers
    FFQ: First fragment qualities
    LFQ: Last fragment qualities

Viewing: bam 파일을 실제로 읽어보고 싶을 때

  1. flags
    flags 는 bam 파일의 FLAG 정보를 설명해줍니다.
    bam 파일에 존재하는 두번째 column 은 flag 값인데요, 숫자 하나로 되어있지만 사실은 2진법 정보의 합으로 되어있는 것입니다. (링크 참조)
    예를 들어, 어떤 read 가 pair 되어 있고, read1 이 unmapped 상태라는 두 가지 상태에 있다고 가정합시다.
    그러면 해당 read 는 1 + 4 = 5 라는 값을 가지게 됩니다.
    이런 FLAG 에 어떤 이진법 조건이 있는지 samtools flags 를 통해 알 수 있습니다.
    다음 명령어로 사용합니다: 
    $ samtools flags
Flags:
0x1	PAIRED        .. paired-end (or multiple-segment) sequencing technology
0x2	PROPER_PAIR   .. each segment properly aligned according to the aligner
0x4	UNMAP         .. segment unmapped
0x8	MUNMAP        .. next segment in the template unmapped
0x10	REVERSE       .. SEQ is reverse complemented
0x20	MREVERSE      .. SEQ of the next segment in the template is reversed
0x40	READ1         .. the first segment in the template
0x80	READ2         .. the last segment in the template
0x100	SECONDARY     .. secondary alignment
0x200	QCFAIL        .. not passing quality controls
0x400	DUP           .. PCR or optical duplicate
0x800	SUPPLEMENTARY .. supplementary alignment
  2. tview
    tview 명령어는 interactive mode 로 특정 genomic region 에서 reads 들의 상태를 볼 수 있게 해줍니다.
    특정 위치를 지정하고 upstream, downstream 으로 이동하면서 reads 들이 어느 위치에 mapping 되어 있는지, 변이가 있는지, read quality는 어떤지 등을 직접 확인할 수 있죠.
    일반적으로는 genome 이 너무 크기도 하고, bam 파일에 있는 read 수도 매우 많기 때문에 잘 사용하지는 않지만, 가끔 variant call 이 제대로 되었는지 확인하는 용도로 사용할 수 있습니다.
    tip. 이 명령어를 사용할 때는 reference genome 을 같이 제공해주는 것이 좋습니다.
    다음 명령어로 사용합니다: 
    $ samtools tview example.bam reference.fa

    명령어를 실헹하고 나면, 터미널 창에 아래와 같은 검은 화면이 뜰 겁니다.
    여기서 slash ‘/’를 입력하고 20:341931 등의 위치를 입력하면 원하는 위치의 view를 볼 수 있습니다.
    키보드 방향키로 움직일 수 있고, 해당 위치에 mapping 된 reads 정보 및 variant 정보를 확인할 수 있죠.
    연속된 글자들이 하나의 read를 나타냅니다.
    다음과 같은 화면이 나타나게 됩니다:

     341931    341941    341951    341961    341971    341981
ccaggctctgtactcagtccatcacatacattaggtcttgcagtcctcataccaccataaggt
....................R..........................................
, ..................G.................
.......................... ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,g,,,,,,
,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,,,
,                   g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
  3. view
    view 명령어는 두 가지 용도로 사용됩니다.
    첫 번째로는 bam 파일의 데이터를 특정 조건(위치, quality 등등)에 맞는 것만 뽑아내는 용도입니다.
    두 번째는 파일 형식을 BAM-SAM-CRAM 으로 바꿔주는 용도입니다.
    BAM -> SAM 파일도 가능하고, SAM -> BAM 으로도 가능합니다.
    기본적인 명령어는 다음과 같습니다: 
    # 1번 용도
$ samtools view example.bam 2:10000-50000 -o example.sub.bam
# 2번 용도 1 - bam 파일을 sam 파일로 바꿀 때
$ samtools view example.bam -o output.sam
# 2번 용도 2 - sam 파일을 bam 파일로 바꿀 때
$ samtools view -b example.sam -o output.bam

여기까지 samtools 설치부터 indexing, editing 그리고 viewing 까지 살펴보았습니다.
명령어가 무척 많아서 한 글에서 다 다루는게 매우 힘들군요.
부족한 부분이 많을텐데, 앞으로 조금씩 보완해나가도록 하겠습니다.
samtools view, flagstat 등은 여러 방면에서 유용하게 쓰이는 명령어이므로 알아두시면 편리하실 거에요.
그럼 다음 시간에 만나요!