[Tutorial] Vibrio cholerae의 비교 유전체 분석

[Tutorial] 건강한 사람과 환자 간의 호흡기 마이크로바이옴 비교
07/25/2017

[Tutorial] Vibrio cholerae의 비교 유전체 분석

이 Tutorial은 BIOiPLUG(http://www.bioiplug.com/)의 비교 유전체(comparative genomics) 분석 이용 방법을 설명하기 위해 작성되었습니다. 여러분의 이해를 돕기 위하여 잘 알려진 병원성 미생물인 Vibrio cholerae (콜레라 원인균) 를 예제로 이용하였습니다.

Vibrio cholerae에 대하여

Vibrio cholerae 는 인간에게 콜레라라는 질병을 감염시키는 박테리아입니다. 콜레라는 급성 설사를 유발시키는 질병이며 적절한 치료를 받지 않을 경우 수 시간 내에 숙주를 사망에 이르게 만듭니다 (자세한 것은 WHO 사이트 참조). 그렇지만 V. cholerae 의 모든 균주가 전부 콜레라의 원인균인 것은 아니며 주로 콜레라 독소(cholera toxin = CTX) 을 포함하고 있는 콜레라 균에 의해 감염됩니다. 이 CTX는 2개의 유전자 (ctxA와 ctxB) 에 의해 만들어지는데, ctxAB 유전자를 보유하고 있는 균주는 “독소유발성”으로 간주됩니다. 흥미롭게도, 이 치명적인 유전자는 ctx 박테리오파지 (CTXφ) 라 불리는 바이러스의 유전자이며, 이는 수평적 유전자 이동(lateral gene transfer)의 주요 메커니즘인 잠재성 박테리오 파지에 의한 유전자이동(Temperate bacteriophage-mediated gene transfer)을 통해 박테리아에서 박테리아로 전파됩니다.

콜레라는 아시아와 아프리카에 걸쳐 매우 널리 확산되었는데, 이것을 콜레라 대유행 (Cholera pandemic) 라고 합니다. 콜레라 대유행은 여전히 진행 중이며 최근의 제7기 콜레라 대유행은 항원형 O1 biovar El Tor (줄여서 O1 El Tor) 라고 알려진 균주에 의해 발생했고, 한 단계 전의 제6기는 항원형 O1 biovar Classical (O1 Classical) 에 의해 발생했습니다. 이러한 클론변위(clonal shift)가 (O1 Classical에서 O1 El Tor로) 발생하는 이유는 잘 알려지지 않았으나 제 6기 콜레라 대유행 당시 어디에서나 분리할 수 있었던 O1 Classical이 더 이상 세계 어느 곳에서도 발견되지 않는다는 것은 매우 흥미로운 사실입니다.

 

시작하기 전에

  • 이 tutorial을 시작하기 전에, 자료의 퀄리티를 향상시키고 싶은 분을 위해 dendrogram과 계통수(phylogenic tree)를 보고 다룰 수 있는 데스크탑 컴퓨터 프로그램을 설치하는 것을 추천합니다. MEGA라는 프로그램이 권장되며, MS Windows와 Mac OS X 모두에서 사용할 수 있습니다.
  • http://www.bioiplug.com 또는 http://www.ezbiocloud.net에 계정을 생성해야 합니다.
  • [My Data]를 클릭하여 “Vibrio cholerae Tutorial Set”을 선택합니다.

이제 선택한 데이터 셋을 이용하여 BIOiPLUG의 comparative genomics 서비스의 주요 기능들에 대해 알아보도록 하겠습니다.

V. cholerae 균주에 대해

균주명 항원형/생물형 대유행
N16961 O1 El Tor 7th Pandemic
TSY216 O1 El Tor 7th Pandemic
IEC224 O1 El Tor 7th Pandemic
2010EL-1786 O1 El Tor 7th Pandemic
B33 O1 El Tor 7th Pandemic
MO10 O139 7th Pandemic
4260B O139 7th Pandemic
O395 O1 Classical 6th Pandemic
ATCC 14035 (Type strain of the species) O1 Classical 6th Pandemic
10432-62 O27
1154-74 O49
TM 11079-80 O1 El Tor
LMA3984-4 O1 El Tor
12129(1) O1 El Tor
HE39 non-O1/non-O139
HE48 non-O1/non-O139

Data set 빠르게 살펴보기

예제로 사용될 data set은 16개의 V. cholerae 유전체를 포함하고 있습니다. 이 16개의 균주들은 각각 임상과 환경에서 분리한 다양한 조합의 항원형, 생물형, 독소 생성성을 나타내는 균주들입니다.

  • [DATE SET]→[Genomes in This Data Set]로 이동해서 각 유전체에 대한 정보 (=Metadata)를 확인합니다. BIOiPLUG는 가능한 많은 metadata를 확보하고 검증하기 위해 노력하고 있습니다.

개별 유전체 살펴보기

개별 유전체를 살펴보기 위해, 사용자는 BIOiPLUG의 “Genome Explorer”를 실행해야 합니다. 여기서는 MO10 균주의 유전체를 열어보겠습니다. 아래 스크린샷의 (a)가 가리키는 [Open] 버튼을 클릭합시다. 이 균주는 항원형 O139에 해당하는 Vibrio Cholerae로, 1992년 인도의 Mandras에서 분리되었습니다. 유전체 염기서열은 완전히 결정되지 않았으며 27개의 컨티그로 구성되어 있습니다. 단일 유전체를 탐색하는 방법을 알아보고 싶다면 Genome Explorer를 위한 tutorial을 이용해 주시기 바랍니다.

유전체 기반 계통 분석 (Phylogenomic analysis)

Comparative genomics 분석을 하기 위한 첫 번째 단계는 균주들 간의 계통학적 관계를 자세히 살펴보는 것입니다. 이를 통해 어떤 균주들이 서로 더 가깝게 연관되어 있는지를 알 수 있습니다.

  • [PHYLOGENOMICS]→[OrthoANI]로 이동합니다.

OrthoANI는 2개의 유전체 염기서열 간의 분류학적 의미를 가지는 유사도를 측정한 값입니다. 여러 개의 유전체/균주들 간의 계통학적 유연관계를 추론하기 위한 가장 좋은 방법이 OrthoANI 값을 이용한 계층적 클러스터링(일반적으로 UPGMA)입니다.

*아래 그림은 MEGA program을 이용해서 편집된 자료임을 유의하여 봐 주시기 바랍니다. 지난 수십 년간 무수한 사람들을 사망에 이르게 만든 병원균이 바로 6기와 7기 유행성 균주입니다.

v_cholerae_orthoani

  • [PHYLOGENOMICS]→[Tetra-nucleotide]로 이동합니다.

Tetra-nucleotide analysis (TNA)는 실제 염기서열 정렬이 아닌 네 개의 뉴클레오타이드 염기서열의 조합을 이용하는 것이기 때문에 정확히는 계통 분석이 아닙니다. 그러나 이것은 유전체들이 어떻게 연관되어 있는지에 대해 유용한 추가 정보를 제공합니다. 아래에 tetra-nucleotide의 구성에 기반하여 만들어진 V. cholerae 균주들의 UPGMA 클러스터링 결과가 있습니다.

 

  1. 이 균주는 tetra-nucleotide 빈도의 조합이 다른 균주와 다소 다르게 결과가 나옵니다. 이러한 결과의 원인에 대해서는 아래에서 알아보도록 하겠습니다.

이 유행성 균주가 전 세계적으로 널리 퍼질 수 있었던 이유는 무엇인가?

아마도 가장 처음 생각해 보게 되는 것은 이러한 질문일 것입니다. 우리는 유행성 균주가 감염과 독소 생성을 가능하게 하는 유전자를 포함하고 있다고 가정할 수 있으며 이 유전자들 덕분에 유행성 균주들은 자연 상태의 다른 V. cholerae 균주와의 경쟁에서 이길 수 있었을 것입니다. 함께 알아봅시다.

  • [PAN-GENOME]→[Identify Differentially Present POGs/Pathways]로 이동합니다.
  • 중앙 패널에서 모든 유행성 균주(제6기, 제7기 대유행에 포함되는)를 선택해서 왼쪽 패널(Group #1)로 옮깁니다. 6기 및 7기를 제외한 나머지 균주들은 오른쪽 패널(Group #2)로 이동시킨 후 [Fisher’s exact test]을 클릭합니다. 이 통계 분석법은 연구자가 충분한 개수의 유전체를 가지고 있을 때 이용할 수 있습니다. 만약 비교할 유전체가 별로 없다면, Differentially present POG를 찾기 위한 단순 알고리즘인 [Exclusive Or (XOR)]을 대신 사용할 수 있습니다.

after selecting DPG

BIOiPLUG는 모든 단백질 코딩 유전자(CDS)를 서로 비교한 결과를 이용해서, 중복되지 않는 상동 CDS 세트(다른 말로 Pan-genome Orthologous Group (POG))를 생성하며, 이 POG들이 모여서 하나의 그룹을 대표하는 “pan-genome”을 생성하게 됩니다.

BIOiPLUG는 pan-genome과 관련된 두 개의 분석 결과를 제공합니다:

  • Differentially present orthologous groups: 두 유전체 그룹 간의 비교를 통해 유의미한 유전자들을 찾습니다.
  • Pathway enrichment: 이러한 orthologous group이 발견되면 이 정보를 이용하여, 상동 유전자가 서로 다르게 존재하는 KEGG pathway를 찾습니다.

[Differentially present orthologous groups] 탭을 클릭하면 아래와 같은 이미지가 나타날 것입니다.

  1. POG의 대표 CDS 이름
  2. Fisher’s exact test를 통해 얻어진 p-value. 값이 작을수록 유의합니다.
  3. 파란색 부분은 유전체에 해당 유전자가 있음을 의미합니다.
  4. 빨간색 부분은 유전체에 해당 유전자가 없음을 의미합니다.

피셔의 정확검정 (Fisher’s exact test)p-value를 사용하면 그룹간에 많은 차이를 보이는 유전자들을 찾아볼 수 있으며 유행성 균주에서만 특이적으로 발견되는 유전자들이 다수 존재하는 것을 알 수 있습니다. 몇 가지 주목할만한 DPG는 ctxA, ctxB (cholera toxin을 암호화하는 유전자) 와 tcpA (toxin과 함께 조절되는 pillin을 암호화하는 유전자) 입니다. 이 유전자들은 콜레라 병원성에 있어서 대단히 중요한 의미를 가지고 있다고 할 수 있습니다. 우리의 분석이 “무엇이 유행성 균주를 만드는가?”에 대한 해답을 발견한 것 같습니다.

우리가 생각해 보아야 하는 다음 질문은 “유행성 그룹과 비유행성 그룹 간에 기능적 측면(e.g. metabolism)의 차이가 존재하는가?” 입니다. 이에 대한 대답은 [Pathway enrichment] 탭을 사용하여 찾을 수 있습니다.

이번에도 p-value가 핵심적인 역할을 합니다. 통계적으로 가장 유의미한 pathway (KEGG database 가운데)로 확인된 것은 아래의 4개이며, 이들은 모두 0.05 미만의 p-value 값을 나타내고 있습니다.

  • Vibrio cholerae pathogenic cycle (질병 연관)
  • Vibrio cholerae infection (질병 연관)
  • Lipopolysaccharide biosynthesis (항원형 및 항체/면역반응 연관)
  • Pyruvate metabolism (대사 과정 연관)

연구가 더 필요한 “Pyruvate metabolism”을 제외하면, 이러한 분석이 분명히 설득력이 있다는 것을 알 수 있습니다.

얼마나 많은 유전자들이 존재하는가?

어떤 두 균주가 동종이라 할 지라도 완벽하게 유전자를 공유하지 않는다는 것은 이미 잘 알려진 사실입니다. 한 쪽의 균주에만 존재하는 유전자가 있는가 하면, 양 쪽 모두에 존재하는 유전자도 있습니다. 유전체 데이터의 이런 양상은 [Pan-genome]→[Venn Diagram]을 통해 확인해 볼 수 있습니다.

  • 4개의 유전체(N16961, TSY216, IEC224, 2010EL-1786)를 선택하고 [Draw] 버튼을 클릭합니다. 이 균주들은 유전체 염기서열 수준에서 매우 가깝게 연관되어 있습니다. (즉, OrthoANI 값이 높습니다.) 값을 직접 확인하고 싶다면, [Phylogenomics]→[OrthoANI]로 가서 ANI 값이 포함되어 있는 CSV 파일을 다운로드 하십시오.

  1. TSY216 균주는 매우 많은 유전자를 추가로 가지고 있다!

4개 균주가 3,413개의 유전자 (CDS) 를 공유하고 있으며 약 100개의 유전자가 단일 유전체에만 존재합니다. 놀랍게도 TSY216은 다른 균주에서는 찾아볼 수 없는 1,153개의 유전자를 추가로 가지고 있습니다. 이 균주들은 유전체 염기서열 수준에서 매우 비슷하기 때문에 이 정도로 많은 추가 유전자를 가지려면 상당한 크기의 염색체외 요소(extra-chromosomal elements)를 가지는 수밖에 없습니다. 이 경우, TSY216가 800,000 bp 크기의 매우 큰 플라스미드를 포함하고 있는 것으로 확인되었습니다(Okada et al., 2015).

전장 유전체 염기서열 정렬

유전체간에 ‘Pairwise’ 정렬을 통해 유전체의 구조적인 변화에 대한 충분한 정보를 얻을 수 있습니다. 앞서서 벤 다이어그램 분석을 통해 TSY216은 1,153개의 추가 유전자를 가지고 있는 것이 확인되었습니다. 아래 그림을 통해서 이 균주와 다른 Vibrio cholerae 7기 대유행 균주들 사이의 구조적인 차이점을 확인할 수 있습니다.

  • [PAIRWISE ANALYSIS]→[Whole genome alignment using NUCMER]로 이동합니다.
  • Reference genome” 으로 N16961균주를 선택합니다.
  • ”Query genome”으로 TSY216균주를 선택하고, [Run] 버튼을 클릭합니다.

  1. TSY216의 추가 유전자들(megaplasmid)은 reference로 비교된 N16961균주와 유사성이 없습니다.

TSY216균주에 있는 1,034개의 고유한 유전자는 최근에 획득한 megaplasmid로부터 얻어진 것이 분명해 보입니다.

유전자 유/무에 대한 분석

유전자가 각각의 유전체에 있는지 없는지에 따른 유전자 조성에 대한 정보는 몇가지 다른 방식으로 분석에 사용될 수 있습니다. 이 정보는 또한 계통수에 반영되기 때문에 계층적 클러스터링을 통해 dendrogram과 heat map으로 표현될 수 있습니다.

  • [Pan-genome]→[Gene presence/absence analysis]로 이동합니다.
  1. TSY216 균주만 가지고 있는 많은 유전자가 존재하는 것이 보입니다. 그러므로 UPGMA dendrogram을 그리면 특이적으로 바깥쪽에 위치하게 될 것입니다.

아래 그림은 구체적인 유전자 정보 없이 동일하게 그린 dendrogram입니다.

  1. TSY216균주가 많은 추가 유전자 때문에 클러스터링 분석에서 특이적으로 바깥쪽에 위치하게 된 것을 확인할 수 있습니다.

 

Differential gene (CDS) information (TSY216의 1,034개의 고유한 유전자와 같은 singleton 제외) 을 사용하면 TSY216 균주는 적절한 이웃들 사이의 위치로 들어가게 됩니다.

  1. TSY216균주에서 추가 유전자들(megaplasmid)을 제외하면 유전자 조성을 비교해 볼때 다른 7기 대유행 균주들과 매우 유사하다는 것을 확인할 수 있습니다.

수평적 유전자 이동(Horizontal gene trasnfer) 추적

박테리아 세상에서, 유전자는 여러 개의 유전자로 구성된 패키지 단위로 이동하곤 합니다. 이것을 “유전자 클러스터 (gene cluster)” 또는 더 흔한 말로 “유전자 섬 (genomic island)”이라고 부릅니다. 유전체에서 이러한 부분을 찾는 가장 좋은 방법은 pairwise gene content matrix를 생성하는 것입니다. 예제 데이터 셋에는 16개의 균주가 존재하므로, 각각을 1대 1로 (16 x 16)으로 비교하여 ortholog를 찾아낼 수 있습니다. 만일 ortholog를 찾아내는 방법에 대해 더 알고 싶다면, 여기를 살펴보시기 바랍니다.

이제 우리는 reference genome과 reference genome을 제외한 모든 다른 유전체를 비교해서 유전자들의 존재 여부를 확인할 것입니다. O 항원 합성에 관여하여 균주를 O 항원형으로 만드는 유전체 영역에 초점을 맞추도록 하겠습니다.

  1. [PAIRWISE ANALYSIS]→[Pairwise Ortholog Matrix]로 이동합니다.
  2. “N16961”을 reference genome으로 선택한 후, [RBH(Protein)] type으로 결과를 확인합니다. RBH는 ‘Reciprocal Best Hit”을 의미하며, 자세한 내용은 여기에 설명되어 있습니다.

  1. “Reference genome”으로 N16961을 균주를 선택합니다. 새로운 reference genome을 적용하기 위해서는 [Apply]를 클릭해야 합니다.
  2. 스크롤을 이용하여 GCF_000006745.1_00231 위치로 이동합니다. 마우스 커서를 GCF_000006745.1_00231의 사각형 위로 움직이면 위치 정보를 나타내는 툴팁이 뜰 것입니다. (1:231:24514; 이 좌표는 1번 컨티그, 231번 유전자, 24514 bp 위치를 뜻합니다.)
  3. 커서를 MO10(O139) 균주의 1_00231를 가리키는 사각형 위로 움직입니다. “23:345:89195”라는 좌표가 뜰 것이며, 이 ortholog가 유전체의 23번 컨티그, 345번째 유전자, 89195bp에 위치한다는 것을 알려줍니다.
  4. 마찬가지로 이 사각형은 N16961 균주의 “1:261:277002”에 위치한 유전자임을 의미합니다.
  5. 이 사각형은 MO10 균주의 “23:382:131391”에 위치한 유전자임을 의미합니다.

위의 ortholog matrix에서 O1 항원형의 균주는 대부분의 유전자를 공유하고 있다는 걸 알 수 있습니다. 그러나 O139 균주(M10와 4260B)에는 상당수의 유전자/ortholog가 누락되어 있는 것으로 보입니다(붉은색/주황색으로 표시됨). N16961의 261-231 영역(d에서 b까지)에는 30개의 유전자가 있습니다. 대조적으로 MO10 균주에서 여기에 대응하는 영역(382-345; e에서 c까지)을 확인해보면 더 많은 37개의 유전자가 포함됩니다. 그러나 N16961을 reference로 하는 이 matrix에서는 유전체 구성이 확인되지 않습니다.

이제 reference를 MO10으로 바꾸어서 MO10 유전체에 어떤 일이 일어나는지 시험해 보도록 하겠습니다. 위의 매트릭스에서는 영역의 시작 부분 유전자가 3번 컨티그의 345번째 유전자에 위치해 있었습니다. [Apply] 버튼을 눌러서 MO10으로 reference를 바꾸고 23번 컨티그의 GCF_000152425.1_00345. 위치로 이동합시다.

아래의 스크린샷을 보면 큰 유전자 세트가 O1 항원형이 존재하는 영역에 삽입되어 있는 것을 쉽게 확인할 수 있습니다. O1 유전자의 대부분이 O139 균주들에 존재하지 않기 때문에 우리는 O139의 이 영역이 (결손 또는 삽입으로) 다른 유전자로 대체되었을 것이라는 가설을 세울 수 있습니다. 이 과정은 퍼즐을 푸는 것과 같습니다.

  1. Reference genome이 MO10 (O139 항원형) 입니다.
  2. O139 균주에만 있는 많은 유전자들이 O1 항원형에 특이적인 영역을 대체하였다는 것을 보여줍니다. 이와 같이 수평적 유전자 이동이 O1에서 O139로의 항원형의 변화를 일으켰다는 것을 확인할 수 있습니다.

만약 이 유전자들의 기능을 확인하고 싶다면 [Product] 컬럼을 보시기 바랍니다. 많은 유전자들이 lipopolysaccharide의 생합성 능력에 관여한다는 것을 알 수 있습니다.

 

유전체와 유전자의 계통학적 이력

BIOiPLUG CG 서비스에서 주어진 세트에 포함된 유전체들의 모든 유전자 (CDS)는 비중복 상동 유전자 그룹 (non-reductant orthologous groups) 으로 묶이게 됩니다 (pan-genome과 core-genome이 어떻게 만들어지는지 보려면 여기를 확인). 이제 유전체와 단일 유전자에 기초한 phylogenetic tree를 그려 보도록 하겠습니다.

유전체 상에서 O 항원을 암호화하는 부분의  “ortholog matrix”를 아래에서 다시 보겠습니다(reference=N16961 O1 El Tor).

  1. “Reference genome”은 N16961로 설정되어 있습니다.
  2. GCF_000006745.1_00227라는 CDS(유전자 이름=VC0236 / waaF,rfaF) 모든 유전체에 존재합니다.

우리는 이 부분, O lipopolysaccharide 항원을 만들기 위해 필요한 ADP-heptose-LPS heptosyltransferase 2를 암호화하는 VC0236 / waaF,rfaF 유전자를 이용해 보도록 하겠습니다. 이 유전자는 모든 16개 유전체에서 발견됩니다. 16개 유전체로부터 모든 유전자를 찾아보려면 [PAN-GENOME]→[Browse Pan-genome Orthologous Groups (POGs)]로 이동해 주십시오. 모든 유전자들이 pan-genome으로 정리된 모습을 볼 수 있습니다.

  1. 빠른 검색을 위해서 “waaF”를 입력합니다.
  2. 해당 POG의 16개 유전체에서의 유전자 정보를 보기 위해 숫자 부분을 클릭합니다.

이제 POG에 포함되는 유전자(waaF 유전자로 명명)들을 볼 수 있는 페이지가 열렸습니다.

모든 waaF 유전자의 DNA 염기서열을 얻으려면 [Download DNA] 버튼을 클릭합니다. MEGA 프로그램이 설치되어 있다면 다운로드된 FASTA 형식의 파일을 이 프로그램으로 열 수 있을 것입니다. Multiple alignment와 phylogenetic analysis를 수행하고 나면 waaF 유전자의 phylogenetic tree를 얻을 수 있게 됩니다. (만일 MEGA7을 사용하고 계신 경우, 아래에 있는 순서대로 실행 하시기 바랍니다.)

  1. [Alignment]→[Align by ClustalW]를 눌러 기본 옵션으로 “multiple alignment”를 수행합니다.
  2. [Data]→[Phylogenetic analysis]를 선택해서 “phylogenetic analysis”로 이동합니다.
  3. MEGA의 메인 윈도우를 선택하고 [Phylogeny]→[Construct Maximum Likelihood Tree] 를 눌러 ML tree를 생성합니다.

아래 그림은 균주의 전장 유전체를 사용하여 구분한 계통 발생(phylogeny)과 병원성 유전자(이번 경우 waaF)으로 구분한 계통 발생 간의 차이를 설명하고 있습니다. 결과를 보면, 독소 비생성 O1 El Tor 균주인 LMA3984-4와 유행성 균주들은 유전체 수준에서는 유사하지 않지만 waaF 유전자를 비교해 보면 매우 높은 유사성을 나타냅니다. 이로 미루어 보아, 전장 유전체를 사용하여 계통발생을 확인을 하는 것이 필요하다는 것을 알 수 있습니다.

gene phylogeny3.PNG

마치며

여기까지 설명된 것은 우리가 BIOiPLUG의 comparative genomics를 통해 제공할 수 있는 많은 기능의 일부분에 불과합니다. 더 구체적인 매뉴얼과 훌륭한 스토리가 계속 제공될 예정이며 Chun et al. (2009)에서 많은 분석 방법을 소개하고 있으니 만일 비교유전체학 방법론과 Vibrio cholera 에 대해 더 알아보고 싶다면 이 논문을 확인해 보시기 바랍니다.

References

  1. Okada, K. et al. Characterization of 3 Megabase-Sized Circular Replicons from Vibrio cholerae. Emerg Infect Dis 21, 1262-1263 (2015).
  2. Chun, J. et al. Comparative genomics reveals mechanism for short-term and long-term clonal transitions in pandemic Vibrio cholerae. Proc Natl Acad Sci U S A 106, 15442-15447 (2009).

알림

이 tutorial은 Jon Jongsik Chun (Seoul National Univ/ChunLab, Inc)과 Suyeon Hong (Yale Univ)에 의해 작성되었습니다.

마지막 갱신일 July 2nd, 2017.