FASTQ 파일에서 목록에 없는 read를 추출하기

수십만 read 이상이 수록된 Fastq 파일에서 특정 read만을 추출하고 싶을 때가 있다. 추출해야 할 read의 ID가 별도의 텍스트 파일(read.list)에 들어있다면 seqtk을 사용하여 다음과 같이 실행하면 된다.

seqtk subseq in.fastq read.list > out.fastq

이번에는 반대의 경우를 생각해 보자. read.list에 없는 read만을 in.fastq 파일에서 추출하려면? read.list를 읽어서 read ID를 key로 삼는 해쉬를 만든 뒤, in.fastq를 읽어들여서 존재하지 않는 read ID를 만나면 이것을 4줄 단위로 출력하는 펄 스크립트를 짜야 하나? 일단은 검색에 의존해 보기로 했다. 역시 누군가는 앞서서 나와 비슷한 고민을 했고, 성실한 네티즌들은 이에 대한 해결책을 이미 제시하였다.

[Biostars] How to remove a list of reads from fastq file?

질문자는 매우 기본에 충실한 방법으로 이 문제를 해결하고 있었다. 우선 fastq 파일로부터 read ID를 추출하여 목록 파일을 만든다.

zcat in.fastq.gz | awk 'NR%4==1' | sed 's/@//' > in.fastq.readsID

그러고 나서 in.fastq.readsID에 해당하지 않는 read의 ID를 추출한다. grep에서 -v 옵션을 주면 match가 되지 않는 라인을 출력한다(all lines but those matching are printed). 나는 부끄럽게도 이 기능을 아직까지 전혀 모르고 있었다! 미련하게도 이런 종류의 일에서는 늘 Perl의 해쉬를 사용하려고 했었으니...

grep -f reads-remove-list in.fastq.readsID -v > remaining.list

마지막으로 seqtk subseq를 사용한다.

seqtk subseq in.fastq.gz remaining.list | gzip - > remaining.fastq.gz

손이 많이 가는 방식이지만 기본에 매우 충실하다. awk를 사용하여 4줄마다 한번씩 줄을 출력하고, sed를 써서 서열 ID 앞부분의 '@'를 제거하고, 마지막으로 grep으로 match되지 않은 라인을 출력하고... 인스턴트 커피 가루와 프리마, 설탕을 티스푼으로 정성스럽게(그러나 취향에 맞게) 떠서 커피를 타듯이 처리를 한 것이다. 질문자는 이렇게 하는 것이 너무 번거롭다고 생각하여 믹스커피 봉지를 뜯듯이 좀더 쉽게 파일 처리를 하는 방법을 물은 것이다.

제시된 해결방법 중 가장 마음에 드는 것은 BBMap의 filterbyname.sh 스크립트를 이용하는 것이다. BBMap은 미국 Joint Genome Institute에서 일하는 Brian Bushnell이 만든 종합선물세트 비슷한 것이다. 지금은 BBTools라는 확장된 패키지로 불리는 것 같다. filterbyname.sh가 좋은 점은 fastq file pair를 동시에 다룰 수 있다는 점이다. seqtk subseq를 사용하려면 infile_1.fastq와 infile_2.fastq에 대해서 각각 작업을 해야 된다. 맨 뒤의 include=f가 바로 목록에 있는 것을 제외하고 출력하게 만드는 스위치인 셈이다.

filterbyname.sh in=infile_1.fq in2=infile_2.fq out=outfile_1.fq out2=outfile_2.fq names=readID.list include=f

단, 오래전에 시퀀싱한 결과물이나 ART로 생성한 simulated read pair를 다룰 때에는 조심해야 한다. read ID의 끝부분에 /1 또는 /2의 tag가 붙어있기 때문이다. 요즘의 일루미나 시퀀싱 결과물에서는 read ID와 공백을 두고 direction tag이 존재하므로 read ID 목록 파일을 이용하여 원본 fastq 파일을 다루기에 좋다.

(old style) @HWI-EAS390_0001:7:120:12431:20961#0/2
(ART style) @chromosome_ST307PT01-2132000/1
(new style) @M00220:56:000000000-AML1R:1:1101:12360:1835 1:N:0:5