차세대염기서열데이터 분석의 가장 첫 단계는 참조유전체(reference genome)에 대한 정렬(alignment)이며, 염기서열 정렬에 사용되는 참조유전체의 정확성과 완결성은 차세대염기서열분석의 유전변이 보고의 정확성에 중요한 영향을 미친다. 2001년 인간 참조유전체가 처음 보고된 이후, 현재까지 Genome Reference Consortium (GRC)에서 5번의 참조유전체 build 업데이트를 발표하였다. 현재 임상유전체분석에 주로 이용되는 참조유전체는 2009년에 발표된 GRC human 37 (GRCh37, hg19)이지만, 2013년에 발표된 GRC human 38 (GRCh38, hg38)이 염기서열 업데이트 측면(염기서열 개수, GC 함량, 엑손 개수)에서 hg19보다 정확성과 품질이 개선된 참조유전체이다[1,2].

GRCh38 임상 검증에 대한 최근 논문에 따르면, hg19보다 GRCh38을 적용할 때 더 정확하고 일관성 있는 유전체 데이터를 생성할 수 있다고 보고되었다[3,4]. 참조유전체 선택의 임상적 중요성을 가장 잘 보여주는 예는 Factor V Leiden allele (FVL) c.1601G>A p.(Arg534Gln)으로, FVL 동형접합형인 경우 정맥혈전의 위험도를 40–80배 증가시키는 변이이다[5,6]. 참조유전체 hg19를 정렬하는 파이프라인에서는 FVL 변이 대립유전자(alternative allele)가 아니라 정상 대립유전자(reference allele)가 검출(calling)되어, FVL을 놓치게 될 위험이 발생한다는 것이 잘 알려져 있다[2,3]. 또한 hg19를 사용할 경우 임상적으로 유의미한 NOTCH2 (Alagille syndrome 2; Hajdu-Cheney syndrome), KCNE1 (Jervell-Lange-Nielsen syndrome 2, long-QT syndrome 5), SIK1 (developmental and epileptic encephalopathy 30) 유전변이를 잘못 선별할 가능성이 높다는 보고도 있었다[3,7]. 이 외에도, 유전질환의 진단에 중요한 유전자 ADAMTSL2 (geleophysic dysplasia 1), NCF1 (chronic granulomatous disease 1), RPS17 (Diamond-Blackfan anemia 4)에서도 hg19와 GRCh38의 염기서열 정보에 차이가 있음이 보고된 바 있다[7]. 염기서열 정렬에 사용되는 참조유전체(hg 19 vs. GRCh38)의 선택은 임상유전체검사의 정도관리, 데이터 공유, 표준화 및 결과의 정확성 확보에 중요하다.

참조유전체 hg19를 GRCh38로 변경하여 적용하기 위해서는 데이터 분석 첫 단계부터 GRCh38로 재정렬을 시행(re-alignment to GRCh38)하거나 리프트오버 도구(liftover tools, 매핑 파일을 이용하여, 염기서열 어셈블리 간 위치를 hg19에서 GRCh38로 변환하는 도구)를 사용할 수 있다. 재정렬 작업을 통해 가장 정확한 변이 정보를 얻을 수 있지만, 시간과 비용에 대한 부담으로 인해 임상검사실에서 쉽게 선택할 수 있는 접근법이 아니다[8]. 최근 이러한 문제를 해결하기 위해 생물정보학적 리소스가 반드시 필요하지 않으면서 임상검사실에서도 쉽게 접근 가능한 리프트오버 도구가 개발되었다[8,9].

본 연구에서는 검증된 외부정도관리물질을 이용하여 참조유전체(hg19 vs. GRCh38) 선택에 따라 검출되는 유전변이 차이가 발생하는지 평가하고자 한다. 또한 리프트오버 도구 CrossMap, UCSCliftOver, NCBIRemap의 변환 정확도를 비교 평가하고자 한다.

대한진단검사정도관리협회(‘대한임상검사정도관리협회’에서 변경)로부터 4개의 FASTQ 파일(2022년 신빙도조사 GNG-22-01, GNG-22-02, GNG-22-04, GNG-22-05; Illumina, San Diego, CA, USA)을 획득하여 bwa-0.7.17(참조유전체 hg19: hs37d5.fa와 hg38: hg38.fa 각각 적용), picard-tools-2.19.0, GenomeAnalysisTK-4.2.0.0를 이용하여 binary alignment map과 variant call format (VCF) 파일을 생성하였다(Fig. 1). FASTQ 파일에는 대한진단검사정도관리협회로부터 사전 검증된 변이 144개(hg19에 정렬할 때)가 모두 포함되어 있다(단일염기변이: n=130, 삽입/중복/결손: n=14). 검증된 변이 144개의 명명은 2022년 신빙도조사사업 공통보고서를 참조하였다(http://eqas.keqas.org/). 리프트오버 도구를 적용하기 위해 VCF로부터 변이 정보를 추출하여 browser extensible data 파일을 생성하였다. 리프트오버 도구는 CrossMap (http://asia.ensembl.org/Homo_sapiens/Tools/AssemblyConverter?db=core), NCBIRemap (https://www.ncbi.nlm.nih.gov/genome/tools/remap), UCSCliftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver)을 이용하였다. FASTQ 파일에 대해 참조유전체 GRCh38을 정렬하여 분석한 유전변이를 정답(ground-truth)으로 하여 리프트오버 도구를 이용하여 hg19에서 GRCh38로 참조유전체를 변경 적용하여 검출한 유전변이를 비교하여 리프트오버 도구의 변환 정확도를 비교하였다(Fig. 1).

Figure 1. Overall workflow. The variants detected based on the realignment to human genome 19 (hg19) were compared to the variants detected based on the realignment to Genome Reference Consortium human 38 (GRCh38). The converted variants from hg19 to GRCh38 using liftover tools were compared to those detected based on the realignment to GRCh38. Abbreviations: BAM, binary alignment map; VCF, variant call format; BED, browser extensible data.

대한진단검사정도관리협회로부터 정보가 검증되지 않은 변이의 경우, 염기서열 깊이(depth of coverage)가 30× 미만인 경우와 예측 전사체(predicted messenger RNA) 및 비암호화 전사체(non-protein-coding transcript)에 있는 변이는 분석대상에서 제외하였고, 암호화 전사체(protein-coding transcript)에 있는 변이만을 대상으로 하였다. 특히 검증된 변이 144개 각각에 대해 참조유전체(hg19 vs. GRCh38) 선택에 따른 변이 검출에 차이가 발생하는지 평가하였다. snpEff (v4_3_GRCh37.p13.RefSeq)를 이용하여 Human Genome Variation Society 권장사항에 따라 염기변이를 기술하여 비교하였다. 변이의 형태는 Integrative Genomics Viewer (IGV)를 통해 시각적으로 확인하였다.

참조유전체 hg19에 정렬한 분석 파이프라인을 통해 사전 검증된 변이 144개를 포함하여 10,966개의 변이를 생성하였다. 참조유전체 GRCh38에 정렬한 분석 파이프라인을 통해서는 검증된 변이 141개를 포함하여 10,725개의 변이를 생성하였다. 대한진단검사정도관리협회로부터 사전 검증된 변이 144개 중에서, 제외된 변이 3개는 NM_000219.6(KCNE1): c.112A>G (p.Ser38Gly), NM_005429.5(VEGFC): c.1256_1258delTGA (p.Met419del), NM_004415.4(DSP):c.741T>G (p.Ala247=)로 확인되었다(Fig. 2). 3개를 제외한 141개 변이는 참조유전체 GRCh38 변경 적용에 따라 유전변이 종류 및 해석의 차이가 발생하지 않았다.

Figure 2. Comparison of alignment to hg19, human genome 19 (hg19) and Genome Reference Consortium human 38 (GRCh38). Alignment to hg19 (Left) and GRCh38 (Right). Three variants, which were not detected in the bioinformatics pipeline using GRCh38: (A) NM_000219.6(KCNE1): c.112A>G p.Ser38Gly; (B) NM_005429.5(VEGFC): c.1256_1258delTGA; and (C) NM_004415.4(DSP):c.741T>G (p.Ala247=).

참조유전체 GRCh38에 정렬한 파이프라인을 통해서 얻은 10,725개의 변이를 리프트오버 도구를 통해 얻은 변이와 비교했을 때, CrossMap과 UCSCliftOver는 전체 변이의 변환에 모두 성공한 반면, NCBIRemap은 1개 변이 NM_007075.3(WDR45): c.-18+2540delA의 변환에 성공하지 못했다(Fig. 3).

Figure 3. Variants that failed conversion by the use of the NCBIRemap. The non-converted variant of NM_ 007075.3(WDR45): c.-18+2540delA by using the NCBIRemap was located in the homopolymer region.

최근 참조유전체 사용 현황에 대한 국외 조사에 따르면, 28곳의 검사실 중 2곳의 임상검사실에서만 GRCh38 참조유전체 적용 파이프라인으로 변경 적용하였다[7]. 절반 이상의 검사실에서 GRCh38로 변경 적용 계획이 없다고 응답하였는데, hg19를 고수하는 주된 이유는 시간과 비용, 인력이 충분하지 않은 반면, 당장은 GRCh38로의 변경에 따른 이점이 크지 않을 것으로 생각하기 때문이라고 응답하였다[7]. 한편, 국내 임상유전체검사실 생물정보학적 파이프라인에 대한 설문에 따르면, 대부분의 검사실에서 hg19를 참조유전체로 사용하여 유전체분석을 진행하고 있으며, 한 임상검사실에서만 hg19가 아닌 GRCh38을 참조유전체로 이용한다는 단면적 현황파악에 대한 발표가 있었다[10]. 하지만 최근 국내에서도 국가바이오빅데이터구축사업(https://www.kobic.re.kr/ngp/process)과 같은 국가지원프로젝트에서는 전장유전체분석에 hg19가 아닌 GRCh38를 참조유전체로 사용할 것을 요청하고 있어 국내 검사실에서도 GRCh38 참조유전체의 임상 도입을 위한 준비가 필요할 것으로 보인다.

현재 사용 중인 참조유전체인 hg19를 지속적으로 사용할 경우, 또는 GRCh38로의 변경 적용하는 경우에는 각각의 장단점이 있다. 예를 들어, hg19를 지속적으로 사용할 경우, 기존 분석과의 일관성이 유지된다는 점에서 안정성이 있고, 시간과 비용 및 분석파이프라인에 대한 추가 투자가 필요하지 않다는 장점이 있다[7]. 하지만 hg19를 지속적으로 사용할 경우, 임상적으로 유의미한 변이를 놓칠 위험과 부정확한 변이 해석의 위험 가능성에 대한 연구논문이 발표되고 있으며, 향후에는 외부에서 제공되는 리소스들에서 hg19 지원을 중지할 가능성이 있다는 단점이 있다[3,5,6]. 반면, GRCh38로의 변경 적용 시에는 참조유전체의 정확성과 완결성 개선에 따른 정확도 향상의 장점이 있으며, 유전 변이 해석에 자주 이용되는 외부 리소스 업데이트 시에도 해당 리소스를 지속적으로 사용 가능하다는 장점이 있다[2,7]. 하지만 GRCh38로의 변경 적용을 위한 생물정보학적분석 파이프라인의 검증이 필요하다는 단점이 있다[7].

본 연구에서는 GRCh38 참조유전체의 임상 도입을 위한 준비의 첫 단계로 참조유전체(hg19 vs. GRCh38) 선택에 따라 검출 여부가 달라지는 변이의 존재를 확인하였다. 검증된 변이 144개 중 3개에서 참조유전체 종류에 따라 검출 여부가 달라지는 현상이 발생하였으며, hg19에 정렬하는 파이프라인에서는 검출된 반면, GRCh38에 정렬한 파이프라인에서는 확인되지 않았다.

첫째, GRCh38를 참조유전체로 정렬한 경우 검출되지 못한 NM_000219.6(KCNE1): c.112A>G p.Ser38Gly 변이는 IGV를 통해 확인하였을 때, 정렬의 정확도(mapping quality)가 감소된 리드에 존재하였다(Fig. 2A). 그리고 리드의 일부는 KCNE1뿐만 아니라 KCNE1 유전자와 염기서열이 비슷한 것으로 알려진 KCNE1B로 정렬한 것을 확인할 수 있었다. KCNB1B는 hg19에는 포함하지 않은 유전자이며, 업데이트된 GRCh38에만 포함된 유전자이기 때문에, NM_000219.6(KCNE1): c.112A>G p.Ser38Gly 변이는 GRCh38을 이용한 분석 파이프라인에서는 검출하기 어려운 변이로 추정된다. 어떤 참조유전체를 사용하든 상관없이 정렬 정확도가 감소된 리드에 변이가 존재할 때, 변이 검출의 정확도가 감소하는 것은 잘 알려져 있다. 참조유전체 자체의 정확도가 증가하였음에도 불구하고, 해당 유전자와 염기서열이 유사한 다른 유전자의 영향으로 인해 유전변이 검출의 정확도가 감소된 경우이다. 따라서 더 정확하고 업데이트된 참조유전체를 사용한다고 해서 반드시 모든 종류의 변이 검출의 정확도가 증가하는 것은 아니라는 사실에 주의할 필요가 있다.

둘째, GRCh38에 정렬한 경우 검출되지 않은 NM_005429.5 (VEGFC): c.1256_1258delTGA 변이의 경우, hg19에서는 존재했던 염기서열 TGA가 GRCh38에서는 존재하지 않은 것을 확인하였다(Fig. 2B). 즉 업데이트된 GRCh38 참조유전체에서는 염기서열 TGA가 삭제되었기 때문에, GRCh38에 정렬한 파이프라인에서는 변이로 검출되지 않아야 하는 변이에 속한다. 직접 재정렬이 아닌 리프트오버 도구를 이용할 때에도 이 변이는 NM_005429.5(VEGFC): c.1253_1255delTGA로 변환하여 검출되므로 변이 데이터 공유 시 소통 오류를 유발할 수 있어 주의가 필요하다. 실제 2022년 1차 국내 신빙도조사 결과보고서(http://eqas.keqas.org/)에 따르면, 이 변이에 대해 c.1253_1255del 또는 NM_005429.5(VEGFC):c.1253_1255delTGA 또는 NM_005429.5(VEGFC):c.1253_1255del로 보고한 검사실이 전체 22개 검사실 중에서 6개 검사실로 확인되었다. 신빙도조사 결과보고서에서는 해당 변이에 대해서는 평가하지 않았으나 hg19 참조유전체를 사용하는 검사실의 경우, NM_005429.5(VEGFC): c.1256_1258delTGA로 보고해야 하며, GRCh38 참조유전체를 사용하는 검사실의 경우에는 이 변이를 보고해서는 안 된다(http://eqas.keqas.org/). 외부정도관리 프로그램에서는 변이 보고 시 참조유전체 hg19를 이용하여 변이 기술을 시행할 것을 안내하였음에도 불구하고, 이와 같은 사무적 오류가 발생할 가능성이 있으므로 특별한 주의가 필요하다. 이 변이는 참조유전체 선택으로 인한 변이 기술 오류의 대표적인 사례임과 동시에, 참조유전체의 염기서열 구성을 고려하지 않고 리프트오버 도구를 이용할 때도 발생할 수 있는 오류이다. 검사실마다 일관성 있게 같은 종류의 참조유전체를 사용하지 않을 경우, 유전변이 데이터 공유 시 소통 오류를 유발할 수 있어 주의가 필요하다.

셋째, NM_004415.4(DSP):c.741T>G (p.Ala247=) 변이는 hg19 참조유전체에서는 변이 대립유전자가 GRCh38 참조유전체에서는 정상 대립유전자로 존재하기 때문에 GRCh38에 정렬하는 파이프라인에서는 변이로 검출되지 않은 경우에 속한다. 신빙도조사 결과보고서에서는 해당 변이에 대해서는 평가하지 않았으나 hg19 참조유전체를 사용하는 검사실의 경우, NM_004415.4(DSP):c.741T>G (p.Ala247=)으로 보고해야 하며, GRCh38 참조유전체를 사용하는 검사실의 경우에는 해당 변이를 보고해서는 안 된다.

본 연구에서는 리프트오버 도구 종류에 따라 참조유전체 GRCh38로 변경 성공률에 대해서도 평가하였으며, 리프트오버 도구 적용을 통해 참조유전체 GRCh38로 변경 성공률은 3개의 리프트오버 도구 모두에서 99.99% 이상으로 상당히 높은 비율을 보였다. NCBIRemap를 통한 변환에 실패한 NM_007075.3(WDR45): c.-18+2540delA 변이가 UCSCliftOver와 CrossMap을 이용한 변환에는 성공적이었으며, 세 개의 리프트오버 도구 모두에서 변환에 실패한 변이는 관찰되지 않았다. 생물정보학적 파이프라인 구축을 위한 추가 투자를 계획하기 어려운 임상검사실에서는 두 개 이상의 리프트오버 도구를 조합하여 사용하는 접근법이 GRCh38 변환 적용에 대한 실용적인 대안이 될 수 있을 것으로 생각된다. NCBIRemap를 통한 변환에 실패한 NM_007075.3(WDR45): c.-18+2540delA 변이는 IGV를 통해 확인하였을 때, 동형중합체가 밀집된 영역(homopolymer region)에서 발생한 변이로 확인되었다(Fig. 3). 동형중합체가 밀집된 영역의 경우 일반적으로 염기서열데이터분석에서 신뢰도가 낮은 영역으로, 리프트오버 도구를 이용한 변환에서도 이러한 영역은 별도로 취급될 필요가 있다. 리프트오버 도구를 통한 변환에 실패하는 유전변이에 대한 내부 데이터베이스를 구축하여 참조유전체 변경 작업 시 보완적으로 활용하는 것이 권장된다. 뿐만 아니라 리프트오버 도구를 이용하여 변환에 성공하더라도, 사용한 참조유전체의 염기서열 구성을 확인하여 위양성이 아닌지 확인하는 추가 절차가 필요할 것이다.

본 연구는 참조유전체 선택(hg19 vs. GRCh38)에 따른 변이 보고의 차이를 평가하기 위한 대상 변이의 개수가 제한적이라는 점에서(검증된 변이 144개) 한계가 있다. 하지만 전반적으로 참조유전체 GRCh38 적용에 대한 논의가 부재한 국내 상황에서 본 연구는 첫 번째 GRCh38 적용 결과라는 점에서 의의가 있다. 또한 본 연구에서는 실제 임상유전체 데이터를 직접 생성한 것이 아니라 외부정도관리 신빙도조사 프로그램에서 제공된 FASTQ 파일을 이용하여 분석하였다. 참조유전체 변경 도입의 필요성 및 리프트오버 도구의 정확성을 검증하기 위해, 추후 실제 임상유전체 데이터를 이용한 대규모 검증 연구가 필요하다.

요약하면, 본 연구에서는 국내에서 처음으로 참조유전체 선택(hg19 vs. GRCh38)에 따른 변이 보고의 차이를 평가하여 GRCh38의 임상 적용을 위한 기초 연구를 시행하였다. 뿐만 아니라 처음부터 GRCh38을 생물정보학적 파이프라인에 적용하여 분석한 유전변이를 정답으로 하여, 리프트오버 도구(CrossMap, UCSCliftOver, NCBIRemap)의 변환 정확성을 비교 평가하였다. 참조유전체 hg19를 이용하여 얻은 유전변이를 리프트오버 도구를 통해 변환할 때, 모든 변이가 GRCh38 정렬을 통한 결과와 일치하는 것은 아니라는 사실을 주지할 필요가 있다. 그럼에도 불구하고 리프트오버 도구는 생물정보학적 리소스에 대한 추가 투자를 계획하기 어려운 임상검사실에서는 선택할 수 있는 실용적인 대안 중 하나가 될 수 있다. 추후 여러 리프트오버 도구의 조합 사용과 참조유전체의 염기서열 영역에 대한 사전 검토를 통해 GRCh38 재정렬 과정의 비용을 절감하면서 유전변이 보고의 정확도 향상에 기여할 수 있을 것으로 기대된다.

이 연구는 대한진단검사정도관리협회의 2023년도 학술연구과제 연구비 지원으로 수행되었다(과제번호: 2023-10).