내가 즐겨사용하는 readstats.py는 FASTA/FASTQ 파일을 읽어들여서 총 sequence 수와 bp를 출력하는 스크립트로서 khmer package에 포함되어 있다. 염기서열 파일을 다루는 많은 유틸리티가 있지만 관성이란 참으로 무서운 것이라서 초기에 접하여 손에 익은 프로그램은 계속 사용하게 된다. 기본적인 사용법은 다음과 같다.
$ readstats.py -o readstats.txt --csv *fastq || This is the script readstats.py in khmer. || You are running khmer version 2.1.1 || You are also using screed version 1.0 || || If you use this script in a publication, please cite EACH of the following: || || * MR Crusoe et al., 2015. http://dx.doi.org/10.12688/f1000research.6924.1 || || Please see http://khmer.readthedocs.io/en/latest/citations.html for details. ... SRR13060952.fastq 0 ... found 37750 bps / 250 seqs; 151.0 average length -- SRR13060952.fastq ... SRR13060953.fastq 0 ... found 37750 bps / 250 seqs; 151.0 average length -- SRR13060953.fastq ... SRR13060954.fastq 0 ...중간 생략... $ head -n 5 readstats.txt bp,seqs,avg_len,filename 37750,250,151.0,SRR13060952.fastq 37750,250,151.0,SRR13060953.fastq 37750,250,151.0,SRR13060954.fastq 37750,250,151.0,SRR13060955.fastq
이틀 전에 소개한 HRGM(Human Reference Gut Microbiome, 블로그 링크)을 공부하기 위하여 NCBI Sequence Read Archive study SRP292575로 등록된 raw sequencing read를 다운로드하였다. 논문에서는 90개의 분변 시료라 하였는데 SRR accession으로는 106개였다.
SRA에서 파일을 다운로드할 때에도 parallel을 이용할 수 있다. 106개의 모든 SRR accession에 대하여 fasta-dump 명령어를 만들어서 다음과 같은 run.sh이라는 텍스트 파일을 만든다.
fastq-dump --split-files SRR13060942 fastq-dump --split-files SRR13060943 fastq-dump --split-files SRR13060944
이를 parallel에 공급하면 된다. Parallel을 사용하는 아주 기본적인 기법에 들어간다.
$ parallel --jobs 16 < run.sh & [1] 17592
모든 파일을 압축하고 싶다면 'parallel gzip ::: *fastq'를 실행하면 된다.
FASTQ 파일을 다운로드하였으니 각 파일이 몇 bp나 되는지를 헤아려 볼 순간이 되었다. 여러 파일의 다운로드 또는 압축 작업은 parallel을 쓰기에 아주 좋다. 왜냐하면 각 파일에 대해서 이루어진 작업에 대하여 출력되는 메시지 따위에 신경을 쓸 필요가 전혀 없기 때문이다. 그러나 readstats.py는 약간 다르다. 인수로 주어지는 FASTA/FASTQ 파일에 대한 집계 수치를 표준 출력/표준에러/'-o filename'으로 지정된 파일에 기록하게 되므로, parallel로 작업을 하려면 세심하게 명령어를 날려야 한다. 예를 들자면 Parallel tutorial 문서를 참조하여 다음과 같이 실행할 수 있다.
$ parallel --results outdir readstats.py ::: *fastq
그러면 새로 만들어진 outdir 디렉토리 아래에 각 fastq 파일에 대한 별도의 디렉토리가 생기고, 그 아래에 seq stderr stdout이라는 파일이 생긴다. 솔직히 말하자면 이것을 정리하는 일이 더욱 귀찮다. 따라서 입력 파일 하나에 대해서 '-o filename'으로 출력 파일을 하나씩 만든 후(이때 parallel을 사용), 가장 마지막 단계에서 모든 출력파일을 집계하여 최종 결과를 만드는 것이 더욱 바람직하다. Parallel 명령어는 다음과 같이 만들어서 실행하라.
$ parallel -j 16 'readstats.py {} > {.}.stat' ::: *fastq
{}는 input line을, {.}는 input line에서 파일 확장자를 제거했음을 의미한다. 매뉴얼을 잘 읽어보면 의미를 알 수 있지만 사용법을 빨리 익히려면 실제 사례를 보는 것이 낫다. 오늘 익힌 이 사례는 Parallelising jobs with GNU Parallel이라는 글에서 참고하였다. 이 글을 쓴 사람은 유전체학 분야에서 일을 하는 것으로 여겨진다. 이 글이 실린 RONIN - our guide to cloud computing's tricky bits이라는 곳이 무슨 목적으로 만들어졌는지도 관심을 가져 보아야 되겠다. RONIN이란, AWS(아마존웹서비스) 위에서 돌아가는 연구용 클라우드 컴퓨팅 서비스 정도로 이해하면 될 것 같다.
댓글 없음:
댓글 쓰기