FASTX-Toolkit의 이상한 행동(fastq_to_fasta)

어제 남세균(Phylum: Cyanobacteria, NCBI Tax ID: 1117)의 MiSeq sequencing read를 다루다가 pair의 수가 맞지 않아서 무척 애를 먹었다. 그 원인은 FASTX-Toolkit에 포함된 fastq_to_fasta 명령어의 작동이 원활하지 않아서 그런 것으로 잠정적인 결론을 내렸었다. 오늘은 이 문제를 심층적으로 다루어 보려고 한다.

이제 고전이 된 FASTX-Toolkit은 2010년에 공개된 버전 0.0.13이 최신판일만큼 업데이트가 되지 않고 있다. 이는 다른 관점에서 보면 더 이상 손을 댈 필요가 없을 정도로 완벽하다는 뜻도 되겠다. 물론 마지막 버전 공개 당시에는 quality encoding이 phred64인 시퀀싱 결과물이 상대적으로 많았기에, 요즘의 표준인 phred33 데이터를 다루려면 -Q33이라는 옵션을 별도로 주어야 한다는 점이 좀 불편하다. FASTX-Toolkit을 능가할 수준의 유틸리티가 요즘 많이 있겠지만, 표준 입출력을 지원하여 파이프 구성이 편리하고, 심플하며 일관성있는 인터페이스를 갖고 있어서 아직도 쓸모가 많다.

어제 사용한 read의 길이 분포를 알아보자. MiSeq에서 유래한 것이라서 HiSeq read처럼 모든 서열의 길이가 일정하지 않다. 내가 원했던 작업은 280 bp 미만의 read는 제거하고, 그보다 긴 것은 앞부분 280 bp만을 취하여 paired read 형태로 정리하는 것이었다. 종종 고정된 길이의 read를 요구하는 후속 프로그램이 있었기 때문이다.

기준 길이 이상의 read만을 선별하고자 할 때 나는 SolexaQA++을 이용한다. CLC Genomics Workbench에서도 이 작업이 가능하지만(Microbial Genomics Module의 Fixed Leng Trimming) read를 불러들여서 작업을 하고 다시 파일로 내보내는 것이 귀찮아서 명령행 인터페이스를 이용하고자 한 것이다.

$ SolexaQA++ lengthsort -l 280 -c MA-KW_1.fastq MA-KW_2.fastq

리포트에 의하면  4244152 + 4244152 read가 만들어졌다고 한다. 실체 출력 파일을 확인해 보면 이 수치가 잘 맞는다. 그러면 다음과 같이 한 줄 명령어를 만들면 원하는 일이 다 이루어질 것으로 생각하였다. '_2' 파일에 대해서도 같은 방식으로 실행을 하였다.

$ fastx_trimmer -Q33 -l 280 –i MA-KW_1.fastq.paired | fastq_to_fasta -Q33 > temp_1.fastq

그런데 최종 파일이 수록한 read의 수가 달라졌다. 이 명령어에서는 필터가 전혀 작동하지 않으므로 최초에 공급한 파일의 read 수(4244152 + 4244152)를 그대로 유지해야 정상이다. 하지만 실제 결과 파일에서는 4244144 + 4244137가 되었다. file 1과 file 2에서 전부 read 수가 조금씩 줄어든 것이었다. 세부적으로 조사를 해 보니 fastq_to_fasta를 거치면서 read가 변했다. fastq_to_fasta 대신에 fastq2fasta.pl(by Brian J. Knaus)를 사용하면 read 수는 전혀 바뀌지 않았다.

좀 더 궁극적인 해결 방법은 없을까? 파일 조작의 마지막 단계에서 pair 여부를 다시 한번 점검하는 것이 필요할 것으로 생각된다. 이런 용도에 딱 맞는 도구가 무엇이 있을지 조사를 해 보았다.  DOE Joint Genome Institute의 BBTools에 포함된 여러 shell script 중에서 repair.sh를 쓰면 될 것 같다. BB는 이 프로그램을 만든 Brian Bushnell의 이니셜이다. 흥미를 가지고서 사용법을 조사해 보니... 쓸 수가 없다. 왜냐하면 repair.sh가 인식하는 입력 파일은 fasta가 아니라 fastq이기 때문이다. BBTools는 정말 쓸모가 많은 다양한 기능의 명령어를 포함하고 있으니 좀 더 시간을 두고 공부를 해 보자.

FASTX-Toolkit의 fastq_to_fasta가 가장 중요한 용의자이므로, 표준 입출력을 지원하는 다른 종류의 포맷 전환기를 사용해 보자. seqtk가 적당하겠다.

$ fastx_trimmer -Q33 -l 280 -i MA-KW_1.fastq.paired | seqtk seq -a > temp_1.fa
$ fastx_trimmer -Q33 -l 280 -i MA-KW_2.fastq.paired | seqtk seq -a > temp_2.fa

4244152 + 4244152 read라는 완벽한 결과가 나와주었다.