k-mer 분석 프로그램인 jellyfish를 미련하게 쓰고 있었다

khmer, jellyfish, KAT... NGS read로부터 k-mer 분석을 할 때 내가 즐겨 사용하는 프로그램들이다. 초창기에는 khmer를 즐겨 사용했었다. 이 프로그램은 특히 k-mer에 대해 어떤 기준을 적용하여 필터링을 할 수 있다. 그러나 모든 작업에서 이렇게 데이터를 잘라내는 '침습적' 조작을 필요로 하지는 않는다. 요즘은 khmer의 핵심 프로그램보다는 부속 파이썬 스크립트를 쓰는 일이 많다. 예를 들어 readstats.py나 interleave-reads.py 같은 것.

그다지 크지 않은 데이터셋에 대해서 간편하게 k-mer abundance 분석을 하려면 jellyfish가 더 편리하다. 최근 K-BDS에 등록할 '묵은지와 같은' 옛날 NGS 데이터(일루미나 HiSeq 2000)을 재정비하다가 jellyfish를 아주 무식하게 쓰고 있다는 것을 깨달았다. 사실 Getting started 문서만 제대로 읽었어도 이런 비효율적인 일은 하지 않았을 것이다. 이 비효율적인 실행 방식은 몇 년이나 작동하고 있었다. 이 모든 실수는 게놈 고물상 영업을 개시하면서 비로소 드러난 것이다.

k-mer spectrum을 시각화하려면 'jellyfish count'로 k-mer를 센 뒤, 'jellyfish histo'로 히스토그램을 만든 다음 gnuplot으로 적당히 그림을 그리면 된다. 그런데 미련하게도 첫 명령어에서 산출된 파일에 대하여 'jellyfish dump'로 모든 k-mer의 수를 수록한 텍스트 파일을 뽑은 뒤 다시 여기에서 awk/uniq/sort 조합으로 히스토그램을 만들고 있었으니... 물론 이것도 나의 순수한 창작은 아니고 인터넷 어디선가 검색을 통해서 알아낸 것이었다. 이 미련한 한 줄 스크립트는 아래와 같다. SAMPLE.kmer21.txt는 'jellyfish count'의 결과 파일이다.

awk ‘{print $2}’ SAMPLE.kmer21.txt | sort -n | uniq -c | awk ‘{print $2 "," $1}' > SAMPLE.jf.hist

텍스트 파일 덤프는 정말 쓸데없는 일이었다. 처음부터 히스토그램을 뽑아내면 되는 것이었다. Fwd 및 rev read를 하나로 합쳐서 interleaved fastq file(*.pe.fq)을 만들어 놓은 다음 k-mer abundance plot을 그리는 방법을 알아보자. jellyfish 명령어로 21-mer를 세어서 히스토그램을 만든 뒤 gnuplot에서 플로팅하는 전체 명령어를 다음과 같이 재정비해 보았다. 

for x in *pe.fq
do
  echo Processing $x...
  x=${x%%.pe.fq}
  jellyfish count -m 21 -s 100M -t 12 -C $x.pe.fq -o $x.counts.jf
  jellyfish histo -o $x.jf.hist $x.counts.jf
  echo Running gnuplot...
  echo set term png > $x.jf.gp
  echo set output \"$x.jf.png\" >> $x.jf.gp
  echo set logscale x >> $x.jf.gp
  echo set logscale y >> $x.jf.gp
  echo set grid >> $x.jf.gp
  echo set xlabel \"21-mer frequency\" >> $x.jf.gp
  echo set ylabel \"number of distinct 21-mer\" >> $x.jf.gp
  echo set key off >> $x.jf.gp
  echo set title \"$x\" >> $x.jf.gp
  echo plot \"$x.jf.hist\" using 1:2 with lines >> $x.jf.gp
  gnuplot $x.jf.gp
done

실행을 마치면 *.jf.png라는 파일이 생겨날 것이다. 히스토그램 파일의 필드 구분자가 만약 콤마인 경우에는 'echo set datafile separator \“,\” » $x.jf.gp'라는 줄이 'echo plot...' 전에 들어가야 한다. 

히스토그램 파일을 넣어주면 분석 그림을 그려주는 GenomeScope라는 웹서비스도 있었다. 2020년에 Nature Communitations에 발표된 논문 제목은 'GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes'. 교싲신저자인 M.C. Schatz의 이름이 아주 낯이 익다. 어디서 봤더라... 논문 목록을 훑어보다가 'Hawkeye and AMOS'라는 2011년도 논문(링크)을 찾아냈다. 내 블로그에도 몇 차례 언급되었던 소프트웨어이고 사용했던 경험도 있다.