Genome Informatics Laboratory at KRIBB

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== E-direct로 명령행에서 Entrez DB 접근하기 ======

===== 기본 개념 및 설치 방법 =====

Entrez는 미국 NCBI가 보유한 biomedical database의 통합적 질의 시스템 혹은 웹 포털을 의미한다. NCBI에서는 유닉스 혹은 리눅스의 명령행 환경에서 Entrez에 직접 접속하여 검색을 실시하고 그 결과를 활용하는데 필요한 유틸리티를 제공하고 있는데, 이것이 바로 [[https://www.ncbi.nlm.nih.gov/books/NBK179288/|Entrez Direct(E-Direct)]]이다. URL 구문을 이용하여 Entrez 시스템의 활용을 도와주는 API의 일종인 [[https://www.ncbi.nlm.nih.gov/books/NBK25501/|Entrez Programming Utilities(E-utilities)]]의 명령행 버전 정도로 이해하면 된다. Entrez의 모든 데이터베이스를 알아보려면  https://eutils.ncbi.nlm.nih.gov/entrez/eutils/einfo.fcgi를 접속해 보라.

웹사이트에서 설명한 설치방법 외에도 bioconda(entrez-direct) 또는 데비안 패키지(ncbi-entrez-direct)로도 제공되므로 이를 활용할 수 있다.

설치와 간단한 사용법에 대한 동영상은 유튜브('[[https://youtu.be/rPDXGsLcTQ0|How to Install EDirect on Mac, Unix or Linux]]')에서 볼 수 있다. E-Direct는 사용자가 제시하는 질의에 대하여 **esearch**로 탐색을 실시하고, **efetch**로는 찾아낸 레코드를 원하는 포맷에 맞게 가져오며, **xtract**를 사용하여 프로그래밍을 하지 않고도 전단계의 명령어에서 반환되는 정보로부터 필요한 필드만을 추출하는 다단계로 이루어진다. 기본적으로 표준 입출력을 사용하므로 각 명령어는 파이프(|)로 연결하여 사용할 수 있다. [[https://ncbi-hackathons.github.io/EDirectCookbook/|EDirect Cookbook]]에서는 직접 복사하여 쓸 수 있는 풍부한 예제가 실려 있으며, 특히 [[https://github.com/NCBI-Hackathons/ncbi-hackathons.github.io/wiki/EDirect-Wiki!|EDirect Wiki!]]에도 참조할만한 정보가 많다. CVR(MRC-University of Glasgow Centre for Virus Research)에서도 [[https://bioinformatics.cvr.ac.uk/ncbi-entrez-direct-unix-e-utilities/|NCBI Enterz Direct UNIX E-utilities]]라는 웹문서를 통해 몇 가지의 예제를 소개하였다. 

===== 유전체 정보 다운로드를 위한 FTP path 알아내기 =====
다음은 EDirect에서 //Clostridium botulinum//의 complete genome 정보를 다운로드할 수 있는 FTP path를 파일로 저장하는 명령어 사례이다. esearch -query '( … )'로 지정해야 하는 질의어가 예상외로 복잡하게 느껴질 수 있으므로 [[https://www.ncbi.nlm.nih.gov/assembly/advanced/|Assembly Advanced Search Builder]]를 이용하여 점검하는 것을 강력하게 권장한다. 이 웹사이트는 NCBI의 assembly DB에 대한 것이므로, 다른 종류의 DB(bioproject, biosample, pubmed, sra, taxonomy… 등)에 대해서는 별도의 Advanced Search Builder를 활용하기 바란다. 다음의 명령어 사례는 //Clostridium botulinum//의 complete genome 서열을 다운로드하기 위한 URL을 full_download_path.txt에 저장하기 위한 것이다. 질의어를 Advanced Search Builder에 넣으면 자동적으로 'derived from surveillance project'와 'anomalous' 레코드를 제외하는 필터가 적용된다. 후속 분석 목적에 따라서는 surveillance project에서 유래한 유전체 정보를 다운로드할 필요도 있으므로, 이러한 상황에서는 질의어를 어떻게 만들어야 하는지 생각해 보라. EDirect Cookbook에서는 질의어를 "species[ORGN] AND latest[SB]"로 매우 단순하게 구성하였다. 추출해야 할 필드가 여러 개라면 xtract의 인수를 '-element GenBank FtpPath_GenBank BioSampleAccn BioprojectAccn'와 같이 만들면 된다. 'efetch -format docsum' 명령어를 제공하면 XML 형태의 자료가 표준출력으로 나온다.

  $ esearch -db assembly -query '("Clostridium botulinum"[Organism] AND \
  (latest[filter] AND (all[filter] NOT "derived from surveillance \
  project"[filter] AND all[filter] NOT anomalous[filter]))")' | \
  efetch -format docsum | \
  xtract -pattern DocumentSummary -element FtpPath_RefSeq | \
  awk -F"/" '{print $0"/"$NF"_genomic.fna.gz"}' > full_download_path.txt

맨 마지막의 awk 명령어는 설명이 필요하다. 직전의 xtract 명령으로 출력되는 라인은 다음과 같은 형식을 갖는다. 이 URL은 파일이 있는 디렉토리이다.
  (1) ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/017/045/GCF_000017045.1_ASM1704v1 

그러나 다운로드할 수 있는 실제 파일의 URL은 다음과 같다.
  (2) ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/017/045/GCF_000017045.1_ASM1704v1/GCF_000017045.1_ASM1704v1_genomic.fna.gz

입력된 값인 (1)에서 슬래시로 구분된 가장 마지막 컬럼을 한번 더 반복한 다음 _genomic.fna.gz를 덧붙인 것이다. 이를 만들어 내기 위해 awk를 사용하였다. 컬럼 구분자는 기본값인 탭이 아니라 슬래시를 쓰게 하였고, [전체 라인 + 슬래시 + 마지막 필드 + _genomic.fna.gz를 출력하여 파일로 리다이렉트한 것이다. $0, $1, $2…는 awk에서 각 컬럼을 나타내는 변수이다. $1은 첫번째,  $2는 두번째 컬럼을 의미하며 $0은 라인 전체를 뜻한다. 그러면 마지막 컬럼은 어떻게 알아낼 것인가? 컬럼의 총 수를 미리 알기는 어려우니 $ 뒤에 구체적인 숫자를 넣을 수는 없지만, NF가 컬럼의 수를 의미하는 변수라는 점을 이용하여 $NF라 표시하면 이는 마지막 컬럼이 된다. 중간 파일(full_download_path.txt)에 기록하는 것을 거치치 않고 곧바로 wget으로 파이프 처리를 하여 다운로드를 해도 된다.

===== 미생물 대표 유전체 염기서열의 일괄 다운로드 =====
이상에서 학습한 것을 응용하여 RefSeq으로부터 prokaryotic reference & representative genome의 GenBank flat file을 한꺼번에 다운로드하는 방법을 알아본다. NCBI Insights에 이에 대한 업데이트 소식이 가끔 올라온다. [[https://www.ncbi.nlm.nih.gov/assembly/?term=(Bacteria%5Borgn%5D%20OR%20Archaea%5Borgn%5D)%20AND%20(reference_genome%5Bfilter%5D%20OR%20representative_genome%5Bfilter%5D)&utm_source=blog&utm_medium=referral&utm_campaign=Prok-rep-genomes-update&utm_term=rep-genome-collection&utm_content=20210513link1|여기]]를 클릭하면 오늘 기준의 미생물 참조 및 대표 유전체 서열 검색 결과가 보일 것이다. 웹페이지 내에서 Donwload Assemblies를 클릭하여 17000개가 넘는 파일을 순조롭게 내려받을 수 있다고 기대할 수는 없을 것이다. 구체적으로 어떤 검색어가 쓰였는지, 작동된 필터는 무엇인지 확인해 보는 것도 좋다.

  ("Bacteria"[Organism] OR "Archaea"[Organism]) AND (reference_genome[filter] OR representative_genome[filter])

검색결과 맨 위에 'Filters activated: Latest, Exclude anomalous'라는 메시지가 보일 것이다. 이 두 개의 필터는 자동으로 적용된 것 같다. 이 검색식을 사용하여 esearch – efetch – xtract 명령어를 만들어 본다. 앞서 보인 사례에서는 FtpPath_RefSeq element를 추출한 뒤 완전한 파일의 명칭까지를 재구성하였다. 그러나 이번에는 FtpPath_RefSeq만을 우선 추출하여 별도의 파일인 download_path.txt로 저장하도록 한다. 왜냐하면 이 디렉토리에 있는 md5checksums.txt 파일을 가져오려면 FtpPath_RefSeq 정보를 그대로 이용하는 것이 더 편리하기 때문이다. 혹은 위에서 만든 full_download_path.txt(full path를 포함한 파일 이름)를 적당히 조작하여 슬래시로 구분된 가장 마지막 필드를 제거해도 된다.

  $ esearch -db assembly -query '(Bacteria[orgn] OR Archaea[orgn]) AND \
      (reference_genome[filter] OR representative_genome[filter])' | \
      efetch -format docsum | xtract -pattern DocumentSummary \
      -element FtpPath_RefSeq > download_path.txt

[[:custom_scripts?s[]=safe&s[]=download#safe_downloadsh|safe_download.sh]] 스크립트에 download_path.txt 파일을 인수로 제공하면  md5checksums.txt와 GenBank flat file을 같이 다운로드한 다음, md5 checksum을 점검하여 만약 오류가 있으면 별도의 파일에 이를 기록해 놓는다. 정상적으로 다운로드된 파일은 미리 만들어진 download 디렉토리로 옮겨진다. 따라서 스크립트가 종료한 뒤 다운로드에 실패한 파일만 별도로 다시 받으면 된다.

  $ safe_download.sh download_path.txt
  

===== 다른 정보를 추출하기 =====

만약 FtpPath뿐만 아니라 여러 가지의 정보를 추출하고 싶다면 '-element Genbank FtpPath_GenBank BioSampleAccn BioprojectAccn((RefSeq에도 등재된 레코드는 두 개의 BioprojectAccn이 출력된다. 두 번째 것은 RefSeq Prokaryotic Genome Annotation Project(PRJNA224116)로서 모든 레코드에 대하여 동일하다.))'와 같이 xtract 명령어의 옵션을 확장하면 된다. 어떤 종류의 element가 있는지 알고 싶다면 xtract 명령 없이 'esearch … | efetch -format docsum' 명령으로 XML 형태의 자료를 출력하여 <DocumentSummary>..</DocumentSummary> 사이의 블록을 살펴보도록 한다.

  $ esearch -db assembly -query GCA_000017025.1
  <ENTREZ_DIRECT>
    <Db>assembly</Db>
    <WebEnv>MCID_5f582c1bd7eb634c162e77b8</WebEnv>
    <QueryKey>1</QueryKey>
    <Count>1</Count>
    <Step>1</Step>
  </ENTREZ_DIRECT>
  $ esearch -db assembly -query GCA_000017025.1 | efetch -format docsum
  <?xml version="1.0" encoding="UTF-8" ?>
  <!DOCTYPE DocumentSummarySet PUBLIC "-//NLM//DTD esummary assembly 20200615//EN" 
  "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20200615/esummary_assembly.dtd">
  
  <DocumentSummarySet status="OK">
  <DocumentSummary>
  <Id>38128</Id>
  	<RsUid>38128</RsUid>
  … 중간 생략 …
  </DocumentSummary>
  </DocumentSummarySet>

이번에는 유전체 정보가 공개된 특정 종 균주가 어떤 곳에서 분리되었는지 알아내는 방법에 대해여 설명하고자 한다. 위에서 설명한 xtract 명령어에서 ‘-element BioSampleAccn’ 옵션을 적용하여 그 결과를 biosample_acc.txt 파일에 저장한다. 이를 esearch + efetch + xtract 명령어에 적용하여 BioSample 데이터베이스를 검색하되 모든 Accession Number에 대해서 적용하도록 반복문을 돌리는 것이 핵심이다.

  $ cat biosample_acc.txt
  SAMN13428652
  SAMN11521245
  SAMN11521242
  SAMN11521244
  SAMN10822424
  $ for i in $(cat biosample_acc.txt)
  > do
  > ll=$(esearch -db biosample -query $i | efetch -format docsum | \
  > xtract -pattern DocumentSummary -block Attribute \
  > -if Attribute@attribute_name -equals isolation_source \
  > -element Attribute)
  > echo -e "$i\t$ll"
  > done
  SAMN13428652    Carrot juice
  SAMN11521245    missing
  SAMN11521242    missing
  SAMN11521244    Homo sapiens
  SAMN10822424    missing
  …(후략)

BioSample 데이터베이스의 Attribute 항목에는 여기에서 소개한 isolation_source 외에도 strain, collection_date, geo_loc_name, sample_type, host_disease 등이 있으니 용도에 맞게 적절히 활용하면 된다. 그러나 BioSample 데이터베이스에 수록된 정보는 자료마다 일정하지가 않아서 특정 attribute의 것만을 취하면 중요한 정보를 놓치는 경우가 많다. 따라서 모든 자료를 (key, value) 쌍의 형태로 추출하는 것이 바람직하다. Assembly database에서 균주와 관련한 정보를 출력하여 이를 data.txt 파일에 저장해 놓았고, assembly accession과 BioSample accession이 이중에서 첫 번째와 세 번째 컬럼에 있다면 다음 스크립트를 이용해 본다. 단, 출력물은 처음 두 컬럼을 제외하면 순서가 일정하지 않으므로 Perl 스크립트 등을 이용하여 후처리를 해야 된다.

  # data.txt 파일은 ‘xtract –pattern DocumentSummary -element Genbank FtpPath_GenBank BioSampleAccn BioprojectAccn’으로 추출된 것으로, 첫 번째와 세 번째 컬럼(굵은 글씨)을 사용함
  $ head -n 1 data.txt
  GCA_014068615.1 ftp://ftp.ncbi.nlm.nih.gov/..(생략).._ASM1406861v1 SAMN13428652     
  PRJNA592543     PRJNA224116
  $ for i in $(awk -vOFS=, '{print $1, $3}' data.txt)
  > do
  > A=${i%,*}
  > B=${i#*,}
  > ll=$(esearch -db biosample -query $B | efetch -format docsum | \
         xtract -pattern DocumentSummary -block Attribute \
         -element Attribute@attribute_name Attribute)
  > echo -e "$A\t$B\t$ll" 
  > done
  GCA_014068615.1 SAMN13428652    strain  CJ0611A1        isolation_source       Carrot juice     collected_by    Toronto Public Health   collection_date 2006-10-18      geo_loc_name    Canada: Toronto lat_lon 45.406270 N 75.743210 W
  …(후략)

위의 사례에서는 data.txt 파일에서 유래한 두 컬럼을 do 반복문 안에서 A와 B 변수로 나누어 저장하기 위해 다소 복잡한 방법1을 동원하였다. 그러나 'awk '{print $1, $2}' data.txt | while read A B' 명령으로는 첫 번째 레코드만 처리되고 종료하는 문제점이 발생하였다. 그리고 awk 문에서 출력용 필드 구분자를 콤마로 바꾸지 않고 기본값(공백 문자) 그대로 사용하면 do 반복문 내부의 실행이 정상적으로 이루어지지 않는다. 만약 이보다 더 합리적인 방법을 찾아낼 수 있다면 바람직할 것이다.

다음의 사례는 출력물의 세 번째 컬럼에 Organism 정보를 포함시켜서 출력하되 모든 결과를 results.txt 파일에 기록하도록 하는 방법이다. BioSample 레코드에서는 Organism이 별도의 block에 들어있지는 않지만, 복수의 -block 구문을 사용하여 서로 다른 섹션에 위치한 정보를 하나의 xtract문으로 추출하는 것이 가능하다([[https://www.ncbi.nlm.nih.gov/books/NBK179288/|NCBI Bookshelf]] 'Exploring Separate XML Regions' 참조). 

  $ for i in $(awk -vOFS=, '{print $1, $3}' data.txt)
  > do
  > A=${i%,*}
  > B=${i#*,}
  > ll=$(esearch -db biosample -query $B | efetch -format docsum | \
         xtract -pattern DocumentSummary -element Organism \
         -block Attribute -element Attribute@attribute_name Attribute)
  > echo -e "$A\t$B\t$ll" >> results.txt
  > done