잠깐, BBMap인가 혹은 BBTools인가? 나도 헷갈린다.
추출해야 할 read가 1로부터 시작하는 번호의 목록으로 주어진다면 다른 방법을 써야 한다. PhyloSift에서는 시퀀싱 리드로부터 marker gene profile에 붙는 일부분만을 떼어내서 quality가 제거된 fasta 파일만을 제공한다. 나는 원본 파일로부터 fastq read의 부분집합을 떼어내고 싶었다. PhyloSift 결과물에서는 원래의 read ID가 사라지고 일련번호가 붙는다.
스크립트를 새로 짜려고 생각했다가 '바퀴를 다시 발명하지 말라'는 금언을 떠올렸다. 누군가 멋진 방법을 이미 구현해 놓았을 것이 뻔하기 때문이다. 약간의 검색 끝에 또다시 Biostars에서 Pierre Lindenbaum의 아주 짧은 스크립트를 발견하였다. 리눅스의 기본 유틸리티를 기가 막히게 조합한 것이다. 약 세 달밖에 되지 않은 신선한(?) 질문과 대답이었다.
Question: Subset Fastq file from list of read numbers
스크립트는 매우 간단하다. 사실상 one-liner이다.
$ join -t $'\t' -1 1 -2 1 \ <(gunzip -c input.fq.gz | paste - - - - | awk '{printf("%d\t%s\n",NR,$0);}' | sort -T . -t $'\t' -k1,1 ) \ <(sort -T . numbers.txt) | \ sort -t $'\t' -k1,1n |\ cut -f 2- |\ tr "\t" "\n" > subset.fastq
이것을 그대로 복사한 다음 input.fq와 numbers.txt 및 subset.fastq를 사용자가 원하는대로 고쳐서 명령행에 붙여넣고 엔터를 치면 그대로 실행이 된다. 단, 스크립트 파일에 내용을 복사해 넣은 뒤 실행하면 안된다. 아마도 process substitution의 복잡한 사정 때문에 그러한 것으로 보인다. join의 두 인수에 해당하는 process substitution을 임시 파일로 대체하도록 스크립트를 만들어서 파일로 저장한 다음, 이를 실행하면 잘 돌아감을 확인하였고 실무에도 아주 유용하게 사용하였다. Thanks, Pierre!
paste와 awk의 신비한 결합에 대해서는 아무리 공부해도 부족하지가 않다. join 명령어도 마찬가지이다. 샘플 파일을 만든 뒤 파이프 기호('|')를 제거한 다음 어떤 출력물이 나오는지를 살펴본다면, 이 명령행이 어떻게 돌아가는지를 알게 될 것이다. Fastq 파일을 분해하여 한 read에 해당하는 네 줄을 한 줄로 표현하고(탭으로 구분), 여기에 다시 숫자를 붙이고, 숫자를 기준으로 정렬하고, 그리고 번호를 참조하여 서로 매칭이 되는 것을 join 명령어로 뽑아내고, 다시 new line 문자를 넣어서 fastq 파일 형태로 환원하고... 정말로 배울 것이 그득한 명령행이다. 오늘도 이렇게 Pierre Lindenbaum에게 한 수를 배웠다. 어쩌면 그는 내 블로그에서 가장 많이 언급하는 인물이 될지도 모르겠다.
awk를 이용한 join, 그리고 join 명령어(join lines of two files on a common field) 자체의 활용은 정말 요긴하게 쓸 수 있는 것이라서 사용법을 잘 익혀두면 파이썬이나 펄로 직접 코딩을 하는 수고를 덜 수 있다.
특별히 이 활용례에서 눈여겨보아야 할 것은 탭문자를 명령어의 인수로 제대로 제시하기 위해 $'\t'라고 표시하는 트릭이다. 이에 대해서는 예전에 블로그에 기록한 일이 있다(Bash에서 탭을 표현하기 어려운 경우의 해결 방법).
댓글 없음:
댓글 쓰기