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.
====== Roary를 이용한 미생물 유전체의 비교 분석 ======
이 문서는 Roary를 처음 접하는 사람을 위한 입문서가 아니므로 Roary및 이와 더불어 필요한 다른 프로그램의 설치 방법 및 기본 사용법에 대해서는 따로 설명하지 않는다. Pan genome pipeline인 Roary의 기능을 아주 간단하게 정의한다면 '(Roary) takes annotated assemblies in GFF3 format and calculates the pan genome'이 된다.
  * (논문) Roary: Rapid large-scale prokaryote pan genome analysis", Bioinformatics, 2015 Nov 15;31(22):3691-3693 [[http://doi.org/10.1093/bioinformatics/btv421|doi]]
  * [[https://sanger-pathogens.github.io/Roary/|Roary 공식 웹사이트]]
  * [[https://github.com/sanger-pathogens/Roary/blob/master/README.md|Roary GitHub]] 더 이상 지원을 할 수 없다는 공지가 맨 위에 보일 것이다. 다시 말하여 더 이상 손을 볼 곳이 없는 소프트웨어라는 뜻도 된다. 
  * [[https://github.com/microgenomics/tutorials/blob/master/pangenome.md|Step-by-step tutorial]]

Core gene은 roary에 사용한 균주 전부에 존재하는 유전자가 아니라 99% 이상의 균주에 나타나는 유전자로 정의된다. 따라서 100개 이상의 균주를 분석에 사용하였다면 특정 core gene을 갖지 않는 균주가 나올 수도 있는 것이다. 다소 어색하게 느껴지겠지만 실용적으로는 별 문제가 없다.
===== 샘플 파일 준비 =====
//Paenibacillus polymyxa//에 속하는 유전체 10개(genome analysis로는 //Paenibacillus// 속의 다른 종으로 보아야 함; [[https://www.frontiersin.org/articles/10.3389/fmicb.2019.00467/full|참조 논문]])에 대한 RefSeq accession을 파일로 만든 뒤 ncbi-genome-download를 이용하여 GenBank flat file을 다운로드한 뒤 BioPerl에 포함된 bp_genbank2gff3.pl 스크립트를 사용하여 [[https://asia.ensembl.org/info/website/upload/gff3.html|GFF3]] 파일로 전환한다.

  $ cat accessions.txt
  GCF_000146875.3
  GCF_000463565.1
  GCF_000520775.1
  GCF_000686825.1
  GCF_000714835.1
  GCF_001707685.1
  GCF_001709075.1
  GCF_001709135.1
  GCF_001956155.1
  GCF_900116035.1
  $ ncbi-genome-download -F genbank -A accessions.txt -o genbank --flat-output bacteria
  $ cd genbank
  $ gzip -d *gz
  # GCF_000463565.1_S6_genomic.gbff.gz -> GCF_000463565.1.gbk로 파일명 변경
  $ for f in $( ls ); do mv $f $( cut -d_ -f1,2 <<<$f ).gbk; done
  $ bp_genbank2gff3.pl *gbk # 또는 'bp_genbank3gff3.pl --dir .'
  $ rename 's/.gbk.gff/.gff/' *gff
맨 마지막의 rename 명령어는 Ubuntu에서 apt로 별도 설치해야 하는 패키지이다. CentOS에 기본적으로 설치된 rename과는 사용법이 약간 다르다.
  $ rename .htm .html *.htm      # CentOS
  $ rename 's/\.htm/.html/' *htm # Ubuntu

==== 다른 방법으로 마련한 통해 마련한 GFF3을 쓰면 안되나? ====
GFF3 파일은 이를 만들어내는 프로그램에 따라서 형식이 조금씩 다르다. Roary가 사용하는 GFF3 파일은 feature에 대한 정보를 지나쳐서 뒷부분의 '##FASTA'라는 선언을 시작으로 전체 염기서열과 각 단백질의 서열이 포함되는 형태이다. Prokka가 만든 GFF 파일은 큰 문제 없이 사용 가능하지만 NCBI에서 직접 내려받은 GFF 파일은 곤란하다. 또한 PGAP(local)이 만든 GenBank 파일을 GFF3로 전환하여 그대로 쓰면 일부 필드가 비어 있어서 에러가 발생할 것이다. 다음은 로컬 머신에서 PGAP을 실행하여 얻은 GenBank 파일의 앞부분 사례이다. Accession과 version 필드가 빈 상태이기 때문에 이를 그대로 GFF3로 전환한 뒤 Roary를 실행하면 제대로 결과가 나오지 않는다.
  LOCUS       chromosome           3744397 bp    DNA     circular BCT 15-MAR-2022
  DEFINITION  Ralstonia pseudosolanacearum strain SL1931 chromosome, complete
              genome.
  ACCESSION
  VERSION
  KEYWORDS    .
  SOURCE      Ralstonia pseudosolanacearum
    ORGANISM  Ralstonia pseudosolanacearum
              Bacteria; Proteobacteria; Betaproteobacteria; Burkholderiales;
              Burkholderiaceae; Ralstonia.
  REFERENCE   1  (bases 1 to 3744397)
따라서 이 필드를 적당히 채워넣어야 한다. 적당한 스크립트를 만들어서 이를 채워야 하며, 만약 multi-record GenBank 파일이라면 일이 더욱 많아진다.
===== roary 실행 =====
gff 서브디렉토리를 만들어서 GFF 파일을 전부 복사한 뒤에 roary를 실행한다. Roary 웹사이트에서는 현 디렉토리에 모든 GFF 파일이 있는 것으로 가정하고 실행하는 사례를 보여주었다('roary *gff'). GFF 파일이 만약 수백 개라면, 결과 파일과 GFF 파일이 한데 뒤섞여서 매우 보기에 좋지 않을 것이다. GFF 파일을 하나의 디렉토리에 전부 모은 뒤 다음과 같이 실행하면 훨씬 깔끔하다. 실행이 잘못되었을 경우 기존의 결과물을 지우고 재실행을 할 때에도 GFF 파일이 별도의 디렉토리에 있는 것이 훨씬 편하다. Roary는 별도의 출력 디렉토리 없이 현 디렉토리에 모든 결과 파일을 써 버리기 때문이다. *.gff 대신 분석에 사용할 GFF 파일 전부를 1.GFF 2.GFF 3.GFF...와 같이 공백을 사이에 두고 나열해도 된다. 
  $ roary –p 8 /path/to/*.gff  # thread 8개 사용; core gene alignment 생략
==== GFF 파일 수정에 따른 주의 사항 ====
십중팔구는 roary 실행과 더불어 다음의 메시지를 보게 될 것이다.
  2022/06/18 22:04:00 Input file contains duplicate gene IDs, attempting to fix by adding a unique suffix, new GFF in the fixed_input_files directory: /home/hyjeong/roary_test/gff/GCF_001956155.1.gff
Roary 실행을 마친 뒤 **query_pan_genome** 스크립트를 이용하여 GFF 파일을 다룰 때에는 수정된 GFF3 파일을 이용해야 된다. 이 수정 과정을 통하여 각 유전자의 고유 ID는 'prefix_숫자'의 형태로 바뀐다. 각 GFF3 파일에 대하여 고유하게 부여되는 이 접두사는 roary를 재실행하면 바뀐다. 실제로는 다음 사례와 같이 상당히 길다.

  f62c63acf1630b40d04ccf54718a60ec
  # 수정된 GFF 파일 내에서의 모습
  ...ID=f62c63acf1630b40d04ccf54718a60ec_19222;Parent=PPE_RS23970.t01
==== Core gene의 multiple sequence alignment 없이 빠르게 실행 ====
Core gene alignment를 이용하여 계통수를 그리지 않을 예정이라면 이 방법을 쓴다.
  $ roary -p 8 gff/*gff
  $ ls
  accessory.header.embl             fixed_input_files/
  accessory.tab                     genbank/
  accessory_binary_genes.fa         gene_presence_absence.Rtab
  accessory_binary_genes.fa.newick  gene_presence_absence.csv
  accessory_graph.dot               gff/
  blast_identity_frequency.Rtab     number_of_conserved_genes.Rtab
  clustered_proteins                number_of_genes_in_pan_genome.Rtab
  core_accessory.header.embl        number_of_new_genes.Rtab
  core_accessory.tab                number_of_unique_genes.Rtab
  core_accessory_graph.dot          summary_statistics.txt 
Core와 accessory gene의 수치를 가장 먼저 확인하도록 하자.
  $ cat summary_statistics.txt
  Core genes      (99% <= strains <= 100%)        3228
  Soft core genes (95% <= strains < 99%)  0
  Shell genes     (15% <= strains < 95%)  3290
  Cloud genes     (0% <= strains < 15%)   3814
  Total genes     (0% <= strains <= 100%) 10332
Core gene의 수가 너무 적거나 0이라면, 한 종으로 볼 수 없는 균주의 유전체를 한꺼번에 roary에 투입했기 때문일 것이다. 물론 필요하다면 일부러 동일 종에 속하지 않는 균주로부터 적은 수의 core gene을 추출하여 활용할 수 있다.
==== Multiple sequence alignment를 포함하는 표준적인 실행 ====
'-e %%--%%mafft'가 아니라 '-e'라고만 입력하면 codon-aware alignment를 하느라 엄청난 시간이 소요된다.
  $ roary -e --mafft -p 8 gff/*gff
  $ ls
  accessory.header.embl             fixed_input_files/
  accessory.tab                     genbank/
  accessory_binary_genes.fa         gene_presence_absence.Rtab
  accessory_binary_genes.fa.newick  gene_presence_absence.csv
  accessory_graph.dot               gff/
  blast_identity_frequency.Rtab     number_of_conserved_genes.Rtab
  clustered_proteins                number_of_genes_in_pan_genome.Rtab
  core_accessory.header.embl        number_of_new_genes.Rtab
  core_accessory.tab                number_of_unique_genes.Rtab
  core_accessory_graph.dot          pan_genome_reference.fa
  core_alignment_header.embl        summary_statistics.txt
  core_gene_alignment.aln

Core gene alignment와 더불어 어떤 파일이 더 생겼는지를 관찰해 보라. '-z' 옵션과 함께 실행하면 중간에 생성된 파일이 삭제되지 않고 고스란히 남는다. 특히 모든 core gene에 해당하는 <gene>.fa와 <gene>.fa.aln이 pan_genome_sequences라는 디렉토리 아래에 남게 되므로 필요한 때에 적절히 활용할 수 있다.

==== 인수(GFF 파일)를 제공하는 다양한 방법 ====
GFF 파일이 모두 현재 디렉토리에 있다고 가정하자.
  $ roary -e --mafft -p 8 *gff
  # 일부 GFF 파일에 대해서만 roary 실행
  $ roary -e --mafft -p 8 GCF_000146875.3.gff  GCF_001709075.1.gff GCF_000463565.1.gff  GCF_001709135.1.gff 
  # Roary를 실행할 일부 GFF 파일(4개)의 목록이 list.txt 파일에 있음
  $ cat list.txt 
  GCF_000146875.3.gff
  GCF_001709075.1.gff
  GCF_000463565.1.gff
  GCF_001709135.1.gff
  $ roary -e --mafft -p 8 $( cat list.txt | paste -sd' ' )
  $ roary -e --mafft -p 8 $( paste -sd' ' list.txt )

==== 결과물 이해하기 ====
[[https://sanger-pathogens.github.io/Roary/|Roary 공식 웹사이트]]의 <output files> 항목에 꽤 상세한 설명이 나온다. 여기에서는 공식 문서에서도 해소되지 않는 궁금증을 풀어 나가도록 한다. 

**clustered_proteins**는 가장 이해하기 쉬운 파일이다. **summary_statistics.txt**에 나타난 total genes는 결국 gene cluster(singleton 포함)이며, 각 클러스터를 구성하는 균주별 유전자의 목록이 clustered_proteins 파일에 라인 단위로 보여진다. 클러스터의 대표 서열은 **pan_genome_reference.fa** 파일에 수록된다. 10개의 genome 중에서 어떤 기준으로 다음의 것이 선정되었는지는 나도 잘 모른다. 아마 Roary 논문에 이에 대한 설명이 있을 것이다. 다음의 사례에서 서열 ID는 실제로 선정된 유전자이고, description 항목의 'dna'는 cluster ID에 해당한다.
  >f62c63acf1630b40d04ccf54718a60ec_3 dnaA
  GTGGACAGCCATACCTCTGAACTATGGCAGCAAATTCTATCCATTATACAAACCAAGCTG
  AGTAAGCCGAGTTACGACACTTGGTTTAAGGCTACCAAGGCAGCGAAACTAAATGACCAC
  ..

Cluster ID는 각 클러스터를 구성하는 유전자에서 가장 빈번하게 나타나는 gene ID를 따른다고 한다. 그러나 나는 이것에 100% 동의하지는 않는다. 예를 들어 이번 분석에서 EL23_RS11730라는 cluster ID가 나왔다. 이는 GCF_000714835.1.gff 하나에서만 존재하는 RefSeq의 locus tag이기 때문이다. 사실 따지고 보면 locust tag은 GFF에서 필수적으로 사용하는 식별자도 아니다. 만약 적당한 것이 없으면 'group_숫자'의 형식을 갖는다.
===== 표준적인 Roary 결과의 사용 =====
Input genome의 수를 늘려가면서 core/accessory gene의 숫자가 어떻게 변화하는지 보여주는 것은 해당 종이 open/closed pan genome 구조를 갖는지를 대략 보여주는 것과 장식적인 효과 말고는 큰 의미가 없다고 본다. Roary 결과의 활용 중에서 가장 중요한 것은 core gene alignment를 이용한 계통수 작성임에는 이의가 없을 것이다.

Roary의 GitHub 사이트에서 제공하는 [[https://github.com/sanger-pathogens/Roary/tree/master/contrib/|contrib]] 스크립트(예: roary2svg.pl)는 우분투의 roary package에는 없으니 필요하면 별도로 설치하여 사용하라.

===== query_pan_genome 스크립트 활용 아이디어 및 심화 과정 =====
==== Core gene 찾기 ====

앞에서도 언급했듯이 query_pan_genome 스크립트에 인수로 제공해야 할 GFF 파일은 roary가 실행되면서 수정을 거친 것이어야 한다. 본 위키 페이지에서 다룬 사례에서는 10개의 GFF 파일 전부 수정되어 fixed_input_files 디렉토리에 저장되었으므로, core gene을 찾으려면 다음과 같이 실행해야 한다.

  $ query_pan_genome -a intersection fixed_input_files/*.gff

그러나 수십 개의 GFF 파일을 사용했다면, 그중에서 일부만 수정되어 fixed_input_files 디렉토리에 지정될 것이다. 따라서 정상적인 GFF 파일과 수정을 거친 GFF 파일을 한 디렉토리(예: dir_all/)에 복사해야 dir_all/*gff 인수 형태로 한번에 query_pan_genome 스크립트를 실행할 수 있다.

  $ mkdir dir_all
  $ cp gff/gff* fixed_input_files/*gff dir_all
  $ query_pan_genome -a intersection dir_all/*gff

만약 roary에 사용했던 10개의 genome 중에서 임의로 지정한 일부에 대해서 core gene을 찾고자 한다면? dir 아래에 원하는 GFF 파일만 복사한 상태에서 'query_pan_genome -a <command> dir/*gff'을 실행하면 된다. 또는 다음과 같이 공백을 구분자로 하여 작업에 사용할 GFF 파일을 나열해도 된다. 목록 파일을 만들어서 사용해도 좋다.
  $ query_pan_genome -a intersection 1.gff 2.gff 3.gff 4.gff 5.gff
 

==== Comma-separated list 만들기 ====
query_pan_genome에서 쉼표로 분리된 유전자 또는 GFF 파일 목록을 인수로 제공해야 할 때가 있다.

  # Extract the sequence of each gene listed and create multi-FASTA files
  # 실제로 출력되는 것은 아미노산 서열이다.
  $ query_pan_genome -a gene_multifasta -n gryA,mecA,abc *.gff
  # Gene differences between sets of isolates
  $ query_pan_genome -a difference --input_set_one 1.gff,2.gff --input_set_two 3.gff,4.gff,5.gff
  
기껏해야 너댓 개의 단어를 나열해야 한다면 손으로 타이핑하면 그만이다. 그러나 열 개를 넘어가면 여간 고통스럽지 않을 것이다. 오타가 발생할 것도 너무나 병백하다. 따라서 나열할 단어를 라인 단위로 텍스트 파일에 저장한 다음, paste 명령어를 이용하면 된다. 텍스트 파일은 물론 각자 목적에 맞추어 만들어야 한다. 실습 목적으로 01.gff, 02.gff...10.gff를 라인 단위로 수록한 list.txt 파일을 만든 뒤, 이를 '01.gff,02.gff,03.gff...,10.gff' 형태의 문자열로 만들어 사용자 환경변수 MY_LIST에 저장해 보자
  $ for i in $( echo {01..10} )
  > do
  > echo $i.gff >> list.txt
  > done
  $ cat list.txt
  01.gff
  02.gff
  03.gff
  04.gff
  05.gff
  06.gff
  07.gff
  08.gff
  09.gff
  10.gff
  $ MY_LIST=$( paste -sd, list.txt )
  $ echo $MY_LIST
  01.gff,02.gff,03.gff,04.gff,05.gff,06.gff,07.gff,08.gff,09.gff,10.gff
이는 내가 매우 즐겨 사용하는 기법이다. 첫 줄에서 01, 02... 10을 생성하는 명령어인 'echo {01..01}'은 'seq -w 2 1 10'으로 대체할 수 있다.

set_one과 set_two가 반드시 상호 배타적일 필요는 없을 것이다. 다시 말하여 set_one + set_two = '전체'가 되도록 설정을 해도 되고, 전체의 부분집합이 되게 설정을 할 수도 있다. set_one이 정의되어 있을 경우 set_two(전체 - set_one)을 만드는 방법을 알아보자. 여기에는 [[https://tldp.org/LDP/abs/html/process-sub.html|process substitution]]과 **comm** 명령어 활용 등 상당한 수준의 기법이 쓰였다.
  $ cat list_all
  GCF_000146875.3.gff
  GCF_000463565.1.gff
  GCF_000520775.1.gff
  GCF_000686825.1.gff
  GCF_000714835.1.gff
  GCF_001707685.1.gff
  GCF_001709075.1.gff
  GCF_001709135.1.gff
  GCF_001956155.1.gff
  GCF_900116035.1.gff
  $ cat list_set_one
  GCF_000146875.3.gff
  GCF_000714835.1.gff
  GCF_001956155.1.gff
  $ comm -13 <( sort list_set_one ) <( sort list_all ) > list_set_two  
  $ paste -sd, list_set_one
  GCF_000146875.3.gff,GCF_000714835.1.gff,GCF_001956155.1.gff
  $ paste -sd, list_set_two
  GCF_000463565.1.gff,GCF_000520775.1.gff,GCF_000686825.1.gff,GCF_001707685.1.gff,GCF_001709075.1.gff,GCF_001709135.1.gff,GCF_900116035.1.gff

마지막 두 개의 paste 명령어는 그 결과를 화면에 보여주기 위한 것이다. 실제로는 다음과 같이 출력을 환경 변수에 넣어서 사용하면 된다('VAR=$( ... )').
  $ set_one=$( paste -sd, list_set_one )
  $ set_two=$( paste -sd, list_set_two )
  
혹시 comma-separated GFF file list를 query_pan_genome 명령어의 *gff 인수를 대체해서 사용할 수 있지 않을까하는 생각이 들 수도 있을 것이다. 그러나 불행히도 그렇게는 작동하지 않는다! 단지 'query_pan_genome -a difference' 명령에서 '%%--%%input_set_one'과 '%%--%%input_set_two'에 바로 이어서 나오는 GFF 파일 목록에만 사용할 수 있다. 또한 GFF file이 현 디렉토리에 있지 않다면 앞에 path가 붙도록 만들어야 한다. 

==== gene_presence_absence.csv 파일을 R에서 뜯어 보자 ====
(작성 예정)
==== Strain/group-specific gene을 찾는 가장 현명한 방법은 무엇일까? ====
'query_pan_genome -a difference' 스크립트를 쓰는 것이 가장 공식적인 방법일 것이다. 그러나 아무리 단순한 set_one이 주어진다 하더라도 이에 맞추어 set_two를 만들어야 하는 것이 다소 성가실 수도 있다. 바로 위에서 설명한 방법을 이용하여 set_one과 set_two 환경 변수에 GFF 파일의 목록을 저장하였다고 가정하자. Path를 붙이지 않은 파일명을 사용하려면 fixed_input_files에서 query_pan_genome 실행 시에 '-g' 옵션으로 clustered_proteins 파일의 위치를 정확히 지정해야 한다.

  $ query_pan_genome -g ../clustered_proteins -a difference --input_set_one $set_one --input_set_two $set_two
  $ ls set_difference*
  set_difference_accessory_graph.dot
  set_difference_common_set
  set_difference_common_set_reannotated
  set_difference_common_set_statistics.csv
  set_difference_core_accessory_graph.dot
  set_difference_unique_set_one
  set_difference_unique_set_one_reannotated
  set_difference_unique_set_one_statistics.csv
  set_difference_unique_set_two
  set_difference_unique_set_two_reannotated
  set_difference_unique_set_two_statistics.csv
set_one은 3개, set_two는 7개의 GFF 파일로 구성된 상태이다. common_set두 그룹에 공통적으로 존재하는 유전자 정보를 담고 있다. 다시 말하자면, 각 그룹을 구성하는 GFF 파일 10개 전부를 통틀어서 존재할 필요는 없다. 예를 들어서 set_difference_common_set의 가장 마지막 줄을 살펴보자.
  $ tail -n 1 set_difference_common_set
  EL23_RS26840: 51172ff78909e9e65f8fb727afd887c6_691      d54d6e1622af26ac05c8a4c1103f23b3_6984
단 두개의 유전자가 보일 뿐이다. 앞의 것은 set_one에 속하는 GCF_000714835.1.gff의 것이고, 뒤의 것은 set_two에 속한다. 이 파일의 앞부분에 위치한 클러스터의 유전자 수는 10개이다. set_difference_unique 계열의 파일도 마찬가지이다. 그러므로 set을 나누었다 하더라도 각 set의 모든 GFF 파일에 존재하는 것을 찾으려면 별도의 작업을 해야 된다.

_reannotated라는 접미어가 붙은 파일과 그렇지 않은 것의 차이가 무엇인지 roary 공식 문서에서는 정확하게 설명하지 않는다. set_difference_common_set에서는 group_XXXX로 표시된 클러스터가 98개인 반면, set_difference_common_set_reannotated에서는 18개로 줄어들었다. 이들은 group_XXXX라는 클러스터 명칭 대신 gene ID를 갖게 된 것으로 확인되었다. 하지만 set_difference_common_set_reannotated 파일에서 새로 얻은 gene ID는 유일한 식별자가 아니다. 예를 들어서 zwf라는 식별자를 찾아보자.
  $ grep zwf set_difference_common_set_reannotated
  zwf: f62c63acf1630b40d04ccf54718a60ec_7917...
  zwf: f62c63acf1630b40d04ccf54718a60ec_15406...
아래의 것은 clustered_proteins와 set_difference_common_set에 원래 있던 것이다. Group ID로 쓰일 수 있도록 gene ID를 새로 부여해 준 노력은 가상하지만, 이전에 존재하던 식별자와 겹친다면 별로 도움이 되지 못한다.

=== 어쩌면 Roary 1차 결과물인 gene_presence_absence.Rtab이 더 쓸모가 있을지도 모른다 ===
Roary 1차 결과물 중 다음의 4개는 지속적으로 연구할 가치가 있다. 맨 마지막 것은 단순한 FASTA file이라서 별로 설명이 필요하지 않다. 그러나 위의 3개 파일은 유전자 클러스터를 매우 다른 방식으로 표현한다.
  * gene_presence_absence.csv
  * gene_presence_absence.Rtab
  * clustered_proteins
  * pan_genome_reference.fa: '-e' 또는 '-e %%--%%mafft' 옵션과 함께 실행해야 생성됨
gene_presence_absence.Rtab은 row에 배열된 각 클러스터가 GFF 파일에 존재하는지를 0과 1의 행렬 형태로 나타내어 준다.
  Gene    GCF_000146875.3 GCF_000463565.1 GCF_000520775.1 GCF_000686825.1 GCF_000714835.1 GCF_001707685.1 GCF_001709075.1 GCF_001709135.1 GCF_001956155.1 GCF_900116035.1
  B439_RS0117690  1       1       1       1       1       1       1       1       1       1
  EL23_RS04475    1       1       1       1       1       1       1       1       1       1
  A7309_RS08725   1       1       1       1       1       1       1       1       1       1
  ..중략..
  L694_RS0123725  0       0       1       0       0       0       0       0       0       0
  L694_RS0115790  0       0       1       0       0       0       0       0       0       0
시험삼아서 set_one(GCF_000146875.3.gff, GCF_000714835.1.gff, GCF_001956155.1.gff; list_set_one 참조)에만 존재하는 유전자를 찾아 보겠다. Bash에서 배열을 생성하고, grep의 성공 여부를 1 or 0(실제의 exit code와는 반대)으로 출력하는 트릭에 주목해 보도록 한다.
  $ gff=( $( head -n 1 gene_presence_absence.Rtab ) )
  $ unset gff[0] # 첫 컬럼인 'Gene'을 제거함
  $ echo ${gff[@]} # 확인용
  GCF_000146875.3 GCF_000463565.1 GCF_000520775.1...
  # temp.txt 파일이 있다면 미리 지우도록 한다.
  $ for i in ${gff[@]}
  > do
  > ! grep -q $i list_set_one; echo $? >> temp.txt
  > done
  $ pattern=$( paste -sd" ' temp.txt ) # TAB is the default delimiter
  $ sed 's/\t/ /g' gene_presence_absence.Rtab > temp.Rtab
  $ grep "$pattern" temp.Rtab
  EL23_RS17190 1 0 0 0 1 0 0 0 1 0
gene_presence_absence.Rtab 파일은 탭을 구분자로 사용한다. 따라서 pattern 변수에도 이를 그대로 흉내내어 탭으로 구분된 1과 0의 패턴을 넣고 싶었으나 echo 문으로는 잘 되지 않았다. 이에 따라서 임시방편으로 필드 구분자를 탭에서 공백으로 바꾸어 버린 것이다. 전체 분석을 통해서 EL23_RS17190이라는 group ID에 해당하는 유전자가 주어진 3개의 GFF 파일에서는 어떤 구체적인 ID로 불리는지 찾는 것은 별로 어렵지 않을 것이다. 아미노산 서열 파일을 찾으려면 'query_pan_genome -a gene_mulifasta -n EL23_RS17190 <GFF file>'을 쓰면 된다. <GFF file> 인수 위치에는 단일 파일 또는 *.gff만이 올 수 있다.

pattern 변수에 '1 1 1...1'과 같이 GFF 파일 수만큼 1을 이어 붙인 값을 넣으면 이를 temp.Rtab에 적용하여 core gene 목록을 추출할 수 있다.

=== R을 이용하여 간단한 일을 복잡하게 수행하기 ===
gene_presenece_absence.Rtab은 본래 R에서 읽어들여서 사용할 것을 전제로 만들어진 것이다. 유전자 cluster가 row로 배열되어 있으며 컬럼은 각 균주(GFF file)에 해당한다. 유전자를 갖고 있으면 1, 없으면 0으로 표시된다. Paralog 형태로 두 벌의 유전자가 있다고 하여 2로 표시되지는 않는다.
  > df = read.table("gene_presence_absence.Rtab",sep="\t",header=T,check=F)
  > total.genomes = ncol(df) - 1 # 분석에 사용한 전체 genome 수 확인
  # core gene을 core 벡터에 저장한다.
  > core = df[apply(df[,c(2:ncol(df))],1,sum) == total.genomes,1]
  # 텍스트 파일로 저장
  > write.table(core,"core_genes.txt",quote=F,row.names=F,col.names=F)

==== clustered_proteins 파일에서 core gene 찾을 때 주의할 사항 ====
이 파일은 cluster_id: gene_1 gene_2 gene_3...의 매우 단순한 구조로 이루어져 있다. 단, cluster_id 뒤에는 탭이 오는 반면, gene 사이의 구분자는 공백이다. 10개의 GFF 파일을 사용했으므로 총 11개의 컬럼으로 이루어진 라인을 찾아서 첫 번째 컬럼을 뽑으면 그것이 core gene의 ID가 될 것이다. 탭과 공백이 섞여 있어도 awk는 이를 알아서 잘 처리한다. 실제 컬럼의 수가 어떠한지 맨 위의 3줄을 뽑아 본다. clustered_proteins 파일을 컬럼 수의 역순으로 라인이 배열되어 있다.
  $ awk '{print $1, NF}' clustered_proteins | head -n 3
  A7312_RS16310: 12
  B439_RS0114220: 11
  A7312_RS15515: 11
이상하지 않은가? 첫 번째 줄은 컬럼이 12개이다. 균주 10개를 비교했는데 11개의 A7312_RS16310으로 대표되는 유전자가 12개  나왔다는 것은 어느 하나(혹은 그 이상)의 균주에서 이 유전자가 2벌 이상 존재(paralog; roary의 클러스터링 기준은 blastp identity > 95%)한다는 뜻이다. 극단적으로 말하자면 어느 한 균주에서 3번, 그리고 7개 균주(나머지 전체가 아님!)에서 1회씩 출현했을 수도 있는 것이다. 그러나 row를 구성하는 유전자의 수를 세어서 core gene 여부를 판별하면 매우 곤란하다.

clustered_proteins의 첫 줄을 면밀히 살펴보니 1e09b63d307cd166cc772781759c09a1_12594와 1e09b63d307cd166cc772781759c09a1_12582 유전자가 동일한 prefix를 가지므로 하나의 균주에서 유래했음이 틀림없다. 확인 결과 GCF_001707685.1.gff였다. gene_presence_absence.csv 파일에서 해당 유전자의 annotation을 확인해 본니 non-ribosomal peptide synthetase였다. Paralog처럼 보이기에 적당한 유전자이다.

==== cluster ID를 이용하여 GFF 파일에서 지정된 유전자의 염기서열 추출하기 ====
query_pan_genome 스크립트는 GFF 파일(하나 또는 여러 개)로부터 주어진 cluster ID의 아미노산 서열을 FASTA file로 추출하는 기능이 있다.
  $ query_pan_genome -a gene_multifasta -n gryA,mecA,abc *.gff
그러나 아미노산 번역본이 아니라 염기서열 원본이 필요하다면 어떻게 할 것인가? 'roary -z' 옵션으로 중간 파일을 지우지 않고 남겨놓도록 해 봐야 pan_genom_sequences 디렉토리에는 갭이 포함된 .fa.aln(aligned FASTA file; 지정된 모든 유전자에 대해서 하나씩 만들어짐)이 생성된다. 여기에는 GFF 파일에 대한 정보는 없지만 roary 실행 단계에 GFF를 수정하면서 만들어진 prefix가 쓰이고 있다. 따라서 GFF와 prefix의 연결 관계를 알면 어떻게 해서든 염기서열을 뽑을 수 있다. gene_presence_absence.csv의 첫 두 줄로부터 GFF 파일과 prefix의 연결 관계를 알아내어 보자.

  $ paste -d, <(head -n 1 gene_presence_absence.csv | tr ',' '\n') <(sed -n 2p gene_presence_absence.csv | tr ',' '\n') | sed 's/"//g' | grep GCF_
  GCF_000146875.3,f62c63acf1630b40d04ccf54718a60ec_16264
  GCF_000463565.1,4c9bd31f831baeca0dfbb87d8830abc7_13550	4c9bd31f831baeca0dfbb87d8830abc7_13554
  GCF_000520775.1,32c86718b06a62051d7fa2065f093252_8846
  GCF_000686825.1,0b99fa4abde2729537f83e788213eb51_20393
  GCF_000714835.1,eb031fc20295abfb065c0b01de926f07_19675
  GCF_001707685.1,1e09b63d307cd166cc772781759c09a1_12594
  GCF_001709075.1,2eb394db95b4da5564b60273e929ddb2_17280
  GCF_001709135.1,47ccc10b2f363d0e62e4554689c165de_18880
  GCF_001956155.1,501382f1bb77eb72358f1cc21fbfba74_12153
GCF_000463565.1에는 두 개의 유전자가 있다. 이는 paralog를 의미한다.

clustered_proteins 파일에서 추출하려는 cluster의 ID를 알면, 이에 해당하는 줄로부터 추출하려는 GFF에 해당하는 prefix를 알 수 있다. 따라서 .fa.aln을 찾아가서 유전자 염기서열을 뽑아낸 뒤 갭을 없애면 된다.

혹은 좀 더 복잡하지만 기본에 충실한 방법을 탐구해 보자. clustered_proteins 파일을 통하여 32c86718b06a62051d7fa2065f093252_8846라는 개별 유전자가 목표임을 알게 되었다고 가정하자. Prefix와 GFF 파일의 연결 관계를 알게 되었으니 fixed_input_file에 가서 이에 해당하는 GFF 파일을 조작하면 된다. 일반적으로 GFF 파일로부터 특정 feature의 염기서열을 뽑아내려면 별도의 FASTA file이 필요하다. 그러나 우리에게는 sequence-embedded GFF 파일이 있지 않은가(##FASTA 다음 줄부터 염기서열 자료에 해당함)? 여기에서 염기서열을 뽑아내어 파일로 만들어 보자. 단, 유전체 염기서열과 이에 부속된 유전자 염기서열이 공존하므로 이를 분리하는 것이 중요하다. [[https://bioinf.shenwei.me/seqkit/|seqkit]] grep 명령은 표준 입력으로 공급되는 서열 자료의 ID에 대한 패턴 매치를 통하여 서열을 필터링하는 막강한 기능이 있다.

  $ f=GCF_000520775.1.gff
  $ f_header=${f%%.gff}
  $ sed '1,/^##FASTA/d' $f | seqkit grep -r -v -p .p01$ | seqret fasta::stdin fasta:${f_header}.fna
  $ sed '1,/^##FASTA/d' $f | seqkit grep -r -p .p01$ | seqret fasta::stdin fasta:${f_header}.ffn

이렇게 실행하면 GCF_000520775.1.fna(유전체 염기서열)와 GCF_000520775.1.ffn(유전자 염기서열)을 얻게 된다. 'ID=32c86718b06a62051d7fa2065f093252_8846;'라는 패턴으로 GCF_000520775.1.gff 파일을 검색하면 GFF 파일 내에서 정의한 유전자 ID(Name=L694_RS0112140)가 나오고, L694_RS0112140.<color #ed1c24>p01</color>을 GCF_000520775.1.ffn 파일에서 찾으면 소기의 목적을 달성한 것이다.
  $ grep 'ID=32c86718b06a62051d7fa2065f093252_8846;' GCF_000520775.1.gff
  NZ_ASSA01000213	GenBank	CDS	5500	7854	.	+	1	ID=32c86718b06a62051d7fa2065f093252_8846;Parent=L694_RS0112140.t01;Name=L694_RS0112140;Note=Derived by automated computational analysis using gene prediction method: Protein Homology.;codon_start=1;inference=COORDINATES: similar to AA sequence:RefSeq:WP_013312094.1;product=condensation domain-containing protein;protein_id=WP_025721859.1;transl_table=11;translation=length.785
GCF_000520775.1.gff 파일에서 염기서열 ID와 위치 정보를 뽑아낸 뒤, GCF_000520775.1.<color #ed1c24>fna</color> 파일에서 이에 해당하는 염기서열을 뽑아내는 것도 얼마든지 가능하다. 나는 이 방법이 더 나은 것 같다.

=== 종합 실습 ===
Core gene 목록을 뽑은 뒤 이에 해당하는 유전자 염기서열을 **GCF_000520775.1.gff**에서 뽑아내어 FASTA file로 만든다. 작업 디렉토리는 fixed_input_files로 한다. 가장 먼저 할 일은 뽑아낼 유전자(core gene)의 cluster ID를 추출하는 것이다.
  $ query_pan_genome -a intersection -g ../clustered_proteins *gff
  $ awk '{print $1}' pan_genome_results | sed 's/://' | sort > core_genes.txt
다음으로 GCF_000520775.1.gff가 gene_presence_absence.csv에서 몇 번째 컬럼에 해당하는지를 알아낸다.
  $ head -n 1 ../gene_presence_absence.csv | tr ',' '\n' | awk '{print NR, $1}' | grep GCF_000520775.1
  17 "GCF_000520775.1"
그러면 gene_presence_absence.csv 파일을 열어서 core_genes.txt 각 라인의 값(추출할 core gene의 ID)이 첫 번째 컬럼에 오는 라인에 대해서 17번째 라인을 추출하면 될까? 그렇게 간단하지는 않다. 왜냐하면 다음의 경우를 고려해야 하기 때문이다.
  - gene_presence_absence.csv는 콤마로 각 컬럼을 구분하고 있으며, 컬럼은 따옴표로 둘러싸여 있다. Annotation 컬럼은 내부적으로 콤마를 포함할 수도 있다. 많은 유전자의 product 이름에는 콤마가 들어있는 경우가 많다.
  - Paralog가 출현하면 하나의 컬럼에 두 개 이상의 유전자가 있을 수 있다. 
1, 즉 따옴표로 둘러싼 컬럼(특히 그 내부에 콤마를 포함하는 경우)을 콤마로 구분하는 csv 파일을 파싱하는 것은 매우 흔한 일이면서도 의외로 까다롭다(([[https://stackoverflow.com/questions/71543341/how-can-i-parse-csv-files-with-quoted-fields-containing-commas-in-awk|How can I parse CSV files with quoted fields containing commas, in awk?]]에 의하면 awk -v FPAT='([^,]+)|(\"[^\"]+\")' file.csv 명령을 써서 필드 구분자가 아니라 필드 자체를 매치하도록 만들면 된다.)). [[https://github.com/BurntSushi/xsv|xsv]], [[http://lorance.freeshell.org/csv/|AWK CSV parser]], [[https://csvkit.readthedocs.io/en/latest/index.html|csvkit]] 등 다양한 솔루션이 있는데, 나는 가장 마지막의 csvkit을 써 보기로 하였다. gene_presence_absence.csv의 1번 컬럼('Gene')이 core_genes.txt와 겹치는 라인의 17번 컬럼('GCF_000520775.1')을 추출하되 하나의 core gene에 대하여 유전자가 1개씩 정확히 있는지를 특히 눈여겨 보아야 한다.
  $ csvcut -c 1,17 ../gene_presence_absence.csv | sort -t, -k1,1 > gene_presence_absence_GCF_000520775.1.csv
  $ join -t, -1 1 -2 1 core_genes.txt gene_presence_absence_GCF_000520775.1.csv > core_gene_info_GCF_000520775.1.txt
**core_gene_info_GCF_000520775.1.txt** 파일을 콤마로 구분했을 때 두 번째 컬럼에 유전자가 정확히 1개씩 있는지, 혹은 탭을 구분자로 하여 2개 이상이 있는지를 확인해 둔다. 
== [1] pan_genome_sequences 디렉토리의 core gene alignment에서 염기서열을 뽑아내기 ==
  # 현 위치는 fixed_input_files
  $ cp core_gene_info_GCF_000520775.1.txt ../pan_genome_sequences/
  $ cd ../pan_genome_sequences/ 
  $ rm core_genes_from_GCF_000520775.1.fa # 있는 경우에만 
  $ awk -F, '{print $1, $2}' core_gene_info_GCF_000520775.1.txt | while read a b
  > do
  > seqkit grep -r -p $b $a.fa.aln | seqret fasta::stdin fasta:stdout | sed "s/>/>$a /" >> core_genes_from_GCF_000520775.1.fa
  > done
염기서열이 추출된 유전자 수는 하나가 적다. core_gene_info_GCF_000520775.1.txt에는 abc-f라는 core gene이 있지만 정작 pan_genome_sequences 디렉토리에 있는 파일은 abc_f.fa.aln이다. 파일명에는 
-'를 쓰지 않게 만드는 것 같다. 다소 불편하지만 실행 전후의 유전자 수를 점검해 보는 것이 좋을 것이다. core_genes_from_GCF_000520775.1.fa 내의 gap을 없내는 것은 EMBOSS의 degapseq 유틸리티를 이용한다.
  $ degapseq core_genes_from_GCF_000520775.1.fa core_genes_from_GCF_000520775.1.fa.nogaps
== [2] GFF 파일과 이로부터 추출한 FASTA file을 이용하는 방법 ==

  # 현 위치는 fixed_input_files
  $ perl test_20220622.pl core_gene_info_GCF_000520775.1.txt GCF_000520775.1.gff > core_genes_from_GCF_000520775.1.fa
Perl 스크립트 'test_20220622.pl'는 그렇게 효율적이지는 못하다. EMBOSS의 seqret 명령어를 추출할 core gene 수만큼 실행하기 때문이다. GCF_000520775.1.fna 파일을 스크립트 실행 개시 부분에서 Bio::Seq 개체로 한 번에 읽어들인 뒤, core gene이 위치한 서열 개체 정보를 이용하여 출력하는 것이 더 빠를 수도 있다. 별로 신통치 않은 스크립트이지만 소개하여 본다.

  #!/usr/bin/perl
  #
  
  open LIST, $ARGV[0];
  while ( <LIST> ) {
          chomp;
          @data = split /,/, $_;
          $seen{$data[1]} = $data[0];
  }
  
  ( $strain = $ARGV[1] ) =~ s/\.gff$//;
  open GFF, $ARGV[1];
  while ( <GFF> ) {
          chomp;
          @temp = split /\t/, $_;
          next unless $temp[2] eq 'CDS';
          $temp[8] =~ m/ID=([^;]+);/;
          $op = '';
          if ( exists $seen{$1} ) {
                  $seq = $strain . '.fna:' . $temp[0];
                  $op = '-sreverse' if $temp[6] eq '-';
                  $idndesc = '>' . $seen{$1} . ' ' . $1;
                  system( "seqret -sbegin $temp[3] -send $temp[4] $op $seq fasta::stdout | sed \"s/>.*/$idndesc/\" " );
          }
  }
두 번째 방법으로 추출한 FASTA 파일에는 gap이 존재하지 않으며, 오로지 위치를 이용하여 염기서열을 추출하게 되므로 core gene의 group ID가 약간 틀려도 문제를 일으키지 않는다는 장점이 있다. 

===== Roary에게 경의를! Scory =====