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.
====== 유전체 주석화(genome annotation) ======
주석화 결과를 외부와 공유하고 관련 유전체와 비교 분석을 할 목적이라면 [[https://rast.nmpdr.org/|RAST(Rapid Annotation using Subsystem Technology) server]]에 회원 등록을 하고 염기서열 파일을 업로드하는 것으로 충분하다. 그러나 로컬 컴퓨터에 직접 프로그램을 깔아서 유전체 주석화를 실시하고 싶은 욕구는 누구나 갖고 있을 것이다. 어떤 프로그램이 있는지 알아보자.
===== Prokka 사용하기 =====
명령행 환경에서 실행되는 미생물 유전체 주석화 프로그램으로 가장 대중적인 것은 [[https://github.com/tseemann/prokka|Prokka]]라 해도 과언이 아닐 것이다. 실제 실행 사례는 Prokka GitHub 페이지의 [[https://github.com/tseemann/prokka|Invoking Prokka]] 항목을 참고한다. Contig ID가 37 문자를 넘어가면 실행이 되지 않으니 CLC Genomics Assembly 등으로 조립한 경우 이를 적절한 길이로 줄여야 한다.

  $ prokka --outdir mydir --prefix mygenome --locustag AYCC --genus Escherichia --species coli --strain mystrain --cpus 16 genome.fasta
  
===== DFAST 사용하기 =====
일본에서 개발한 DFAST(DDBJ Fast Annotation and Submission)의 stand-alone version인 [[https://github.com/nigyta/dfast_core|dfast_core]]는 functional annotation을 위해 사용하는 데이터베이스 용량이 prokka보다는 좀 더 크고, pseudogene에 대한 정보를 주기 때문에 유용하다. DFAST는 [[https://dfast.ddbj.nig.ac.jp/|웹 버전]]으로도 사용 가능하다. 

===== PGAP 사용하기 =====
Prokaryotic Genome Annotation System(PGAP, [[https://www.ncbi.nlm.nih.gov/genome/annotation_prok/|NCBI]] or [[https://github.com/ncbi/pgap|GitHub]])은 세균 유전체의 자동 주석화를 위하여 NCBI에서 공식적으로 사용하는 프로그램이다. 여러 개 유전체 서열에 대하여 신속하게 주석화를 하려면 Prokka가 매우 편리하지만, 대용량의 DB를 참조하여 주석화를 실시하는 PGAP이 더욱 양질의 결과를 산출하게 된다. 원래 PGAP은 RefSeq genome의 주석화용으로 내부적으로만 쓰이다가 누구나 설치할 수 있는 형태로 배포되기에 이르렀다. 유튜브에는 사용자의 유전체를 PGAP으로 직접 주석화하는 방법을 소개하는 [[https://youtu.be/pNn_-_46lpI|동영상]]이 올라와 있다. 설치와 사용 방법에 대한 상세한 설명은 PGAP 위키 사이트의 [[https://github.com/ncbi/pgap/wiki/Quick-Start|Quick-Start]]를 참고하도록 한다. Standalone 버전이 처음 나왔을 떄에 비하면 설치 방법이 훨씬 간단해진 것 같다. PGAP 버전 번호는 ‘YYYY-MM-DD.build####’의 형식을 따른다. 2023년 6월 22일에 설치한 input-2023-05-17.build6771 버전의 설치 후 용량은 32GB 정도이다. PGAP은 conda와 상관이 없다.

  # 현재 배포 중인 PGAP의 최신 버전 확인하기
  $ curl --silent "https://api.github.com/repos/ncbi/pgap/releases/latest" | grep -Po '"tag_name": "\K.*?(?=")' > VERSION
  $ cat VERSION 
  2024-07-18.build7555

PGAP은 docker 환경을 사용하므로, 사용자는 관리자이거나 sudo 권한을 갖고 있어야 한다. PGAP 배포판에 포함된 샘플 유전체 서열을 대상으로 주석화를 실행하는 방법은 다음과 같다. pgap.py 스크립트는 /data/apps/pagp에 있다고 가정한다.

  # docker가 실행 중인지 확인
  $ ps aux | grep -i docker | grep -v grep # 또는 ps –ef | grep docker
  root     25146  0.8  0.0 2681660 107976 ?      Ssl  01:13   2:33 /usr/bin/dockerd -g /data/docker -H fd:// --containerd=/run/containerd/containerd.sock
  # PGAP 설치
  $ cd /data/apps/pgap
  $ curl -OL https://github.com/ncbi/pgap/raw/prod/scripts/pgap.py
  $ chmod +x pgap.py
  $ ./pgap.py --update
  The latest version of PGAP is 2023-05-17.build6771, you have nothing installed locally.
  installation directory: /home/hyjeong/.pgap
  /bin/df: /home/hyjeong/.pgap: No such file or directory
  Downloading (as needed) Docker image ncbi/pgap:2023-05-17.build6771
  Downloading and extracting tarball: https://s3.amazonaws.com/pgap/input-2023-05-17.build6771.tgz
  Installing PGAP test genomes
  /home/hyjeong/.pgap/test_genomes-2023-05-17.build6771
  https://s3.amazonaws.com/pgap-data/test_genomes-2023-05-17.build6771.tgz
  Downloading and extracting tarball: https://s3.amazonaws.com/pgap-data/test_genomes-2023-05-17.build6771.tgz
  ...
  Status: Downloaded newer image for ncbi/pgap:2023-05-17.build6771
  docker.io/ncbi/pgap:2023-05-17.build6771
  # test genome(L43967.2 Mycoplasma genitalium G37, complete genome, 580076 bp) 주석화
  $ cp  $HOME/.pgap/test_genomes/MG37/ASM2732v1.annotation.nucleotide.1.fasta .
  $ ./pgap.py -r -o mg37_results -g ASM2732v1.annotation.nucleotide.1.fasta -s 'Mycoplasmoides genitalium'
  PGAP version 2023-05-17.build6771 is up to date.
  Output will be placed in: /data/apps/pgap/mg37_results
  ...

PGAP은 $HOME/.pgap에 설치되며, pgap.py 스크립트의 실행 위치는 별로 중요하지 않다. -r 또는 %%--%%report-usage-true 옵션은 PGAP을 실행할 때 NCBI에 보고를 하기 위한 것이다. 이를 원치 않으면  -n 또는 %%--%%report-usage-false 옵션을 지정하면 된다. 

==== GenBank 제출을 위한 실행(선택) ====

PGAP에 필요한 [[https://github.com/ncbi/pgap/wiki/Input-Files|input file]]은 사용자가 준비한 실제 유전체 서열(FASTA)과 두 개의 YAML 형식 메타데이터 파일(generic and submol)이다. Generic YAML 파일은 pgap.py 실행에서 직접적인 인수로 주어진다(위에 보인 test genome의 PGAP 실행에 대해서는 input.yaml 파일). 만약 유전체 서열 파일이 Ecoli1_genomic.fna라 하면, generic YAML file의 내용은 다음과 같다. YAML 파일에서는 탭문자가 아니라 공백을 써야 함에 유의하라. 

  report_usage: true 
  fasta: 
      class: File
      location: Ecoli1_genomic.fna
  submol:
    class: File
    location: E_coli1.yaml

FASTA file의 서열 ID는 너무 길게 작성하지 않도록 한다. 이것이 결과물(GenBank file)에서 첫 줄의 LOCUS 값으로 쓰이게 되므로, 서열 ID가 너무 길면 다음과 같이 잘리는 불상사가 발생한다.

  $ head -n 1 pilon-3rd.fasta
  >Akkermansia_muciniphila_strain_Pendulum
  # PGAP 실행
  $ head –n 4 annot.gbk
  LOCUS       Akkermansia_muciniphila_strain_Pendul> 2664055 bp    DNA
              circular BCT 22-APR-2021
  DEFINITION  Akkermansia muciniphila strain Pendulum chromosome, complete
              genome.

Submol file의 명칭은 generic YAML file에서 지정한다. 필수 필드가 아닌 것이 상당히 많으므로 위키 페이지를 참조하여 적절히 취사 선택하기 바란다. 혹은 test_genomes/MG37 디렉토리에 있는 샘플용 input.yaml과 submol.yml 파일을 가져다가 필요한 부분만 수정하여 사용하는 것도 좋을 것이다. 가장 간단하게는 organism의 genus_species 정보만 있어도 실행이 된다. 다음은 submol YAML file의 예제이다. 이 파일의 bioproject/biosample/sra accession은 유효한 것이 아니므로 이를 수정 없이 그대로 사용하면 에러가 발생한다.

  topology: 'circular'
  organism:
      genus_species: 'Escherichia coli'
      strain: 'my_strain'
  contact_info:
      last_name: 'Doe'
      first_name: 'Jane'
      email: 'jane_doe@gmail.com'
      organization: 'NIH'
      department: 'NCBI'
      phone: '301-555-0245'
      fax: '301-555-1234'
      street: '9000 Rockville Pike'
      city: 'Bethesda'
      state: 'MD'
      postal_code: '20850'
      country: 'USA'
  authors:
      - author:  
          last_name: 'Doe'    
          first_name: 'Jane'
          middle_initial: 'A'
      - author:  
          last_name: 'Doe'    
          first_name: 'John'
  consortium: 'E. coli genome group'
  bioproject: 'PRJ9999999'
  biosample: 'SAMN99999999'      
  locus_tag_prefix: 'pgaptmp'
  sra:
      - accession: 'SRR9999999'
      - accession: 'ERR9999999'
  publications:
      - publication:
          pmid: 29112715

Contact_info를 이루는 필드 중 전화번호 등 어느 하나라도 임의로 넣지 않으면 PGAP 실행 중에 에러가 발생한다. Contact_info를 전부 생략하는 것은 상관이 없다. Locus_tag_prefix를 지정하지 않으면 ‘pgaptemp’가 기본값으로 쓰인다. submol YAML에서 정의한 topology(linear or circular)는 FASTA file에 여러 염기서열이 존재할 경우 모두에게 똑같이 적용된다. 따라서 염기서열에 따라서 topology를 다르게 지정하려면 FASTA file의 definition line에 ‘>seq1 [topology=circular]’의 형태로 기재하면 된다. Topology 정보가 없으면 linear로 간주한다.

==== 개별 사용을 위하여 PGAP을 실행하기 ====
GenBank에 제출할 용도가 아니라 개인적으로 쓰기 위한 것이라면, 부수적인 정보를 담고 있는 YAML 파일은 필요하지 않다. 단지 genome FASTA file과 미생물의 이름(organism name)만 있으면 된다. Organism name은 'genus' 또는 'genus species' 형식으로 제공한다. 이 기능은 비교적 최근에 추가되었다.
  $ ./pgap.py -r -o <results> -g <fasta> -s '<organism_name>'



===== GenBank flat file에서 정보 추출하기 =====
많은 유전체 주석화 도구가 GenBank 파일([[https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html|샘플 레코드 및 설명]]) 형태로 결과물을 제공한다. GUI JAVA 프로그램인 [[https://www.sanger.ac.uk/tool/artemis/|Artemis]]를 사용하면 유전체 염기서열 문맥 안에서 [[https://www.insdc.org/submitting-standards/feature-table/|feature]] 또는 염기 단위의 다양한 분석과 시각화 작업을 할 수 있다. 염기서열 및 주석화 정보의 편집도 artemis 내에서 자유롭게 할 수 있다. 때로는 후속 작업을 위해 GenBank 파일에서 locus tag, protein, amino acid sequence) 등의 형태를 tab-delimited file로 추출하는 것이 필요하지만 기성 프로그램은 염기/아미노산 서열의 일괄 출력이나 GFF([[https://asia.ensembl.org/info/website/upload/gff.html|GFF/GTF]] or [[https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md|version 3]]) export 기능 정도만 제공하는 것이 일반적이다. [[:custom_scripts#gbkInfo.pl|gbkInfo.pl]] 스크립트는 바로 이런 상황에서 사용하기 위해 만들어졌다. [[https://github.com/lskatz/lyve-SET/blob/master/scripts/getgbk|getgbk]]는 원래 VirtualBox용으로 배포된 [[https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0060120|CMG-Biotlools]]에 포함되어 있던 Perl script로서, accession number를 이용하여 GenBank 포맷 파일을 다운로드한다.

  $ gbkInfo.pl 
  Usage : gbkInfo.pl <GenBank file> [-seq]
          script last modified at Thu Jun 29 09:48:07 2023
  $ getgbk -a CP000727 > CP000727.gbk
  $ gbkInfo.pl CP000727.gbk 
  [STDERR] LOCUS: CP000727 (VERSION: 1)
  [STDERR] Organism: Clostridium botulinum A str. Hall
  [STDERR] CP000727.gbk has 1 sequence(s)
  [STDERR] CP000727.gbk has 3571 gene features
  [STDERR] CP000727.gbk has 3403 active CDS features (not marked as 'pseudo')
  [STDERR] CP000727.gbk has 56 pseudo genes
  [STDERR] Feature information is being written to [ CP000727.gbk.txt ] (pre-existing file was overwritten!)

INPUT.gbk 파일을 인수로 넣으면 INPUT.gbk.txt라는 파일이 만들어지는데, 여기에는 탭으로 구분된 16개의 컬럼이 출력된다. 컬럼의 이름은 출력 파일의 첫 줄(‘#’로 시작)에서 확인할 수 있다. 각 유전자의 염기서열 및 아미노산 서열을 출력하려면 스크립트 실행시 '-seq' 옵션을 맨 뒤에 주어야 한다. 원본 GenBank 파일에서 '/pseudo' qualifier를 지닌 유전자는 CDS feature를 동반하지 않으므로 product 정보도 존재하지 않는다. 따라서 이러한 유전자는 isPseudo? 컬럼이 pseudo로 표기되고 위치와 strand를 제외한 나머지 정보는 빈 상태가 된다. 

하나의 GenBank file에 여러 염기서열이 담겨 있어도 작동은 정상적으로 이루어지며, 출력 파일에는 각 locus(=sequence)에 따라 순서대로 결과가 정돈된다.

  $ gbkInfo.pl GCF_000063585.1_ASM6358v1_genomic.gbff
  [STDERR] LOCUS: NC_009495
  [STDERR] Organism: Clostridium botulinum A str. ATCC 3502
  [STDERR] LOCUS: NC_009496
  [STDERR] Organism: Clostridium botulinum A str. ATCC 3502
  [STDERR] GCF_000063585.1_ASM6358v1_genomic.gbff has 2 sequence(s)
  [STDERR] GCF_000063585.1_ASM6358v1_genomic.gbff has 3709 gene features
  [STDERR] GCF_000063585.1_ASM6358v1_genomic.gbff has 3630 active CDS features (not marked as 'pseudo')
  [STDERR] GCF_000063585.1_ASM6358v1_genomic.gbff has 75 pseudo genes
  [STDERR] Feature information is being written to [ GCF_000063585.1_ASM6358v1_genomic.gbff.txt ] (pre-existing file was overwritten!)

gbkInfo.pl 스크립트는 단순한 유전자 구조를 지닌 prokaryotic genome의 GenBank 파일에 대해서만 정상적인 작동을 보장한다. 따라서 spliced gene이나 translation exception 등의 정보가 수록된 GenBank 파일은 처리하지 않기를 권한다.
==== 특정 위치(1-bp lentgh)의 유전자, 염기, 코돈 및 아미노산 알아내기 ====
일루미나 시퀀싱 데이터를 이용하여 참조 서열을 기준으로 변이가 발생한 위치와 내역을 알아냈다 하여도 분석 프로그램에 따라서는 어느 유전자에 변이가 발생하였는지, 그리고 어느 코돈이 어떻게 바뀌어서 실제로 아미노산 잔기가 바뀌었는지를 친절하게 보여주지는 않는다. 여기에서는 바로 직전 섹션에서 설명한 gbkInfo.pl의 부가 기능을 사용하여 문제를 해결하는 방법을 알아보도록 하자. 우선 변이가 발생한 염기가 어느 유전자에 위치하는지를 확인하여 data.txt 파일에 기록해야 한다. 이를 위해서는 (1) GenBank 파일로부터 만든 GFF 파일과 변이 위치를 수록한 BED 파일(position.bed; length는 반드시 1 bp여야 함)을 만든 뒤 [[https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html|bedtools intersect]]에 인수로 공급하면 공통으로 존재하는 위치가 출력된다. 기본 조건에서는 (1)과 (2)의 컬럼이 전부 합쳐진 상태로 출력이 되므로, 서열 ID와 유전자의 locus tag 및 염기의 위치만 출력되도록 다소 복잡한 awk one-liner를 이용하였다.

  # 다음의 sed 명령은 GFF 파일에서 ‘##FASTA’로 시작하는 줄부터 파일 끝까지를 제거한다.
  $ bp_genbank2gff3.pl -r CP000727.gbk -o - | sed -n '/^##FASTA/q;p' > CP000727.gff
  # 변이 위치 BED 파일(position.bed)은 단일염기변이 기준으로 작성되어야 하므로 시작과 끝 위치는 동일하다. 또한 시작 위치를 기준으로 정렬이 되어 있어야 한다.
  $ cat position.bed 
  CP000727	758733	758733
  CP000727	867070	867070
  CP000727	1672839	1672839
  CP000727	2233512	2233512
  CP000727	2687719	2687719
  CP000727	3340681	3340681
  $ bedtools intersect -wa -wb -a CP000727.gff -b position.bed | awk -F"\t" -vOFS="\t" '$3~/gene/{p=$11; gsub(/ID=/,""); gsub(/;.*$/, ""); gsub(/\.gene*$/, ""); print $1, $9, p}' > data.txt
  $ cat data.txt 
  CP000727	CLC_0737	758733
  CP000727	CLC_0840	860707
  CP000727	CLC_1601	1672839
  CP000727	CLC_2104	2233512
  CP000727	CLC_2562	2687719
  CP000727	CLC_3207	3340681
  
data.txt 파일이 준비되었으면 gbkInfo.pl을 debug mode로 실행한다. gbkInfo.pl 스크립트를 적당한 편집기로 열어서 $debug 변수를 1로 바꾼다. INPUT.gbk.txt 파일은 정상적으로 만들어진다.

  $ ./gbkInfo.pl CP000727.gbk > results.txt
  [STDERR] Entering debug mode...[ $debug = 1 ; data.txt file will be used ]
  [STDERR] LOCUS: CP000727 (VERSION: 1)
  [STDERR] Organism: Clostridium botulinum A str. Hall
  [STDERR] CP000727.gbk has 1 sequence(s)
  [STDERR] CP000727.gbk has 3571 gene features
  [STDERR] CP000727.gbk has 3508 active CDS features (not marked as 'pseudo')
  [STDERR] CP000727.gbk has 56 pseudo genes
  [STDERR] Feature information is being written to [ CP000727.gbk.txt ] (pre-existing file was overwritten!)
  $ cat results.txt 
  #SeqID	locus_tag	strand	feature_start	SNP_position	NT_position (gene)	AA_position (gene)	position_in_codon	codon	AA	nucleotide (genome)
  CP000727	CLC_0737	+	757560	758733	1174	392	1	AGA	R	A
  CP000727	CLC_0847	+	866364	867070	707	236	2	GAG	E	A
  CP000727	CLC_1601	+	1672354	1672839	486	162	3	AAT	N	T
  CP000727	CLC_2104	+	NA	2233512	NA	NA	NA	NA	NA	C
  CP000727	CLC_2562	-	2687587	2687719	546	182	3	TTA	L	T
  CP000727	CLC_3207	-	3340611	3340681	1493	498	2	GTA	V	A

맨 마지막 컬럼에 보이는 염기는 유전체 서열을 기준으로 한 것이므로, 만약 해당되는 유전자가 reverse strand에 있다면 상보적인 것을 택해야 한다. 이상에서 소개한 방법은 오직 **단일 염기의 치환** 결과물만을 정확하게 다룰 수 있음에 유의하라.
===== InterProScan 활용하기 =====
핵산 혹은 단백질 서열의 주석화 도구로서 BLAST가 널리 사용되지만 이것은 local alignment를 이용하므로 기능 정보의 transfer를 위한 엄격한 기준점을 제시하기 힘들고 DB의 정확성에 크게 의존한다. InterPro는 단백질 패밀리의 데이터베이스 및 검색도구('classification of protein families')로서 [[https://www.ebi.ac.uk/interpro/search/sequence/|웹서버]]에 단백질 혹은 CDS 서열을 입력하여 분석을 실행할 수 있다. InterPro에서는 protein signature 및 family 정보를 이용하므로 훨씬 정확한 기능 예측 결과를 제공한다. Entry type(homologous superfamily, family, domain, repeat 및 site)에 대한 개요는 [[https://www.ebi.ac.uk/interpro/help/faqs/|FAQ 웹사이트]]를 참조하라. 유전체에서 얻어진 프로테옴 세트 전체에 대한 분석을 실시하려면 local server에 InterProScan 소프트웨어와 데이터베이스를 설치하여 이용하는 것이 바람직하다([[https://interproscan-docs.readthedocs.io/en/latest/UserDocs.html#installation-requirements|설치 요구 사항]] [[https://interproscan-docs.readthedocs.io/en/latest/HowToDownload.html|복사본 다운로드]] [[https://interproscan-docs.readthedocs.io/en/latest/HowToRun.html|실행 방법]]).

동시에 여러 CPU를 사용하려면 -cpu <CPU> 또는 --cpu <CPU> 옵션으로 원하는 수치를 입력하면 된다. 기본 조건은 interproscan 설치 디렉토리에 있는 interproscan.properties 파일의 number.of.embedded.workers 파라미터에 정의되어 있다.

  $ /usr/local/apps/interproscan-5.36-75.0/interproscan.sh –i INPUT.faa –f tsv

-f <OUTPUT-FORMATS> 옵션은 [[https://interproscan-docs.readthedocs.io/en/latest/OutputFormats.html?highlight=format#output-formats|결과 파일의 포맷]]을 지정한다. 단백질 서열의 경우 TSV, XML 및 GFF3 파일이 기본 생성된다. Analysis의 범위를 원하는 수준으로 한정하려면 -appl TIGRFAM,PANTHER,CDD,PRINTS와 같이 사용하려는 member database의 목록을 제공하면 된다. 


==== InterProScan 결과물의 시각화 ====
여러 유전체 서열로부터 추출한 단백질 서열을 대상으로 InterProScan을 실시한 뒤 특정 InterPro entry의 분포가 species에 따라서 어떠한 차이를 보이는지를 R(ggplot2)로 시각화하는 방법을 알아보자. 기본 가정은 각 균주에 대하여 GCA_000235785.2.faa([[https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/235/785/GCF_000235785.1_ASM23578v2/|ftp 사이트]]) 형태의 아미노산 서열 파일을 얻은 뒤 이를 InterProScan으로 처리하여 GCA_000235785.2.faa.tsv 형태의 결과 파일을 얻었다는 것이다. Assembly accession(‘GCA_000235785.2’)과 species 및 strain의 관계는 acc_species_strain 파일에 탭으로 구분하여 수록해 두도록 한다. 집계할 InterPro entry는 family 파일에 넣어둔다. 이상의 두 개 파일은 수작업으로 적절히 만들어야 한다. Strain은 EzBioCloud에서 각 genome에 할당한 것을 기준으로 하였다.

  $ cat acc_species_strain
  GCA_001742205.1	Lactobacillus fermentum	strain=NCC2970
  GCA_001748065.1	Pediococcus pentosaceus	strain=LP28
  GCA_001886915.1	Leuconostoc mesenteroides	strain=DRC1506
  …
  $ cat family
  IPR039697
  IPR014182
  IPR002347
  IPR012079
  IPR012394

각 strain의 아미노산 서열 파일에 대하여 얻은 InterProScan 결과(.tsv)을 다음의 test.sh로 처리하여 family file에서 정의한 InterPro entry에 맞는 유전자의 수를 집계한 뒤 .tsv.txt 파일에 저장한다.

  $ cat test.sh
  #!/bin/sh
  
  for i in $(cat family)
      do 
          grep $i $1 | awk '{print $1}' | uniq > $i
          wc -l $i | awk -vOFS="\t" '{print $2, $1}'
          cat $i >> all
          rm $i
      done
  sort all | uniq > all.nr
  mv all.nr all
  wc -l all | awk -vOFS="\t" '{print $2, $1}'
  rm all
  $ ls *faa.tsv | while read f; do sh ./test.sh $f > $f.txt; done
  $ cat GCA_000568875.1.faa.tsv.txt
  IPR039697	2
  IPR014182	0
  IPR002347	4
  IPR012079	1
  IPR012394	1
  all	7

이제부터의 작업은 R에서 실시한다. 모든 *.tsv.txt 파일을 모아서 하나로 병합한 뒤 results.txt 파일에 저장한다.

  rm(list=ls())
  list.filenames = list.files(pattern="*faa.tsv.txt")
  # create an empty list
  list.data=list()
  length(list.filenames)
  # create a loop to read in your data
  for (i in 1:length(list.filenames))
  {
  list.data[[i]] = read.table(list.filenames[i],sep="\t",header=F)
  colnames(list.data[[i]]) = c("markers",list.filenames[i])
  }
  
  # full outer join
  for (i in 1:(length(list.filenames)-1))
  {
  df = merge(x=list.data[[i]],y=list.data[[i+1]],by="markers",all=TRUE)
  list.data[[i+1]] = df
  }
  df = list.data[[length(list.filenames)]]
  df[is.na(df)] = 0
  rownames(df) = df[,1]
  df = df[,-1]
  colnames(df) = gsub(".faa.tsv.txt","",list.filenames)
  
  # 최종 확인
  dim(df)
  View(df)
  
  # accession을 species로 치환
  data = read.table("acc_species_strain",sep="\t",stringsAsFactors=F)
  key = data[,1]
  names(key) = data[,2]
  
  x = c()
  for (i in names(df)){
  temp = names(key)[key==i]
  x = append(x, temp)
  }
  
  df = rbind(df, x)
  rownames(df)[7] = "species"
  View(df)
  
  # row를 정렬하고 transpose
  sorted = c("species","IPR039697","IPR014182","IPR002347","IPR012079","IPR012394","all") 
  df = df[sorted,]
  df.2 = as.data.frame(t(df))
  df.2[,2:7] = sapply(df.2[,2:7],as.numeric)
  str(df.2)
  write.table(df.2,"results.txt",sep="\t",quote=F) 

results.txt 파일을 이용하여 다양한 변형 및 탐색을 실시해 보자. 이때 dplyr 패키지가 매우 유용하게 쓰인다. 각 species에 대하여 family에 정의된 InterPro entry에 해당하는 유전자의 최댓값은 얼마인가? 그것에 해당하는 균주 또는 accession은 무엇인가? 다음 R 코드에서 그 방법을 알아보자.

  # R 세션을 끝내고 다시 시작한다고 가정함
  # 혹은 rm(list=ls())를 실행하여 변수를 초기화한다.
  # reading from pre-existing file
  df = read.table("results.txt",sep="\t")
  # accession을 추출하여 첫번째 컬럼으로 삽입한다.
  accession = rownames(df)
  df = cbind(accession, df)
  write.table(df,"results_with_accession.txt",sep="\t",row.names=F,quote=F)
  
  # 각 species에 대하여 family에 해당하는 유전자가 가장 많이 검출된 균주(assembly accession; row)를 찾는다.
  # all 대신 IPR039697 등 다른 컬럼을 선택할 수 있다.
  # [1] 단순한 방법
  df.agg = aggregate(all ~ species, df, max)
  df.max = merge(df.agg, df)
  
  # [2] 또는 dplyr을 사용한 방법
  library(dplyr)
  x = df %>% group_by(species) %>% filter(all==max(all)) %>% arrange(species)
  x # 필요하면 파일로 저장한다.
  
  # quick descriptive information  https://ademos.people.uic.edu/Chapter11.html
  with(df,summary(species))
  with(df,summary(IPR039697))
  with(df,summary(all))
  df[df$species=="Lactobacillus brevis",]
  df[df$species=="Lactobacillus brevis",]$all
  summary(df[df$species=="Lactobacillus brevis",])
  
  # dplyr을 사용하면 데이터 프레임에 대한 탐색적 분석을 쉽게 할 수 있다.
  # 데이터 프레임을 tbl_df(‘tibble’) 형태로 변환하면 좀 더 쉽게 다룰 수 있다.
  df_df = tbl_df(df)
  select(df_df, accession, species, all) # 선택된 column을 출력
  filter(df_df, all > 10) # filter()는 조건에 맞는 row를 출력
  df_df %>% group_by(species) %>% summarise(number_of_strains=n()) %>% arrange(species)
  df_df %>% group_by(species) %>% summarise(max(all)) %>% arrange(species)

마지막으로 각 species에 대하여 family에 정의된 InterPro entry 검출 결과를 집계하여 error bar(표준편차)를 포함한 막대그래프를 그려 보자. ggplot2 라이브러리를 로드해야 한다.

  library("ggplot2")
  rm(list=ls())
  # https://stackoverflow.com/questions/32984974/add-error-bars-to-a-barplot 
  # reading from pre-existing file
  data = read.table("results.txt",sep="\t")
  nrow(data)
  species_list = levels(data$species)
  ipr_entries = colnames(data)
  ipr_entries = ipr_entries[-1]
  
  myData = aggregate(data[,2:7],by=list(data$species),FUN=function(x) c(mean=mean(x), sd=sd(x)))
  
  df = data.frame()
  j=1
  for (i in species_list) {
  x = matrix(data=as.vector(unlist(myData[j,2:7])),ncol=2,byrow=T)
  df.tmp = data.frame(species=rep(i,length(ipr_entries)),ipr=ipr_entries,mean=x[,1],sd=x[,2])
  df = rbind(df,df.tmp)
  
  j = j + 1
  }
  
  write.table(df,"species_mean_sd.txt",sep="\t",quote=F)
  
  # https://www.r-graph-gallery.com/4-barplot-with-error-bar.html  
  # https://www.bioinformatics.babraham.ac.uk/training/ggplot_course/Introduction%20to%20ggplot.pdf
  
  # geom_bar(color="black")을 설정하면 검정색 선으로 테두리가 생긴다.
  plot1 = ggplot(df,aes(species,mean,fill=ipr)) + geom_bar(position="dodge",stat="identity",size=1)
  # axis label을 세로로 세우기.
  plot1 = plot1 + theme(axis.text.x=element_text(angle=90,hjust=1)) + geom_hline(yintercept=28, linetype="dashed", size=0.3)
  plot1
  plot2 = plot1 + geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd),width=0.5,position=position_dodge(width=0.9),size=0.2)
  plot2
  plot3 = plot2 + xlab("Species") + ggtitle("Distribution of ADH/ALDH across probiotic bacterial species") + theme(plot.title=element_text(hjust=0.5))
    
  pdf("final_plot.pdf",width=11,height=7)
  plot3
  dev.off()

===== eggNOG-mapper를 사용한 orthology assignment 기반의 기능 주석화(functional annotation) =====

eggNOG ([[https://academic.oup.com/nar/article/51/D1/D389/6833261|v6.0]]는 상동성 관계와 유전자의 진화 역사 및 기능 주석 정보의 데이터베이스이다. [[http://eggnog6.embl.de/|웹사이트]]에서는 DB 자체에 대한 검색 및 query 서열을 입력하여 주석화를 실시할 수 있다. Ortholgy assignment에 의한 기능 주석화를 수행하는 도구는 eggNOG-mapper [[https://academic.oup.com/mbe/article/38/12/5825/6379734?login=false|v2.0]]이다. [[http://eggnog-mapper.embl.de/|eggNOG-mapper 웹사이트]]에서 단일 서열 혹은 파일 업로드를 통한 배치 주석 작업을 실시할 수 있으며, 좀 더 빠른 실행을 위해 로컬 서버에 eggNOG-mapper를 설치하여 사용하는 것도 가능하다. 매우 빠른 BLAST 호환 검색 프로그램인 DIAMOND가 있어야 eggNOG-mapper를 실행할 수 있다. 최신 버전인 eggNOG-mapper v2의 상세한 설명은 [[https://github.com/eggnogdb/eggnog-mapper/wiki/|위키 사이트]]를 참조하라. Query protein의 수가 >100M라면 FASTA 서열 파일(single-line FASTA)을 잘게 나누어서 처리하는 방법을 권장한다([[https://github.com/eggnogdb/eggnog-mapper/wiki/eggNOG-mapper-v2.1.5-to-v2.1.11#user-content-Setting_up_large_annotation_jobs|'Setting up large annotation jobs']]).

eggNOG-mapper의 기본적인 사용법은 다음과 같다. CPU 수를 지정하지 않으면 2개를 사용한다.

  $ python /data/apps/eggnog-mapper/emapper.py -i INFILE.faa --output INFILE_eggNOG -m diamond --cpu 16
  …
  Done
     INFILE_eggNOG.emapper.seed_orthologs
     INFILE_eggNOG.emapper.annotations
  Total time: 1551.87 secs

Query 서열에 대한 주석 정보는 22개의 컬럼으로 이루어진 .annotations 파일에 기록된다. 각 컬럼에 대한 설명은 다음과 같다.

  -  query_name
  - seed eggNOG ortholog
  - seed ortholog evalue
  - seed ortholog score
  - Predicted taxonomic group
  - Predicted protein name
  - Gene Ontology terms
  - EC numberKEGG_ko
  - KEGG_Pathway
  - KEGG_Module
  - KEGG_ReactionKEGG_rclass
  - BRITE
  - KEGG_TC
  - CAZy 
  - BiGG Reaction
  - tax_scope: eggNOG taxonomic level used for annotation
  - eggNOG OGs
  - bestOG (deprecated, use smallest from eggnog OGs)
  - COG Functional Category
  - eggNOG free text description

eggNOG mapper에서 출력한 annotation 파일에서 KEGG ko 번호(KEGG orthology)를 추출하면 [[https://www.genome.jp/kegg/mapper/reconstruct.html|KEGG Mapper – Reconstruct]]에 업로드할 수 있는 gene list file을 만들 수 있다. 하나의 단백질에 부여된 복수의 ko 번호를 분리하여 여러 라인으로 출력하는 것이 핵심이다. 다음의 awk one-liner를 사용하면 된다.

  $ awk -F"\t" -vOFS="\t" '$9~/^ko:/{id = "\n"$1"\t"; sub("ko:","",$9); gsub(",ko:",id,$9); print $1, $9}' eggNOG.emapper.annotations > gene_list_for_kegg_mapper.txt
===== Genomic island 예측 =====
[[https://www.pathogenomics.sfu.ca/islandviewer/upload/|IslandViewer 웹사이트]]에 annotation이 끝난 유전체의 GenBank 파일을 업로드하여 예측한다. CDS primary tag 내부에 translation 정보가 필요하며, prokka 및 dfast_core는 이를 모두 충족하는 GenBank 파일을 제공한다. Draft genome sequence을 이용한 예측은 원래 권장되지 않으나 사용자들의 요청에 의하여 그 기능이 추가되었다. 이 경우에는 사용자가 제시한 reference genome에 맞추어 먼저 정렬을 실시하여 pseudochromosome을 만들어서 분석을 개시한다.
===== 항생제 내성 유전자(AMR determinant)의 예측 =====
===== ResFinder 사용하기 =====
항생제 내성 유전자 DB 및 분석 툴로서 대표적인 것은 맥마스터 대학교의 [[https://card.mcmaster.ca/|The Comprehensive Antibiotic Resistance DB(CARD)]]와 덴마크 공대 Center for Genomic Epidemiology(CGE)의 [[https://cge.cbs.dtu.dk/services/ResFinder/|ResFinder]]가 있다. 여기에서는 웹 서비스를 사용하지 않고 Resfinder DB와 검색 도구(ResFinder)를 사용하여 획득된 유전자에 의하여 발생하는 내성을 예측하는 방법을 알아본다. 덴마크 공대의 PointFinder는 약제 내성 결핵균(Mycobacterium tuberculosis)의 경우에서처럼 chromosomal gene의 자발적 돌연변이에 의한 내성을 찾는 도구이다. resfinder.py 실행시 -db <DATABASES> 옵션을 별도로 주지 않으면 모든 항생제 내성 DB를 사용하며, 쉼표를 이용하여 복수의 내성 DB를 지정할 수 있다. Threshold와 min_coverage의 기본값은 각각 0.9와 0.6이다. 

  $ mkdir out_all out
  $ resfinder.py -i test.fsa -p /nas/DB/resfinder_db -o out_all
  $ resfinder.py -i test.fsa -p /nas/DB/resfinder_db -d aminoglycoside,macrolide -o out
  
CGE 웹사이트에서 서비스하는 것은 resfinder.pl 스크립트이다. 이것은 실행 옵션이 다르고(예: 출력물 디렉토리를 사전에 생성할 필요가 없음) blastall을 사용하며, 모든 내성 DB에 대한 검색을 할 수가 없다. resfinder.pl에 직접 옵션을 주면 스크립트에 내장된 blastall 옵션(-p blastn -a 5 -F F)에 우선하여 적용된다. 그러나 -a <num_threads> 옵션 이외의 것을 바꾸는 것은 권장되지 않는다. 
====ABRicate를 이용한 항생제 내성 병원성 인자(virulence factor) 유전자 예측 ====
[[https://github.com/tseemann/abricate|ABRicate]]는 contig 혹은 유전체 서열을 대상으로 NCBI, CARD, ARG-ANNOT, Resfinder, MEGARES, EcOH, PlasmidFinder, Ecoli_VF 및 VFDB 등 다양한 DB에 대한 검색을 실시하여 항생제 내성과 병원성 인자를 예측하는 프로그램이다. Raw sequencing read(FASTQ)는 처리하지 못한다. 뒤에서 다룰 TORMES pipeline에 포함되어 있지만, 이 프로그램은 raw sequencing read을 입력물로 받아서 트리밍 등 전처리와 조립부터 출발하는 용도로 쓰인다. 따라서 이미 다른 방식으로 조립이 완료된 contig 서열 파일이 있다면 ABRicate를 사용하는 것이 편리하다. 

  $ abricate --list
  DATABASE	SEQUENCES	DBTYPE	DATE
  ecoh	597	nucl	2018-Oct-20
  card	2237	nucl	2018-Oct-20
  ncbi	4579	nucl	2018-Oct-20
  vfdb	2597	nucl	2018-Oct-20
  plasmidfinder	263	nucl	2018-Oct-20
  resfinder	3021	nucl	2018-Oct-20
  ecoli_vf	2701	nucl	2018-Oct-20
  argannot	1749	nucl	2018-Oct-20
  $ abricate -db card contigs.fa 
  Using nucl database card:  2237 sequences -  2018-Oct-20
  #FILE	SEQUENCE	START	END	GENE	COVERAGE	COVERAGE_MAP	GAPS	%COVERAGE	%IDENTITY	DATABASE	ACCESSION	PRODUCT
  Processing: contigs.fa
  Found 5 genes in contigs.fa
  contigs.fa	Enterococcus	514111	515397	efmA	1-1287/1287	===============	0/0	100.00	100.00	card	AB467372.1:285-1572	efmA is an MFS transporter permease in  E. faecium.
  contigs.fa	Enterococcus	1922156	1923199	efrB	45-1086/1086	========/======	6/6	95.76	75.05	card	HG970103.1:1-1087	efrB is a part of the EfrAB efflux pump and both efrA and efrB are necessary to confer multidrug resistance.
  contigs.fa	Enterococcus	1923857	1925292	efrA	1-1434/1457	========/======	6/8	98.22	75.12	card	HG970100.1:1-1458	efrA is a part of the EfrAB efflux pump and both efrA and efrB are necessary to confer drug resistance.
  contigs.fa	Enterococcus	2040199	2040747	AAC(6')-Ii	1-549/549	===============	0/0	100.00	99.82	card	L12710:1-550	AAC(6')-Ii is a chromosomal-encoded aminoglycoside acetyltransferase in Enterococcus spp.
  contigs.fa	Enterococcus	2426678	2428156	msrC	1-1479/1479	===============	0/0	100.00	94.46	card	AF313494:1-1480	msrC is a chromosomal-encoded ABC-efflux pump expressed in Enterococcus faecium that confers resistance to erythromycin and other macrolide and streptogramin B antibiotics.

여러 샘플에 대하여 ABRicate를 실행하였다면 abricate --summary 옵션을 이용하여 gene presence/absence matrix를 만들 수 있다. Contig 서열 파일에 대하여 각각 검색을 실시한 결과를 따로 갖고 있어도 되고, 한 파일에 여러 샘플에 대한 결과가 수록되어 있어도 무방하다.
===== 이차대사물 생합성 유전자(biosynthetic gene cluster, BGC) 예측 =====

[[https://antismash.secondarymetabolites.org/#!/start|antiSMASH bacterial version]]에 GenBank 파일을 업로드하여 분석을 실시한다. 웹서버에서는 동시에 돌릴 수 있는 작업의 수에 한도가 있으므로, local server에 설치하여 사용하는 것도 좋다. 공개된 미생물 유전체에 대하여 미리 예측된 BGC 정보를 검색하려면 [[https://antismash-db.secondarymetabolites.org/#!/start|antiSMASH database]]를 이용하라.