NCBI에서 미생물(bacteria) 유전체 서열을 한번에 다운로드하기
지금 와서 돌이켜보니 설명이 많이 부족하였다. 또한 이미 받은 서열 파일 이외에 NCBI에서 추가로 등재된 파일만을 가져오려면 어떻게 하는 것이 좋을까? 지난 4월 초에 총 6만개가 넘는 유전체 서열 파일(*_genomic.fna.gz)를 받는데 거의 5일 가까이 걸렸었다. 만약 ftp 세션을 여러개 열어서 나름대로 최적화를 시도했다면 조금 더 빨리 받을 수도 있었겠지만 ftp 서버의 관리자가 이를 좋아하지는 않았을 것이다. 데이터 전체를 매번 새로 받을 것이 아니라 새로 등록된 것만 받고, 아울러서 무슨 이유인지는 몰라도 삭제된 것은 제거하는 방식으로 데이터베이스를 유지하는 것이 바람직할 것이다. 오늘의 포스팅에서는 그 실제적인 방법을 정리하고자 한다.
박테리아의 경우 NCBI RefSeq ftp 사이트에 있는 파일 assembly_summary.txt가 가장 핵심적인 파일이다. 박테리아가 아닌 곰팡이나 식물 등 다른 계열의 것을 원한다면 다음 위치에서 해당 디렉토리로 들어가서 assembly_summary.txt를 찾아라.
ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/
이 파일을 이루는 각 컬럼의 설명이 필요하면 여기를 확인해 보라. assembly_summary.txt의 20번째 필드에 ftp를 위한 기본 path가 수록되어 있다. 이를 awk로 뽑아내어 별도의 파일 bacteria_dir.txt를 만든다. 시험삼아 Paenibacillus polymyxa E681 assembly version 2의 RefSeq ftp 기본 위치를 가 보자. 10여 개의 파일이 있다. 나는 이 중에서 유전체 염기서열 파일(*_genomic.fna.gz)만을 받을 것이다. 그러려면 bacteria_dir.txt 파일을 적절히 가공하여 유전체 염기서열 파일의 전체 다운로드 경로 파일(bacteria_genomic_file)을 만든 뒤 wget으로 받으면 된다. 전체 과정을 하나의 스크립트로 표현하면 다음과 같다.
$ curl 'ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt' > bacteria_assembly_summary.txt
$ awk '{FS="\t"} !/^#/{print $20}' bacteria_assembly_summary.txt > bacteria_dir.txt
$ sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/)(GCF_.+)|\1\2/\2_genomic.fna.gz|' bacteria_dir.txt > bacteria_genomic_file
$ wget --input bacteria_genomic_file
첫번째 명령은 wget ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt라고 단순하게 입력해도 된다.
RefSeq가 아니라 GenBank 레코드를 받고 싶다면 파일 경로와 이름이 약간 달라진다. 바뀌어야 하는 곳을 위에서 빨강색으로 표시했으니 이는 숙제로 남긴다. 그러면 약간의 응용을 해 보자. 만약 genus Vibrio의 완성된 유전체에 해당하는 것만을 전체 다운로드하려면 어떻게 하면 좋을까? assembly_summary.txt에서 해당되는 조건의 row만을 뽑아서 좀 더 특화된 다운로드 경로 파일을 만들면 된다. 위에 소개한 명령 세 줄에서 가운데 것만 다음과 같이 바꾸면 된다. awk에서는 따옴표로 둘러친 명령어 내에서 field separator를 탭("\n")으로 지정해도 되고, 명령행 옵션으로(awk -F "\t") 지정해도 된다.
$ awk '{FS="\t"} $8 ~ /^Vibrio/ && $12 ~ /Complete Genome/ {print $20}' bacteria_assembly_summary_05-24-2016.txt > bacteria_dir.txt
Incremental download 요령
우선 로컬 파일시스템에 다운로드한 기존의 파일 이름을 수록한 목록 파일을 만든 뒤 sort를 해 둔다.
$ ls | grep _genomic.fna.gz | sort > DOWNLOADED_FILES_SORTED
다음으로는 최산 assembly_summary file에서 파일 목록을 만들어 둔다. 다운로드 링크를 붙이지 않은 파일 이름만의 목록을 만들어야 하고, 마찬가지로 sort도 거쳐야 한다.
$ sed -r 's|(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/)(GCF_.+)|\2_genomic.fna.gz|' bacteria_dir.txt | sort > ALL_FILES_SORTED
이제 comm이라는 놀라운 유틸리티를 구동하자. comm은 두 개의 파일을 라인 단위로 비교하여 첫번쨰 파일에만 있는 라인, 두번쨰 파일에만 있는 라인, 그리고 모든 파일에 공통으로 있는 라인을 각각 별도의 컬럼으로 출력한다. 두 파일은 전부 사전에 sort가 되어있어야 하고, 같은 내용의 라인을 중복해서 갖고 있어서는 안된다. 인수(-1, -2, -3)는 출력이 되지 않도록 하는 컬럼을 지정하는 것이다. 즉 comm -23 file_1 file_2은 첫번째 파일에만 있는 라인을, comm -13 file_1 file_2는 두번째 파일에만 있는 라인을 출력한다.
$ comm -23 DOWNLOADED_FILES_SORTED ALL_FILES_SORTED > FILES_TO_DELETE$ comm -13 DOWNLOADED_FILES_SORTED ALL_FILES_SORTED > FILES_TO_DOWNLOAD
$ sed -r 's|(^.*)_genomic.fna.gz$|ftp://ftp.ncbi.nlm.nih.gov/genomes/all/\1\/\1_genomic.fna.gz|' FILES_TO_DOWNLOAD > FILES_TO_DOWNLOAD_WITH_LINKS
$ wget --input FILES_TO_DOWNLOAD_WITH_LINKS
$ cat FILES_TO_DELETE | while read f
> do
> rm $f
> done
파일 이름에 생성 날짜를 넣고 싶다면
$ mv bacteria_assembly_summary.txt bacteria_assembly_summary_`date +%m-%d-%Y`.txt
펄 스크립트로 무엇인가를 구현하려고 머리를 쥐어짜다가 조금만 검색을 거치면 이미 그런 용도로 만들어진 리눅스 유틸리티를 발견하고 늘 감탄을 하게된다. 세상은 넓은데 내가 새로이 할 일은 별로 없다!
댓글 없음:
댓글 쓰기