readstats.py는 너무 느리다

Fastq 파일이 얼마나 많은 read와 base를 갖고 있는지 확인하기 위해 습관적으로 readstats.py를 쓰고 있었다. 문제는 너무 느리다는 것. 그런데 어느날 BBMap의 stash.sh를 read statistics 산출에 쓰면 어떨까하는 생각을 하게 되었다. 여러 파일을 다루어야 하니 statswrapper.sh를 쓰는 것이 옳다.

프로그램 실행 시간을 비교해 보았다. Fastq 파일의 QC(단지 statistics을 추출하는 것으로부터 quality trimming/adapter clipping 등의 적극적인 조작까지)를 위한 소프트웨어는 한두 가지가 아니겠지만, 그래도 내가 가끔 쓰는 FastQC까지를 포함해서 비교해 보기로 했다. 입력물로 사용할 fastq file 두 개는 총 811 Mbp 정도의 read를 갖는다. 화면 출력은 생략한다.

$ time readstats.py KCTC3164_1.fq.gz KCTC3164_2.fq.gz
real 1m41.304s
user 1m40.721s
sys 0m0.339s

$ time statswrapper.sh format=6 in=KCTC3164_1.fq.gz,KCTC3164_2.fq.gz
real 0m9.528s
user 0m9.888s
sys 0m0.329s

$ time /usr/local/apps/FastQC/fastqc KCTC3164_1.fq.gz KCTC3164_2.fq.gz
real 0m49.393s
user 0m57.077s
sys 0m2.492s

비교가 되지 않는 실행 속도이다. 당연히 statswrapper.sh를 쓰지 않을 이유가 없다. 그러나 출력 포맷에서 필요하지 않은 컬럼이 많으니 이를 적당히 정리하면 좋다. 여기에 awk를 이용하자. Contig가 아니라 scaffold에 대한 데이터를 출력하게 만들었다. 만약 contig에 대한 수치를 다루게 되면, N이 들어있는 read를 둘로 나눌지도 모른다는 생각이 들었기 때문이다. stats는 평균 길이에 대한 정보는 주지 않으니 awk 명령어 안에서 계산을 시키면 된다. statswrapper.sh를 실행할 때 나타나는 첫 줄(java -ea 어쩌고...)은 표준출력이 아니라 표준에러이므로 파일로 리다이렉션을 하면 파일에는 들어가지 않는다. 다음의 awk 사례에서는 file name(전체 경로를 제고), read 수, 총 bp 수, 평균 bp, 그리고 평균 GC 함량이다. format=4로 실행하면 scaffold에 대한 수치만 나오게 되고 '#'로 시작하는 헤더 라인도 나오지 않으니 아래의 사례(format=6)을 쓰는 것보다 awk 명령을 더 단순하게 작성할 수 있을 것이다.

$ statswrapper.sh format=6 in=KCTC3164_1.fq.gz,KCTC3164_2.fq.gz | awk -vOFS="\t" '!/^#/{sub("^.*/", "", $20); print $20, $1, $3, $3/$1, $18}' > readstats.txt
java -ea -Xmx200m -cp /usr/local/apps/bbmap/current/ jgi.AssemblyStatsWrapper format=6 in=KCTC3164_1.fq.gz,KCTC3164_2.fq.gz
[hyjeong@tube Fastq_46_final]$ cat readstats.txt 
KCTC3164_1.fq.gz 2687642 405833942 151 0.35005
KCTC3164_2.fq.gz 2687642 405833942 151 0.35068

타성에 젖어서 늘 하던 방식을 고집할 것이 아니라 때로는 살짝 일상에서 벗어나 보는 것도 필요하다.

내가 가고픈 그곳으로만 고집했지 - '변해가네' 가사에서