khmer 유틸리티와 파이썬

요즘은 deep metagenome sequencing data로부터 미생물 유전체를 재구성하는 방법에 대해 공부를 하는 중이다. 내년에는 MinION 플랫폼도 경험을 좀 해봐야 하는데... 하루가 멀다하고 새로운 기술이 나오는 것이 그렇게 반갑지는 않다. 더 많은 프로젝트를 할 수 있다는 것은 분명한 사실이지만 그것이 항상 가치있는 일인지는 여전히 고민스럽기 때문이다.

khmer란 무엇인가? 이는 k-mer에 기반을 둔 서열 데이터의 분석 및 전환(transformation) 도구이다. 주로 metagenome이나 transcriptome의 de novo assembly를 위한 각종 전처리를 해 주는 파이썬 라이브러리와 스크립트의 모임이라고 보면 된다. 과연 어떤 기능들이 있는지를 홈페이지를 통해 알아보았다.

• normalizing read coverage ("digital normalization")
• dividing reads into disjoint sets that do not connect ("partitioning")
• eliminating reads that will not be used by a de Bruijn graph assembler;
• removing reads with low- or high-abundance k-mers;
• trimming reads of certain kinds of sequencing errors;
• counting k-mers and estimating data set coverage based on k-mer counts;
• running Velvet and calculating assembly statistics;
• optimizing assemblies on various parameters;
• converting FASTQ to FASTA;

나는 파이썬에 대해 잘 모르기 때문에 가끔 파이썬 스크립트를 실행하려면 당혹스러울 때가 많다. 파이썬의 버전, pip, setuptools, virtualenv 등에 관한 체계적인 지식이 없어서 늘 애를 먹는다. khmer의 설치 역시 그러하였다.  파이썬이 특히 혼동스런 것은 서로 다른 버전을 동시에 운용할 필요가 있기 때문이다. 내가 사용하는 리눅스 머신은 CentOS 6.7이고, 기본적으로 설치된 것은 파이썬 2.6.6이다. 하지만 khmer는 파이썬 2.7을 요구한다. 그러면 어떻게 해야 할까? 구글링을 통해서 유용한 자료를 하나 찾았다.

Python 2.6.6이 설치된 CentOS 6.2에서 Python 2.7.3을 설치하여 활용하기

distribute란 파이썬 패키지의 빌드, 설치, 업그레이드 및 삭제를 용이하게 하는 패키지이다(Setputools에게 자리를 내어주었다는데?). virtualenv란 파이썬 개발 프로젝트 단위로 서로 다른 버전의 파이썬을 유지하게 하는 도구인 것으로 생각된다. 

간혹 gcc의 버전이 너무 낮아서 문제가 되는 경우도 있다. 현재의 CentOS에 기본적으로 딸려오는 gcc는 4.4이다. 이를 4.7로 업그레이드하려면? 리눅스 배포판에 맞추어진 패키지의 버전을 임의로 올리는 것은 그렇게 간단하지는 않다. 다음 사이트를 보면 이러한 gcc 업그레이드 방법이 상세하게 나와있다. 별도의 리포지토리를 등록하여 devtoolset-1.1을 설치한 뒤 scl command를 실행하는 것이다.

How to upgrade gcc on CentOS

khmer 설치의 실제

khmer의 git clone을 만들어서 설치하려면 오히려 더 혼동스럽다. 관리자 권한으로서 devtoolset-1.1이 설치된 상황이라면 다음과 같이 하라. /usr/local/apps 아래에 적당한 디렉토리를 만들어서 virtualenv 환경을 만든 다음에 pip2로 khmer를 설치하는 것이 좋다. 관리자 권한이 없다면 OS X에서 설치하는 방법을 따르면 된다.

# cd /usr/local/apps
# mkdir khmer
# cd khmer
# python2.7 -m virtualenv khmerEnv
New python executable in khmerEnv/bin/python2.7
Also creating executable in khmerEnv/bin/python
Installing setuptools, pip, wheel...done.
# scl enable devtoolset-1.1 bash
# gcc --version
gcc (GCC) 4.7.2 20121015 (Red Hat 4.7.2-5)
Copyright (C) 2012 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

# source khmerEnv/bin/activate
(khmerEnv)# pip2 install khmer
Collecting khmer
Collecting bz2file (from khmer)
Collecting screed>=0.9 (from khmer)
Installing collected packages: bz2file, screed, khmer
Successfully installed bz2file-0.98 khmer-2.0 screed-0.9
(khmerEnv)#

[Next-Generation Sequencing] Mapping과 관련한 수치 구하기

SAM/BAM 파일에서 mapping과 관련한 수치를 뽑는 방법을 one-liner 철학에 입각하여 정리해 보도록 하겠다. 물론 다른 프로그램을 사용하면 좀 더 고차원적인 수치를 뽑을 수도 있을 것이다. 이 글을 쓰면서 추구하는 것은 되도록 적은 수의 기본 프로그램과 UNIX/LINUX에 설치된 유틸리티를 최대한 활용하여 (옵션이 길어짐은 피할 수 없음) 원하는 결과를 얻는 것이다. 이 문서를 작성하면서 인터넷 상의 여러 정보를 참조하였지만 특히 GitHub의 SAM and BAM filtering oneliner에서 많은 도움을 얻었다.

오늘의 작업을 위해서 쓰일 가장 핵심적인 도구는 samtools임은 부정할 수 없다. samtools는 1.x 버전으로 업그레이드되면서 공식 웹사이트도 변경되었음을 먼저 알아두도록 하자.

• SAMtools 0.1.19까지 (링크)
• SAMtools 1.x부터 (링크 http://htslib.org/)

1. 가장 일반적인 mapping statistics 출력하기(samtools flagstat)

$ samtools flagstat BL21.bam
869658 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
842935 + 0 mapped (96.93%:-nan%)
869658 + 0 paired in sequencing
434829 + 0 read1
434829 + 0 read2
791630 + 0 properly paired (91.03%:-nan%)
831272 + 0 with itself and mate mapped
11663 + 0 singletons (1.34%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Mapping statistics 수치를 뽑는 가장 간단한 방법은 samtools flagstat를 이용하는 것이다. 위에 든 사례(E. coli BL21)는 duplicate 하나 없이 너무나 깔끔하게 mapping이 된 결과물이라서 별로 공부할만한 것이 없다. IS가 잔뜩 박혀있는 Shigella boydii의 genome 서열에 대한 일루미나 매핑 결과물을 가지고 다시 분석을 해 보자. Bowtie2를 사용하되 시간을 줄이기 위하여 100만개의 read(50만 read pair)로 입력 데이터를 제한하였다.

$ bowtie2 -x SB_ref -1 SBpart_1.fastq -2 SBpart_2.fastq -S SB.sam
500000 reads; of these:
  500000 (100.00%) were paired; of these:
    55609 (11.12%) aligned concordantly 0 times
    402919 (80.58%) aligned concordantly exactly 1 time
    41472 (8.29%) aligned concordantly >1 times
    ----
    55609 pairs aligned concordantly 0 times; of these:
      11135 (20.02%) aligned discordantly 1 time
    ----
    44474 pairs aligned 0 times concordantly or discordantly; of these:
      88948 mates make up the pairs; of these:
        72973 (82.04%) aligned 0 times
        7334 (8.25%) aligned exactly 1 time
        8641 (9.71%) aligned >1 times
92.70% overall alignment rate

bowtie2.2.4 매핑 과정에서 약간의 수치가 출력된다. 이제 결과물을 BAM 파일로 전환하여 samtools로 다시 수치를 뽑아보자.

$ samtools view -b -S -o SB.bam SB.sam
[samopen] SAM header is present: 1 sequences.
[hyjeong@tube 01_mapping_test]$ /usr/local/apps/a5_miseq_linux_20141120/bin/samtools flagstat SB.bam
1000000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
927027 + 0 mapped (92.70%:-nan%)
1000000 + 0 paired in sequencing
500000 + 0 read1
500000 + 0 read2
888782 + 0 properly paired (88.88%:-nan%)
919416 + 0 with itself and mate mapped
7611 + 0 singletons (0.76%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

 어라? duplicate가 하나도 나타나지 않는다. 같은 데이터를 가지고 bwa 0.6.2-r126을 돌려 보았지만 mapped read의 수치가 아주 조금 줄어들었고 이번 포스팅에서 설명을 하려 했었던 duplicates 수치가 역시 하나도 나오지 않았다. "duplicates" 수치가 의미하는 것은 read가 반복하여 매핑된 횟수들의 합을 의미한다.

2. Unmapped reads

이제부터는 alignment 파일(SAM or BAM) 내에 있는 FLAG를 참조하여 적절히 필터를 하는 기법이 필요하다. samtools view 명령으로 특정 조건에 맞는 read alignment만을 출력하거나(-f INT) 혹은 그 조건에 맞는 것을 걸러서 제거하여 출력하는 것(-F INT)이 가능하다. 예를 들어 0x04(0x가 앞에 붙은 것은 16진수를 의미한다. 따라서 0x04는 십진수의 4와 같다) FLAG은 "segment unmapped"를 의미한다. 다음의 두 명령은 전혀 다른 수치를 출력한다. 

samtools view -f 0x04 -c align.bam # unmapped read의 수
samtools view -F 0x04 -c align.bam # mapped read의 수

좀 더 정확히 말하자면 mapping이 일어난 경우 mapped read의 수라기보다는  mapped position이라고 하는 것이 나을 것이다. 하나의 read가 reference genome 내의 여러 위치에 정렬하는 경우 SAM/BAM 파일 내에 각각 별개의 레코드로 존재하게 되기에 그렇다. -c는 matching record의 카운트만을 출력하는 것이다. 따라서 -b -o out.bam이라고 하면 read/alignment를 BAM 파일로 출력함을 의미한다. 또한 -f 0x04는 -f4 -f0x04 등의 형태로 사용함을 허락한다.

3. Uniquely mapped read의 수 세기(이 항목은 잘못되었음 2019-08-06)

samtools view -F 0x40 aligned_reads.bam | cut -f1 | sort | uniq | wc -l
samtools view -f 0x40 -F 0x4 aligned_reads.bam | cut -f1 | sort | uniq | wc -l #left mate
samtools view -f 0x80 -F 0x4 aligned_reads.bam | cut -f1 | sort | uniq  | wc -l #right mate

4. Paired end read의 정렬 결과로부터 insert size 추정하기

이것은 one-liner로서는 뽑기가 쉽지 않다. 다행히 getinsertsize.py(GitHub link)라는 훌륭한 공개 파이썬 스크립트가 있다. 사용법은 스크립트 저자의 블로그(2012년 4월 6일 포스팅)에 잘 설명되어 있다. SAM/BAM 파일을 표준 입력으로 받을 수 있어서 매우 깔끔한 처리가 가능하다. getinsertsize.py SAM_file 또는  samtools view BAM_file | getinsertsize.py - 어느 것이든 실행 가능하다.

$ samtools view aln-pe.bam | getinsertsize.py -
1M...  # read의 전체 수를 의미하는 것으로 생각된다
Read length: mean 101.0, STD=0.0
Read span: mean 394.062581504, STD=31.8156907927

5. Coverage 관련 수치 구하기

Reference genome 전체에 대한 평균 coverage는 얼마일까? zero coverage region의 위치를 추출할 수는 없을까? 솔직히 말해서 one-liner로 이런 수치를 계산한다는 것은 다소 무리가 아닐까 싶다. 구글에서 bam sam read depth coverage라고 검색을 해 보면 Biostars와 SEQanswers에 이에 대한 몇 가지 질문과 답이 있으니 참고하기 바란다.

[Next-Generation Sequencing] Variant calling 작업의 오류 해결과 Bowtie 2

다음 달에 있을 미생물 유전체 정보 분석 교육에 사용할 실습 자료를 만들고 있다. 워낙 많은 양의 시퀀싱 결과물이 있어서 적당한 것을 골라서 정해진 시간 내에 실습이 무난히 이루어지도록 분량을 줄이고 프로그램 버전을 확인하는 작업을 한창 진행 중에 있다. 어쩌다보니 이번 가을에는 개인적으로 몇 사람의 연구원에 대한 교육을 시작하여 이제 마무리 단계에 있고, UST 강의까지 겹쳐서 그동안 축적한 정보를 정리하는데 아주 좋은 기회가 되고 있다.

다만 아쉬운 점이 있다면 나는 실무에서 상용 tool인 CLC Genomics Workbench를 써서 대부분의 분석 작업을 실시하지만, 교육에서는 그럴 수가 없다는 것이다. 따라서 사용이 다소 불편하더라도 공개된 소프트웨어를 이용하는 방식으로 교육을 하지 않을 수 없다. 상용 툴이 모든 것을 다 해주는 것은 아니기에 공개 툴을 조금씩 쓰기는 하지만 이를 남에게 전수하려면 더욱 많은 사전 경험이 필요하므로 실수가 없도록 철저하게 매뉴얼을 숙지하고 점검에 점검을 거듭하고 있다.

어제는 미생물 NGS 데이터를 이용한 variant calling 작업을 하다가 오류를 만났다.

$ samtools mpileup -uf ref.fa BL21.sorted.bam | bcftools view -bvcg - > var.raw.bcf
[mpileup] 1 samples in 1 input files
Set max per-file depth to 8000
[bcf_sync] incorrect number of fields (0 != 5) at 0:0
[afs] 0:0.000

빨강색으로 표시한 것이 오류 메시지이다. 구글을 열심히 뒤져 보았으나 이것과 딱 맞는 사례는 없었다. 도무지 이유를 알 수가 없었다. 필드의 수가 5개가 아니고 6개라는 오류 메시지를 보고한 경우는 매우 많았으나 나처럼 아예 아무것도 나오지 않은 사례는 보이질 않았다. 문제 해결을 위해 이것저것 한참을 알아보다가 내 시스템에 몇 가지 프로그램이 버전이 맞지 않은 상태로 뒤섞여 존재함을 알 수 있었다.

문제는 바로 de novo assembly pipeline인 A5-miseq 때문이었다. 내가 사용하는 A5-miseq에는 0.1.99-4428cd 버전의 samtools, bcftools 및 vcfutils.pl이 포함되어 있다. 그런데 별도로 설치한 samtools 1.2와 뒤섞여 실행이 되면서 오류가 난 것이다. samtools 0.1.18에는 bcftools/vcftools.pl이 같이 포함되어 있었지만 1.* 대의 버전으로 올라오면서 각자 따로 배포되고 있다. 이에 대해서는 최신 배포판 사이트를 참조하면 된다. Scaffolding tool인 SSPACE에도 bwa와 bowtie가 같이 들어있다. 이를 설치할 때에는 시스템에 먼저 깔려있던 것들과 서로 혼동되지 않도록 주의가 필요하다. 아무 생각 없이 bwa를 실행하였을 때 도대체 어느 경로에 설치된 바이너리가 돌아가고 있는지는 알고 있어야 한다는 뜻이다.

그렇다면 /usr/local/bin에 위치해서 나를 혼동스럽게 했던 samtools 1.2를 지우고 A5-miseq에 포함된 samtools와 bcftools만을 이용하여 variant calling 작업을 해 보자. 이제는 아무런 오류 없이 성공적으로 수백개의 variant가 출력되었다. 이런 성가신(?) 작업이 CLC Genomics Workbench에서는 일관적인 GUI 환경 안에서 편리하게 수행되지만 편리한 만큼 적지 않은 돈이 드는 선택이라 강요할 수는 없는 노릇이다. Command line interface 위주의 공개형 도구는 비록 사용은 불편할지라도 분석 과정의 기본을 익히기에는 좋으니 이것과 담을 쌓고 살 수는 없는 노릇이다. 

Variant Calling 과정의 상세

글을 작성한 김에 변이 추출(variant calling) 과정에 대하여 좀 더 자세히 기술해 보고자 한다. 변이 추출은 유전체 resequencing의 가장 기본적인 목적이다. copy number variation이나 structural variation처럼 까다로운 것은 다음으로 미루고, SNV의 추출에 대한 것만을 정리하겠다.

매핑 도구는 bowtie2를 쓰기로 한다. 따라서 bowtie2에 포함된 공식 매뉴얼을 참조하는 것이 바람직할 것이다.

• Using SAMtools/BCFTools downstream
• Fast gapped-read alignment with Bowtie 2 (Nature Methods 2012년)

bowtie2로 만들어진 SAM 파일을 정렬 및 인덱싱을 하여 samtools pileup 명령으로 처리하는 것이 기본이다. ref.fa 파일은 미리 samtools faidx 명령으로 인덱스를 생성해 두어야 한다.

samtools mpileup -u -f ref.fa sample.sorted.bam > sample.bcf
bcftools view -v -c -g sample.bcf > sample.raw.vcf 

파이프를 사용하면 한 줄의 명령어로 줄여서 실행할 수도 있다.

samtools mpileup -uf ref.fa sample.sorted.ban | bcftools view -bvcg - > sample.raw.bcf

이쯤에서 SAMtools의 공식 문서를 참조해 보자.

Calling SNPs/INDELs with SAMtools/BCFtools

SAM format specification

위 단계에서 추출한 변이를 100% 전부 믿고 받아들일 수는 없는 노릇이다. 이제 vcfutils.pl로 적절히 필터링하여 정확한 결과만을 남기도록 한다. read depth가 100 이상인 곳에 대해서만 필터링하려면 다음과 같이 한다.

bcftools view sample.raw.bcf | vcfutils.pl varFilter -D100 > sample.flt.vcf

Short read aligner의 최강자는 누구인가?

여기에 대해서는 꽤 많은 논문이나 포스팅이 있으니 몇 개를 간추려서 소개하는 것으로 갈음하도록 하자.

• (쉬운 글) BWA, SOAP, MAQ, BLAT...OMG Alignment Methods! (2014년)
• Benchmarking short sequence mapping tools (BMC Bioinformatics 2013년 7월 7일)
• DNA aligner accuracy: BWA, Bowtie, Soap and SubRead tested with simulated reads (2014년 11월 4일)