서 론
미생물법의학(Microbial forensics)은 생물 범죄 및 누출 예방 차원에서 생물학작용제의 기본 특성을 파악하고 출처를 추정하는 학문이다. 2001년 탄저균(Bacillus anthracis)을 이용한 테러 이후 미생물법의학 연구의 필요성이 크게 증가하였으며 특히 해당 종에 대한 계통 및 발생 추적이 활발히 이루어져 왔다[1]. USAMRIID (United States Army Medical Research Institute of Infectious Diseases)는 테러에 사용된 생물학작용제가 Ames 아류형(Subtype)으로 포자 형성 유전자의 부위의 인위적 변이가 있음을 밝혀냈으며, 2016년에는 1979년 Sverdlovsk 유출 사건의 샘플을 재분석하여 유전자 조작의 증거가 없음이 밝혀졌다[2].
미생물법의학은 과거 DNA 지문 분석(DNA fingerprinting)이나 다중 유전자좌 가변수일렬반복 분석(Multiple Locus Variable number of tandem repeat Analysis, MLVA) 등 짧은 서열을 중심으로 이루어졌다[3]. 이는 Sanger sequencing을 통해 얻을 수 있는 서열의 길이가 최대 1 Kb이며, megabase(Mb)당 $500의 비용이 소요됐기 때문이다[4]. 차세대염기서열분석법(Next Generation Sequencing, NGS)은 대량으로 전체 유전체(Whole-Genome Sequence, WGS) 정보를 얻음으로써 소요 비용과 시간을 혁신적으로 단축했으며, 전통적인 분석을 in silico 수준에서 가능케 했다. Illumina 사의 short read 기반 Genome Analyzer의 경우 최대 17 gigabases(Gb)의 서열 확보가 가능하며 Mb당 소요 비용은 $6 수준으로 감소하였다[4].
k-mer 분석법은 in silico 분석 중 하나로, 전체 서열을 길이가 k인 단편 서열로 나눈 후 자카드 계수(Jaccard similarity index)를 계산하여 전체 유전체 간 Distance를 판단하는 방법이다[5]. 일반적인 정렬 기법의 경우 서열 길이에 따라 계산 시간이 기하급수적으로 늘어나지만, k-mer 분석은 연구자의 k 값과 prefix(단편 서열의 첫 부분) 선택에 따라 소요 시간 단축이 가능하다. 박테리아 유전자에서 서열 ATGAC 는 코딩영역에 분포해 있으며, 기존의 문헌에서 다른 서열에 비해 높은 Average Nucleotide Identity(ANI)와의 상관계수를 가져 유의미한 그룹 분리가 가능함을 입증하였다[6]. k-mer 분석법은 이와 같은 prefix를 가진 단편 서열로 분석대상을 제한함으로써 소요 시간을 줄일 수 있으며, 넓은 분해능을 가져 역 단계(superkingdom domain; 고균역, 세균 역, 진핵생물 역)부터 Bacillus 속의 계통 추정까지 가능하다고 알려져 있다[5,7].
k-mer 기반 분석은 연산량을 조정할 수 있고 넓은 범위의 분류군을 구분하기 좋지만, 삽입, 삭제 및 재배열과 같은 변이가 유사도로만 판단되는 단점이 있다. 단일염기다형성(Single Nucleotide Polymorphism, SNP)은 유전적 다양성을 나타내는 또 다른 척도 중 하나로 대상 집단 내 통상 1 % 이상의 빈도로 발생하는 염기서열 변이이다[8]. 참조서열에 샘플 서열을 정렬한 후 Nucleotide Difference(ND) 값을 계산하여 신뢰구간 안에 있는 SNP를 판별한 후 Neighbor-joining 알고리즘을 통해 트리를 완성한다[9]. 전체 유전체 SNP 분석(whole genome SNP analysis; wgSNP)은 정렬 및 매핑 과정에 많은 소요 시간을 요구하지만, 진화적 패턴을 감지하고 특정 종에 대하여 계통군(clades)이나 하위분기군(subclade) 수준의 정밀한 분류가 가능하다[9,10].
본 논문에서는 젤상 고형 배지에 혼합되어있는 B. anthracis의 핵산을 분리하여 NGS 방법으로 전장유전체서열을 획득했다. NCBI(National Center for Biotechnology Information) 데이터베이스로부터 획득한 844개의 B. anthracis 서열 및 964개의 NGS 미가공 데이터(Sequence Read Archive, SRA)를 이용한 k-mer 분석 결과, 미지의 핵산은 유럽에서 2000년부터 2013년 사이 헤로인 투여 후 집단 감염된 Injectional Anthrax (IA)에 속하는 것으로 확인되었다. k-mer 분석 결과를 바탕으로 NCBI 데이터베이스로부터 획득한 IA 그룹의 SRA 데이터 58개에 대한 wgSNP 분석을 진행한 결과, 미지의 핵산은 B. anthracis의 Ames Ancestor [txid: 261594]의 염색체에서 236078, 507075, 3770887이 각각 G, A, T로 치환된 특징을 갖는 A4564임을 특정했다.
실험방법
분석 시료는 RKI(Robert Koch Institute)로부터 전달받았으며 0.85 % NaCl(#23-1-2)와 고형배지(#23-3-1)에 핵산이 혼합된 형태였다. Qiagen 사의 DNeasy Blood & Tissue Kits를 이용해 핵산추출을 하였고, 농도가 낮은 샘플은 효과적인 유전체 염기서열 확보를 위하여 TruPrime Whole Genome Amplification KitⓇ(4basebio™)을 이용하여 핵산 증폭을 수행하였다. 증폭된 핵산은 Covaris sonicator를 사용하여 Illumina sequencing에 맞는 크기(300 bp)로 분절화하였다. 그 후에는 Apollo324 라이브러리 자동화 시스템(TaKaRa)과 Illumina 호환 PrepX DNA Libray Kit(TaKaRa)를 사용하여 shotgun 라이브러리를 생성했다. 라이브러리 생성에는 Illumina sequencing에 적용되는 dual index adaptor을 사용하였고, DNA 서열 데이터는 Illumina 사의 양방향 서열 분석을 사용하여 생성했다. Illumina 사의 TruSeq Library Prep Kit를 이용해 시료 전처리 과정을 수행하였다. 모든 cleanup은 Beckman Coulter의 AMPure Xp bead를 이용했으며, 생성된 라이브러리에 대한 농도 검증은 일차적으로 Agilent Technologies의 TapeStation 4200과 D1000 ScreenTape, 이차적으로 KAPA Biosystems의 KAPA Library Quantification Kit와 Thermo Fisher Scientific 의 QuantStudio 6 Flex Real-time PCR system을 이용하였다. Illumina 사의 10 % PhiX control v3와 함께 NextSeq 500 benchtop platform을 통해 염기서열분석을 진행했다.
NCBI로부터 획득한 B. anthracis[txid: 1392] SRA 데이터 964개는 CLC genomics workbench 24.0을 이용하여 adaptor trimming 및 quality trimming(>0.05)을 진행했으며 Ames Ancestor를 참조서열로 삼아 공통서열(Consensus sequence)을 확보했다[12]. 마찬가지의 방법으로 미지 시료 Unknown Biological Sample(UBS) #23-1-2 (#23-1)과 #23-3-1(#23-3)로부터 분리된 핵산의 NGS 결과에서 공통서열을 확보했다. 이와 함께 조립(assembly)이 완료된 GenBank 데이터 844개를 사용하여 k-mer 분석을 수행했으며 k는 21 mer로, 참고 문헌에 따라 prefix는 ATGAC를 사용했다[6].
IA 그룹의 NGS 미가공 데이터 58개에 대한 wgSNP 분석은 Ames Ancestor를 참조서열로 삼아 기존의 문헌을 참고하여 수행되었으며, SNP calling은 SAMTools v0.1.18을 사용하여 정렬 후 위치당 coverage 10 %, prune distance 10을 기준값으로 잡아 필터링했다[9,13]. 매핑된 reads에 대해 z-score>=1.96 이상만 SNP로 인정했으며, UPGMA(Unweighted Pair Group Method with Arithmetic mean)을 사용해 bootstrap을 1000 회 시행했다. 각 시료의 wgSNP 분석에 사용된 유전체 내 염기서열의 위치(position)는 Ames Ancestor의 유전체를 기준으로 분석했다. 염기서열 분석 결과로부터 k-mer 분석 및 SNP 분석, 그리고 strain 특정의 전체 과정은 Fig. 1의 순서로 수행되었다.
결 과
미지시료에 포함된 B. anthracis 핵산의 출처를 밝히기 위하여 시료로부터 분리한 핵산의 NGS 결과 두 건(#23-1, #23-3)과 NCBI 데이터베이스에서 확보한 964개의 SRA 데이터로부터 Ames Ancestor를 참조서열로 하여 공통서열을 확보했다(Fig. 1). 이와 함께 NCBI GenBank로부터 획득한 조립된 서열 844개를 포함한 총 1,809개 데이터에 대해 k-mer 분석을 시행했다. 전체 데이터에 대하여 k = 21로 k-mer 분석을 진행했을 때 연산 자원 상의 문제로 Out of Memory (OOM) 에러가 발생했으며, prefix 값을 ATGAC로 지정해 이를 해결하였다.
k-mer 분석 결과 그룹 군은 7개의 그룹으로 분류되었으며 이는 기존에 알려진 아류형과 유사성을 보였다(Fig. 2). B. Anthracis 분류는 크게 clade A(∼90 %), clade B(∼10 %), clade C(∼1 %)로 이루어져 있으며[2], clade B에 속하는 B branch가 하나의 그룹으로 묶이는 것을 확인했다(Fig. 2, green). 이외의 6개 그룹은 모두 clade A에 속한 subclade로 이루어져 있었으며 여기에는 Ancient A, TransEurAsia(TEA), V770, Ames, Austrailia94(Aust94), Sterne 등이 속해있다[2]. 이들 중 TEA subclade는 가장 넓은 다형성(polytomy)을 가지며 Tsiankovskii, Injectional Anthrax(IA) 등 여러 그룹을 포함한다고 알려졌지만 k-mer 분석에서는 분산된 형태를 띠었다(Fig. 2, purple and ivory). 또한 sterne 그룹은 pXO가 결손된 백신 균주로 pXO를 가지고 있지 않은 모든 B. anthracis는 k-mer 분석 결과 sterne 그룹에 속하는 것을 관찰했다(Fig. 2, yellow). 미지 시료 #23-1과 #23-3은 이 중 IA 그룹에 속해있었다. IA 그룹은 2000년도 노르웨이에서 헤로인 사용자들로부터 처음 발병이 보고되었으며 공통적으로 사용한 주사기를 통해 감염 및 전파가 일어났다[14-16]. 가장 근접한 strain은 NCBI 데이터베이스로부터 획득하여 본 연구에서 공통서열을 추출한 SRR1999714(A4564)였으며, GenBank 로부터 획득한 IA 그룹의 조립된 서열과도 유사한 결과를 관찰했다(Fig. 2).
GCA로 시작하는 GenBank 조립서열은 기존 논문에서 SKESAv2.2를 통해 조립된 것이다. 실험을 통해 추출한 raw data를 가공하여 하나의 유전체 혹은 다량의 contig를 만드는 조립 과정은 다양한 알고리즘을 사용하지만, 조립 시 선택한 변수 조건(gap cost, reads 간 유사도 등)에 따라 tandem repeat이 반복되는 위치에서 조립오류(assembly error)를 동반하기 때문에 결과가 달라질 수 있다[18,19]. 그럼에도 SRR1999714로부터 유래된 GCA_29701425.1과 미지시료 #23-1/#23-3이 가장 가장 유사한 서열임을 k-mer 분석상에서 관찰할 수 있었다. 이는 첫째로 SKESA 알고리즘이 k = 21을 기본값으로 하는 raw reads의 k-mer 유사도 값을 이용하기 때문이며 prefix를 기반으로 한 k-mer 분석이 coding region에 주목함으로써 조립오류에 영향을 받지 않았으며, SKESA 알고리즘을 사용한 조립의 퀄리티가 충분히 좋았다고 해석할 수 있다. 해당 알고리즘은 SRA 데이터 및 약한 오염 시료에 대해 SPAdes 및 Megahit 같은 타 조립 알고리즘 보다 적은 조립오류를 보였다[17].
k-mer 분석법은 전체 서열을 길이가 k인 단편 서열로 나눈 후 자카드 계수(Jaccard similarity index)를 계산하는 방식이고, 본 연구에서는 연산량을 줄이기 위해 prefix를 ATGAC로 설정했다[5,6]. 따라서 k-mer 분석 결과 미지시료 #23-1/#23-3이 strain A4564와 가장 비슷함을 관찰했더라도 유전체에 대한 wgSNP 분석을 통해 염기서열을 다형성을 관찰해야만 미지시료의 출처 특정이 가능하다. 따라서 우리는 IA 그룹의 전장유전체서열을 획득한 프로젝트 PRJNA282011의 SRA 데이터 58개와 미지시료 #23-1/#23-3 두 시료에 대해 Ames Ansestor를 참조서열로 하여 wgSNP 분석을 수행했다.
SNP 분석 결과, IA 그룹 내에서 참조서열과 비교했을 때 6개 이하의 적은 SNP를 지닌 그룹과 16개 이상의 많은 SNP를 지닌 두 그룹으로 나누어지는 것을 확인했으며 각 서브 그룹(group-I, group-II) 내에서 SNP는 9개 이내였다(Fig. 3). 이는 IA 그룹이 2009∼ 2010년 스코틀랜드 발병 사건(group-I)과 2000년도 노르웨이로부터 유래한 2012-2013 유럽 전역의 발병 사건(group-II)으로 나누어진다는 기존의 연구에 부합했다[20].
Fig. 3.
Phylogenetic tree based on wgSNP alignment with 1000 bootstrap replications and SNP heatmap matrix of IA subclade
미지시료 #23-1과 #23-3이 속한 IA group-I의 48개 표본에 대한 심층분석 결과 확인된 SNP는 총 30개였으며 참조서열인 Ames Ancestor의 유전체를 기준으로 염기서열의 위치를 나타내었다(Fig. 4). A4575와 같은 SNP를 가진 14개의 strain은 A4575로 대표되었으며 마찬가지로 A4614와 같은 9개의 strain 또한 A4614로 표시했다. 미지 시료 #23-1/#23-3과 A4564는 236087의 아데닌(A)이 구아닌(G)으로, 507075의 시토신(C)이 아데닌(A)으로, 3770887의 시토신(C)이 트레오닌(T)으로 변이된 특징을 가지는 변이주로서 그 특성이 완벽하게 일치함에 따라 A4564로부터 유래한 시료임을 특정했다.
결 론
본 연구에서는 k-mer 분석 및 wgSNP 분석을 통하여 특정 매트릭스가 포함된 미량의 미지시료에 대한 법의학적 분석을 수행했다. 미지시료 #23-1/#23-3의 전장유전체서열 분석 결과로부터 공통서열을 추출하고 NCBI 데이터베이스로부터 확보한 1,808개의 B. anthracis 유전체 서열을 k-mer 분석을 통해 하위분기군 수준의 분류를 수행한 결과, 미지시료는 헤로인 사용자들에서 발병한 IA 그룹으로 분류됨을 확인했다. 이를 토대로 IA 그룹과 미지시료에 대한 wgSNP 분석 결과 미지 시료 #23-1과 #23-3은 A236087G, C507075A, C3770887T 변이를 특징으로 하는 B. anthracis strain A4564로부터 유래되었음을 밝혀냈다.
세균 출처의 법의학적 조사를 위해 NCBI 데이터베이스로부터 획득한 1,808개의 서열에 대한 wgSNP분석을 수행하는 것은 연산 자원상의 문제가 생길 가능성이 높을 뿐만 아니라, 제출된 공통서열의 조립 방식이나 분석에 사용한 NGS 결과의 품질에 따라 부정확한 SNP가 관찰될 가능성이 있다. 이러한 문제를 해결하기 위해 본 연구에서는 prefix를 설정한 k-mer 분석을 통해 미지시료의 출처를 IA 그룹의 58종으로 선별했으며, 추가적으로 wgSNP 분석을 수행하여 미생물 출처를 특정했다. 또한 고순도의 실험실 시료가 아닌 실제 법의학적 증거를 모사한 매트릭스가 포함된 미량의 시료로부터 핵산추출 및 NGS를 활용한 wgSNP 분석이 미생물법의학적 증거 획득에 유의미함을 증명했다.
PDF Links
PubReader
ePub Link
Full text via DOI
Download Citation
