Constructing an ethnic-specific variant calling workflow based on a systematic comparison of multiple pipelines

Abstract

Next Generation Sequencing (NGS) is an experimental method that can read genome at high speed with low cost, while the sanger sequencing method used for the Human Genome Project (HGP), which took 13 years. NGS amplifies fragmented DNA by PCR and reads a sufficient amount of sequence through sequencing. The bases read by the sequencer, called "reads" are analyzed through three steps. The first step is the Quality Control (QC), which shows the quality of the reads and filters out the low-quality reads. The second is the alignment, which maps the reads to the reference genome. In the final step, variant calling analysis is performed to extract bases that differ between the mapped reads and the reference genome. There have been developed many methods for the alignment and the variant call of human with different algorithm. Their performance, however, has been studied using the reference genome based on Caucasian and Caucasian data called “NA12878”. Considering genetic differences, questions have been raised whether the methods optimized for Caucasians to show the same results for East Asians genetic differences. In this study, I analyzed the performance difference of each method using Whole Genome Sequencing (WGS) data of NGS. The WGS data used for comparing are Caucasian (NA12878) and Korean data. In case of NA12878 data, the performance is compared with the known correct answer (gold standard), and for Korean data, performance is compared with the data of the Korean microarray chip. Among the various methods that were used in the previous report, I selected 2 alignment methods (BWA-mem, NovoAlign) and 4 variant call methods (GATK4, Strelka2, DeepVariant, Samtools) that are known to perform well. The analyzable 8 combinations were compared. As a result of analyzing NA12878, the run time required for variant call in Samtools was long, and the recall was relatively low. Comparing the required time, the run time for the combination of BWA-mem and Strelka2 was the shortest. For Single Nucleotide Polymorphism (SNP) calling, which was regarded as a criterion for evaluating performance, the run time for the combination of BWA-mem and Strelka2 was the highest, and for recall, the combination of NovoAlign and GATK4 was the highest. In INDEL calling, precision was highest in the BWA-mem and Strelka2 combination showed very high precision, and recall precision was highest in the NovoAlign and Deep Variant combination. Therefore, in the combination for analyzing Korean genome, the variant call was performed with the combination excluding Samtools. The Markduplicate acting between alignment and variant call caused an overall time difference. The Korean sequencing data did not perform the PCR step, so there was no significant difference when comparing the performance with and without the Markduplicate step. Therefore, if the PCR step has not been performed, Markduplicate step can be skipped. When comparing the variants resulting from each pipeline with the variants of chip data using Korean sequence data, there is no large difference in the number of matching variants for each pipeline. However, among the pipeline combinations, the number of variants that matched the NovoAlign and GATK4 combinations was the highest from the data of all Koreans. Conclusively, the BWA-mem and Strelka2 combination is fastest to call variants and the NovoAlign and GATK4 combination is highly concordance with chip data and call many variants. When using the Korean reference genome instead of the existing reference genome in looking only at the BWA-mem and Strelka2 combinations, the SNP recall and INDEL precision and recall values were low, and there were few variants that matched the chip data. Therefore, it took a result that even if using Korean sequence data, alignment and variant call using existing reference genome would perform better.|차세대 염기서열 분석법(Next Generation Sequencing, NGS)라고 불리는 시퀀싱 방법은 13년이 걸렸던 인간 게놈 프로젝트(Human Genome Project, HGP)에서 사용한 생어 시퀀싱(sanger sequencing) 기법과 달리 빠르고 저렴한 비용으로 인간 유전체를 읽을 수 있는 분석법이다. NGS는 조각 낸 DNA를 중합효소 연쇄 반응(PCR)을 통해 증폭시키고 충분한 양이 된 서열들을 시퀀싱 기법을 통해 읽게 된다. 읽은 염기서열을 “reads” 라고 부르며 크게 3단계의 분석 단계를 거친다. 읽은 서열의 퀄리티를 나타내고, 퀄리티가 높은 reads만을 선별하는 품질관리(Quality Control, QC) 단계와 선별된 reads를 참조 유전체에 맵핑 시켜주는 정렬(alignment) 단계, 그리고 맵핑 된 reads와 참조 유전체에서 차이가 나는 염기를 추출해주는 변이 검출(variant call) 단계로 분석이 진행된다. 정렬과 변이 검출 단계는 사용할 수 있는 많은 방법들이 존재하며, 이 방법들은 저마다 다른 알고리즘을 가지고 있다. 정렬과 변이 검출 단계에서 사용되는 참조 유전체는 서양인을 기반으로 만들어져 있으며, NGS 데이터 분석 방법의 성능에 대한 여러 연구들도 “NA12878” 서양인 데이터를 사용하였다. 하지만 서양인에게 최적화 되어있는 방법들이 유전적 차이가 존재하는 동양인에서도 같은 결과를 보일지에 대해 의문이 제기되었다. 본 연구에서는 NGS를 이용한 전장유전체 시퀀싱(Whole Genome Sequencing, WGS) 데이터를 이용하여 각 방법의 성능 차이에 대한 분석을 진행하였다. 비교에 사용된 WGS 데이터는 서양인(NA12878) 데이터와 한국인 데이터이며 NA12878 데이터의 경우에는 알려진 정답(gold standard)과 성능을 비교하였고, 한국인 데이터의 경우에는 한국인을 대상으로 한 마이크로어레이 칩 데이터를 사용하여 비교했다. 다양한 방법들 중 이전 보고에서 성능이 좋았던 2개의 정렬 방법(BWA-mem, NovoAlign)과 4개의 변이 검출 방법(GATK4, Strelka2, DeepVariant, Samtools)을 선택하여 8개의 분석 가능한 조합을 비교했다. NA12878을 분석한 결과 Samtools로 변이 검출을 했을 때 소요시간이 길었고, 재현율이 상대적으로 떨어졌다. 소요시간을 비교하면 BWA-mem과 Strelka2 조합의 소요시간이 가장 짧았다. 성능을 평가하는 기준으로 삼았던 정밀도와 재현율에서 단일 염기 다형성(Single Nucleotide Polymorphism, SNP)은 BWA-mem 과 Strelka2의 조합이 가장 높았고, 재현율은 NovoAlign과 GATK4의 조합이 가장 높았다. 염기의 삽입과 결실 돌연변이(Insertion Deletion, INDEL) 정밀도는 BWA-mem과 Strelka2 조합에서 가장 높았고, 재현율은 NovoAlign과 DeepVariant 조합에서 가장 높았다. 따라서 한국인 유전체를 분석하는 조합으로는 변이 검출 단계에서 Samtools를 제외한 조합으로 분석을 진행하였다. 정렬과 변이 검출 과정 사이에 중복되는 reads를 표시하는 단계(Markduplicate)에 의해 소요시간의 차이가 생겼다. 한국인 시퀀싱 데이터는 PCR 단계를 수행하지 않았기 때문에 Markduplicate 단계를 수행하지 않은 경우와 수행한 경우의 성능을 비교했을 때, 큰 차이가 나지 않았다. 따라서 PCR 단계를 수행하지 않는 경우에는 Markduplicate 과정을 하지 않아도 된다. 한국인 시퀀싱 데이터를 사용하여 각 파이프라인 결과로 나오는 변이와 칩 데이터의 변이를 비교했을 때, 파이프라인 별로 일치하는 변이의 개수 차이가 크게 나지 않았다. 하지만 파이프라인 조합들 중에서도 모든 한국인 데이터에서 NovoAlign과 GATK4의 조합이 일치하는 변이의 개수가 가장 많았다. 결론적으로 빠른 속도로 변이를 확인하고 싶다면 BWA-mem과 Strelka2 조합을 사용하고, 많은 변이를 검출하며 칩 데이터와의 일치성이 높은 파이프라인을 사용하고 싶다면 NovoAlign과 GATK4 조합을 사용하는 것이 좋다. 기존의 참조 유전체 대신 한국인 참조 유전체를 사용한 경우 BWA-mem과 Strelka2 조합만을 봤을 때, SNP의 재현율과 INDEL의 정밀도, 재현율 값이 낮았고, 칩 데이터와 일치되는 변이가 적었다. 따라서 한국인 시퀀싱 데이터를 사용하더라도 기존의 참조 유전체를 사용하여 정렬과 변이 검출을 하는 것이 더 성능이 좋다는 결론이 도출되었다.