2022년 11월 11일 금요일

Khmer recipe의 뒤늦은 발견

생명정보 분야에서 Khmer(Read the Docs)는 '크메르'(802년부터 1431년까지 캄보디아에 존재한 제국)가 아니고 k-mer 분포를 이용하여 NGS read를 분석하는 도구이다. 한동안 파이썬 2.7 대에 머물러 있어서 설치하기가 상당히 귀찮았었지만 지금은 bioconda 채널에도 등록되어 있고 심지어 우분투 패키지로도 제공된다.

Khmer와 경쟁하는 소프트웨어라면 KAT을 들 수 있다. 실행 속도도 빠르고, k-mer abundance spectrum과 GC content를 한꺼번에 보여주는 효과적인 plot을 제공하는 면에서는 우수하다. 그러나 k-mer abundance를 이용하여 read를 필터링하는 방법은 명확하지가 않다. 내 의견으로는 khmer가 이 기능에서 압도적으로 우수하다. Jellyfish는 단순한 k-mer counter이다.

최근 아주 까다로운 장내 미생물의 NGS data를 다루면서 khmer와 smalt(read mapper) 및 tablet(sequence alignment & assembly의 시각화 도구)의 기능을 철저히 복습하게 되었다. khmer의 부속 스크립트를 활용하여 read를 필터링하거나 k-mer abundance spectrum을 플로팅하는 자작 스크립트를 약 6년 전에 만들어서 이따금 사용해 왔다. 그런데 khmer 명령어의 상세한 사용법을 인터넷에서 검색해 보다가 아주 훌륭한 레시피가 있음을 뒤늦게 발견하게 되었다.

Welcome to the khmer-recipes site!

예를 들자면 'load-into-counting.py'로 countgraph를 만든 뒤 이로부터 'abundance-dist.py'로 histogram을 생성한 다음, 그림을 그리기 위해 gnuplot을 수작업으로 돌렸었다. 그런데 공식적으로 제공되는 스크립트 중에 'plot-abundance-dist.py'가 있었다. 사용자가 지정한 k-mer coverage의 구간 정보를 이용하여 read를 추출하는 'slice-reads-by-coverage.py'라는 똑똑한 스크립트가 있다는 것도 오늘 알았다. 이전에는 k-mer coverage threshold를 주고 이를 만족하는 read를 추출하는 'filter-abund.py'만 사용했었다.

KAT의 read filtering 방법은 조금 다르다. 필터를 적용할 k-mer hash를 먼저 확정해 놓은 다음, 이를 지정된 퍼센트 이상으로 포함하는지 여부를 가지고서 필터링을 한다. 그래서인지 칼로 무를 썰듯이 깨끗한 필터링이 잘 되지 않았다. 물론 내가 매뉴얼을 정확히 이해하지 못해서 깊숙한 곳에 숨어있는 고급 기능을 모르고 있을 수도 있다.

NGS 자료를 다루는 사람이라면 k-mer를 이용한 소프트웨어를 다루지 않을 수 없는 것이 현실이다. De novo assembly, sequence search, genome 특성의 파악 등 쓰이지 않는 곳이 없다. 그러나 이를 정확히 이해하고 다른 사람에게 설명할 수 있는가? 그것은 별개의 문제이다. 생화학을 모른다고 해서 음식을 씹어 먹고 소화할 자격이 없는 것은 아니지 않는가? 유전학을 몰라도 아이를 낳는데 아무런 문제가 없듯이 말이다. 가만있자... 이는 적당한 비유가 아니로구나. 소화가 잘 되지 않고 아이가 잘 생기지 않는 원인을 알고 해결책을 마련하는 데에는 전문 지식을 하는 것이 도움이 되기는 하겠구나.

2022년 11월 12일 업데이트

khmer recipes site에서 소개한 일부 스크립트(예: slice-reads-by-coverage.py, plot-abudance.py)는 우분투 deb 패키지에 포함된 스크립트(/usr/lib/khmer/bin/)에 더 이상 존재하지 않는다.

댓글 없음: