PDF
Citation
Print
서 론
염기서열분석기술의 향상으로 인간참조유전체(human reference genome)의 완결성이 점차 증가되면서, Genome Reference Consortium (GRC)에서는 주기적으로 인간참조유전체 build 업데이트를 발표하고 있다[1]. 현재 임상유전체분석에 주로 이용되는 참조유전체는 GRC human 37 (GRCh37/human genome 19; hg19)이며, 이를 기반으로 염기서열을 정렬(alignment)하여 유전변이를 명명하고 해석하여 보고하고 있다. 하지만, GRCh37에 비해 정확성과 완결성이 개선된 GRC human 38 (GRCh38, hg38)을 적용할 때 임상유전체분석의 정확성을 증가시킬 수 있다는 최근 연구 결과가 다수 보고되고 있다[2, 3].
GRCh38을 기반으로 명명되는 유전변이를 보고하기 위해서는 원칙적으로 염기서열 정렬 단계부터 GRCh38을 이용한 생물정보학적 파이프라인에서 재분석을 시행해야 한다. 재정렬을 통한 재분석을 시행하는 경우 유전변이를 가장 정확하게 변환할 수 있으나, 이를 위해서는 분석파이프라인, 시간, 비용, 인력에 대한 투자가 요구된다[4]. 재정렬에 대한 비용효과적인 대안으로, 리프트오버 도구(liftover tools)가 개발되어 연구 목적으로 활용되고 있다[5, 6]. 리프트오버란 기존에 GRCh37에 매핑된 파일을 그대로 이용하여 GRCh37과 GRCh38 염기서열 어셈블리 상호 간 위치를 변환하는 방법이다. CrossMap (http://asia.ensembl.org/Homo_sapiens/Tools/AssemblyConverter?db=core) 또는 UCSC liftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver)는 이러한 목적으로 개발된 대표적인 리프트오버 도구 중 하나이다[6]. 그런데 이 도구들은 단일염기서열변이(single nucleotide variants, SNVs)의 변환에는 적합하지만, 삽입/결실/중복(insertion/deletion/duplication, indels)의 변환에는 적합하지 않아, 기존 연구에서는 단일염기서열변이의 변환에 대한 평가 결과를 주로 보고하였다[6, 7]. 뿐만 아니라, 기존 리프트오버 도구에서는 유전체 변환이 되지 않거나 잘못 변환되는 빈도가 높고, 이는 리프트오버의 임상 적용에 대한 한계로 보고되고 있다[6, 7]. 또한, 삽입/결실/중복 변이는 반복서열(repetitive sequence) 구간에 위치하는 경우가 많고, 반복서열 구간은 리프트오버의 정확성을 감소시키는 주된 요인으로 보고된 바 있다[7]. 삽입/결실/중복 변이에 대해 효과적으로 작동할 수 있는 리프트오버 도구(예를 들어, Genozip/DVCF, GenomeWarp)도 개발되었으나, 삽입/결실/중복 변이의 리프트오버의 정확성은 단일염기변이의 정확성보다 낮은 결과를 보였다[8, 9].
재료 및 방법
ClinVar 데이터베이스에서 GRCh37과 GRCh38에 정렬된 2개의 variant call format (VCF) 파일을 획득하였다(https://ftp.ncbi.nlm.nih.gov/pub/clinvar/). ClinVar 데이터베이스에서 제공된 변이 평가기준(review status)에 따라, 가이드라인에 등재된 유전변이(practice guideline designation; 4 gold stars), 전문가 패널에 의해 검토된 유전변이(reviewed by expert panel; 3 gold stars), 다수 보고된 유전변이(variants provided by multiple submitters with assertion criteria and evidence; 2 gold stars) 262,156개(GRCh38 기준, No. of SNVs=238,735, No. of indels/duplications=23,421)를 선별하였다. 23,421개의 삽입/결실/중복 변이에서 분자 기능에 영향이 없는(“molecular consequence: no sequence alteration”) 2개[NM_001378454.1(ALMS1):c.36_74=(p.Leu12_Glu25=), (ClinVar variation ID: 193377); NM_001105206.3(LAMA4):c.848_849=(p.Asp283=), (ClinVar variation ID: 177976)]를 제외하고, 삽입/결실/중복 변이 23,419개를 대상으로 분석하였다.
GRCh37에 정렬된 VCF 파일의 유전변이에 대해 BCFtools/liftover를 이용하여 참조유전체 GRCh38로 리프트오버를 수행하였다. UCSC chain file을 이용하여 리프트오버를 수행하였다(http://hgdownload.cse.ucsc.edu/goldenpath/hg18/liftOver/hg18ToHg38.over.chain.gz). GRCh38에 정렬된 VCF 파일에서 얻은 유전변이 정보를 정답(ground-truth)으로 하여, BCFtools/liftover를 통해 리프트오버된 변이의 변환 정확도를 비교하였다(GRCh38-realigned variants vs. liftovered variants from GRCh37 to GRCh38). ClinVar에서 제공된 ClinVar variation ID를 기준으로, 리프트오버된 변이와 재정렬된 변이를 일치시켜(match) 비교하였으며, 각각의 변이에 대해 염색체(chromosome), 위치(genomic position), 정상대립유전자(reference allele), 변이대립유전자(alternative allele) 정보를 비교하였다. 동일한 조건에서 CrossMap을 구동하여 BCFtools/liftover 변환 결과와 비교 평가하였다.
Ensembl Variant Effect Predictor (VEP release-106)를 이용하여 Human Genome Variation Society (HGVS) 권장사항에 따라 염기변이를 기술하여 비교하였다. Mutalyzer 3 (https://mutalyzer.nl/)를 이용하여 변이 명명의 정확도를 재확인하였다. 변이가 반복서열 구간에 위치하는지 확인하기 위해 UCSC의 RepeatMasker track information을 이용하였다.
결 과
GRCh38에 정렬된 VCF 파일에 포함된 분석 대상 삽입/결실/중복 변이 23,419개 중에서 8개의 변이(8/23,419, 0.03%)가 BCFtools/liftover를 통한 리프트오버 변이에서 확인되지 않았다: NM_001378454.1(ALMS1):c.1571CTC[1] (p.Pro525del), (ClinVar variation ID: 210124); NM_001267550.2(TTN):c.100766-10dup, (ClinVar variation ID: 196651); NM_000492.3(CFTR):c.1210-12T[9], (ClinVar variation ID: 161188); NM_006446.4(SLCO1B1):c.359+10[16], (ClinVar variation ID: 440285); NM_000527.5(LDLR):c.314-446_1187-386dup, (ClinVar variation ID: 251140); NM_006031.6(PCNT):c.8751+23dup, (ClinVar variation ID: 1256693); NM_000381.4(MID1):c.661-7dup, (ClinVar variation ID: 196261); NM_007325.5:c.-2dup, (ClinVar variation ID: 167157). BCFtools/liftover를 통해 변환되었으나, GRCh38에 정렬된 변이와 비교하였을 때 위치 및 대립유전자 정보에서 차 이를 보이는 변이는 20개(변환 성공률 23,399/23,411, 99.94%)가 관찰되었다(Table 1). 이 변이 중 16개 변이의 경우, 리프트오버와 재정렬 결과 변이의 위치는 모두 일치하였으며, 단 4개의 변이(Var1, Var10, Var16, Var20)에서는 위치와 대립유전자 정보 모두에서 차이를 보였다(Table 1). 반면, CrossMap을 통해 리프트오버를 수행한 결과, GRCh38에 정렬된 변이와 비교하였을 때 위치 및 대립유전자 정보에서 차이를 보이는 변이는 557개(2.38%, 557/23,411)로 확인되었다(Table 2).
반복서열 구간에 위치한 삽입/결실/중복 변이 1,646개 중 (1) NM_003172.4(SURF1):c.54+13TGCGGGG[5] (Var8, ClinVar variation ID: 496202, microsatellite, 14 bp), (2) NM_000527.5(LDLR):c.695-67_1586+371del ((Var14, ClinVar variation ID: 251409, large deletion(7,636 bp)), (3) NM_001291088.2(WDR87):c.6889GAGGAGGAGGAGGAAAGG[1] (p.2285EEEEER[3]), (Var16, ClinVar variation ID: 546322, microsatellite, 18 bp)를 제외한 대부분의 변이(99.82%, 1,643/1,646)는 BCFtools/liftover를 통해 성공적으로 변환되었다(Table 1). 반면, 반복서열 구간에 위치한 1,646개 중 98.18% (1,616/1,646)의 변이가 CrossMap을 통해 성공적으로 변환되었다. 또한, BCFtools/liftover를 통해 성공적으로 변환되지 않은 변이 중 NM_003172.4(SURF1):c.54+13TGCGGGG[5] (Var8, ClinVar variation ID:496202, microsatellite, 14 bp), (2) NM_000527.5(LDLR):c.695-67_1586+371del ((Var14, ClinVar variation ID: 251409, large deletion(7,636 bp))은 CrossMap을 통해 성공적으로 변환되었다.
고 찰
본 연구에서는 삽입/결실/중복 변이의 유전체 변환을 위한 리프트오버 도구로 BCFtools/liftover의 임상 적용 가능성을 평가하였다. 임상적으로 중요한 유전변이의 집합인 ClinVar 데이터베이스에 등록된 삽입/결실/중복 변이를 대상으로, BCFtools/liftover을 통해 GRCh37에서 GRCh38로 참조유전체를 변환한 변이(liftovered variants) 정보와 GRCh38에 정렬된 변이(aligned variants) 정보를 비교하였다. GRCh38에 정렬된 변이의 위치, 정상대립유전자, 변이대립유전자 정보를 비교하였을 때, BCFtools/liftover를 통한 삽입/결실/중복 변이의 변환 성공률은 99.94% (23,399/23,411)로 확인되었으며, 특히 반복서열 구간에 위치한 삽입/결실/중복 변이의 변환 성공률은 99.82% (1,643/1,646)로 확인되었다. 반복서열은 참조유전체 정렬뿐만 아니라 리프트오버의 정확성을 떨어뜨리는 주된 요인 중 하나이다. 본 연구에서 평가한 BCFtools/liftover는 2개의 microsatellite (Var8, Var16)와 7,636 bp deletion 변이(Var14)를 제외하고 반복서열 구간에 위치한 모든 삽입/결실/중복 변이의 유전체 변환을 성공적으로 수행하였다.
본 연구에서 CrossMap을 통한 삽입/결실/중복 변이의 변환 성공률은 97.62% (22,854/23,411), 반복서열 구간에 위치한 삽입/결실/중복 변이의 변환 성공률은 98.18% (1,616/1,646)로 확인되었다. CrossMap의 유전체 변환 정확도와 비교하였을 때, BCFtools/liftover를 통한 변환 성공률은 CrossMap을 통한 변환 성공률보다 개선된 것을 확인할 수 있었다[6, 10]. 1000 Genome Project variant callsets을 이용하여 리프트오버 도구 5개의 성능을 비교 평가한 이전 연구 결과에서도, CrossMap은 가장 낮은 변환 성공률을 보였으며, BCFtools/liftover는 가장 높은 변환 성공률을 보였다[10]. 하지만, Var8 (microsatellite)과 Var14 (7,636 bp large deletion)에서와 같이 CrossMap은 성공적으로 변환시키지만, BCFtools/liftover은 변환시키지 못한 변이도 있는 것을 본 연구에서 확인하였다. 추후 더 많은 형태의 유전변이 데이터를 이용하여 여러 리프트오버 도구의 정확도의 비교 평가를 시행하여, 리프트오버 도구별로 차이를 보이는 conversion-unstable regions를 공유하고, 여러 리프트오버 도구를 상호 보완적으로 활용한다면, 리프트오버를 통한 유전체 변환의 정확도를 향상시킬 수 있을 것으로 판단된다.
리프트오버와 재정렬 결과를 비교하였을 때, 대립유전자 정보에서 차이를 보인 변이 20개(20/23,411, 0.09%)는 allele swap (변이대립유전자가 정상대립유전자로 뒤바뀌는 경우 또는 새로운 정상대립유전자가 생성되는 경우)이 발생한 변이로 확인되었으며, 그중 19개 변이는 multiallele 변이(단일 유전변이가 리프트오버를 통해 2개 이상의 변이로 변환되는 경우)로 변환된 것을 확인하였다(Table 1). 예를 들어, NM_014249.4(NR2E3):c.951del (Var11, ClinVar variation ID: 445726)는, GRCh37 기준 15: 72105929 위치에서, GRCh38 기준 15: 71813589 위치에서 염기서열 C의 결실이 발생된 변이이다(정상대립유전자: AC, 변이대립유전자: A). 그런데, BCFtools/liftover를 이용하여 리프트오버를 수행하면, 정상대립유전자 A가 변이대립유전자 AC로 바뀌면서 염기서열 C의 중복이 발생되는 것으로 잘못 변환되었다(NM_014249.4(NR2E3):c.951dup). 리프트오버를 진행하는 동안 allele swap을 통해 유전체 변환의 정확도를 높이는 것은 기존 CrossMap과 같은 리프트오버 도구와 차별화되는 BCFtools/liftover의 중요한 특징 중 하나이지만, Var11의 경우와 같이 잘못 변환되는 경우가 있어 주의가 필요하다[10]. BCFtools/liftover를 통한 리프트오버 결과를 살펴보면, 해당 변이에서 allele swap이 발생했는지 여부를 확인할 수 있도록 추가 정보가 제공되고 있다. 따라서, BCFtools/liftover를 통한 리프트오버의 정확성을 검토할 때, allele swap 발생 정보를 flag로 활용할 수 있을 것으로 판단된다.
Multiallele 변이로 변환된 변이 19개 중 변이의 위치와 대립유전자 정보가 모두 일치하지 않은 변이(Var1, Var10, Var16, Var20)는 4개로 확인되었다. 이 4개를 제외하면, 리프트오버를 통해 변환된 변이 2개 중 하나의 변이는 재정렬을 통한 변이와 일치하는 것을 확인할 수 있었다. 예를 들어, Var2는 GRCh37에 정렬되었을 때 NM_000348.4(SRD5A2):c.81_91del로 명명된 변이로 BCFtools/liftover를 적용하면 NM_000348.4(SRD5A2):c.90del과 NM_000348.4(SRD5A2):c.81_94del과 같이 2개의 변이로 변환되었다. 이 변이는 GRCh38에 정렬하였을 때 NM_000348.4(SRD5A2):c.81_94del로 명명되었으며, 이는 리프트오버를 통해 변환된 2개의 변이 중 하나에 속한다. 해당 변이는 GRCh37에서 2: 31805877, GRCh38에서는 2: 31580807 위치에서 시작하는 변이로, 해당 유전변이가 위치하는 참조유전체 GRCh37과 GRCh38 염기서열을 비교하면, GRCh38에서 2: 31,580,813 부위에 염기서열 G가 추가되어 있다(정상대립유전자: CCGGAGGGCTTCGCG). GRCh38로 유전체를 변환하는 과정에서는 기존 GRCh37에서의 정상대립유전자(CCGGAGGCTTCGCG)까지 변이대립유전자로 간주하게 되고, 2개의 변이대립유전자(CCGGAGGCTTCGCG, C)를 갖는 것으로 변환되게 된다. 이런 경우는 GRCh37에서의 정상대립유전자 염기서열(CCGGAGGCTTCGCG)을 제외하고, 염기서열 C만을 변이대립유전자로 선택하면, 리프트오버를 통한 변환 결과와 재정렬을 통한 결과가 일치하는 것을 확인할 수 있었다. Multiallele 변이 등으로 변환되는 경우에는 추가적으로 변환 전후 참조유전체의 염기서열 구성을 확인하여 변이의 위치와 대립유전자 정보를 확인하는 절차를 통해 리프트오버의 정확성을 검토하는 것이 권장된다.
본 연구에서 반복서열 구간에 위치한 삽입/결실/중복 변이 1,646개 중 98.18% 변이가 CrossMap을 통해 성공적으로 변환된 반면, BCFtools/liftover를 통해서는 99.82% 변이가 성공적으로 변환된 결과는 상당히 흥미롭다. CrossMap과 UCSC liftover를 평가한 이전 연구에 따르면, conversion-unstable region의 특징을 참조유전체 어셈블리 간 간극(gaps), 참조유전체 build 간 contig 차이, segmental duplication을 포함한 반복서열 구간으로 보고한 바 있다[7]. 본 연구에서 위치와 대립유전자 정보 모두에서 차이를 보인 4개의 변이 중에서 Var16만이 반복서열 구간에 위치하였으며, 모두 segmental duplication도 아닌 것으로 확인되었다. BCFtools/liftover를 개발한 Broad Institute에서 conversion unstable region 정보를 제공하고 있지는 않으나, BCFtools/liftover의 conversion unstable region은 기존 리프트오버 도구와는 다른 특징이 있을 것으로 추정된다. BCFtools/liftover의 임상 적용을 위해서는 BCFtools/liftover에 의한 conversion-unstable region의 특징을 규명하기 위한 추가연구가 필요할 것으로 판단된다.
본 연구에서는 GRCh38에 정렬된 변이에는 존재하지만, BCFtools/liftover를 통한 변환된 변이 리스트에 존재하지 않는 변이가 8개로 확인되었다. 이 변이들은 모두 GRCh37에 정렬된 VCF 파일에는 존재하지 않는 변이로 확인되었다. 즉, 이 변이 검출에 대한 리프트오버의 위음성은 BCFtools/liftover의 한계 때문이 아니라, 참조유전체가 GRCh38로 업데이트되면서 기존 GRCh37에는 변이로 명명되지 않았던 염기서열이 새롭게 유전변이로 등장한 경우에 해당한다. 리프트오버 도구 적용과 상관없이, 참조유전체 차이로 인해 GRCh38에 정렬된 파이프라인에서는 유전변이로 명명되는 반면, GRCh37에 정렬된 파이프라인에서는 변이로 검출되지 않은 경우가 있다. 이와 반대로 GRCh37에 정렬된 파이프라인에서는 변이로 검출되지만, GRCh38에 정렬된 파이프라인에서 변이가 아닌 경우도 가능하다. 예를 들어, Factor V Leiden allele인 NM_000130.4(F5):c.1601G>A (p.Arg534Gln)는 GRCh37에 정렬하는 경우, 변이대립유전자가 아니라 정상대립유전자로 검출되는데, 이런 종류의 변이의 경우 리프트오버 도구의 성능을 개선한다고 해도 검출할 수 있는 변이가 아니다. 따라서 리프트오버 도구를 적용하기 전에, 참조유전체의 염기서열 구성에 따라 변이 여부가 달라질 수 있다는 것을 고려하여 리프트오버 도구의 성능을 평가할 필요가 있다.
본 연구에서는 최근 개발된 BCFtools/liftover를 이용하여 Clin-Var 데이터베이스에 등록된 삽입/결실/중복 변이의 유전체 변환을 평가하였다. BCFtools/liftover는 기존 리프트오버의 한계로 지적된 반복서열 구간에 위치한 삽입/결실/중복 변이의 유전체 변환을 성공적으로 수행하였다. 추후 BCFtools/liftover의 임상 적용을 위해서는, 실제 임상유전체 데이터를 이용한 대규모 검증 연구가 필요하다.
XML Download