NGS DNA-seq pipeline: GATK Best Practice Code – Part2. Bam to Vcf

지난 글에 이어서, 생성된 Bam 파일로부터 변이들을 읽어 들이고, haplotype call을 하는 두번째 파트의 code를 정리해보겠습니다. 아래 코드는 GATK 4.1.3 버젼을 기반으로 작성되었습니다. GATK 버젼에 따라서 조금씩 Tool과 명령어에 차이가 있습니다. 전반적인 흐름은 아래와 같습니다.

관련 포스팅 보기>

NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam

NGS 분석 파이프 라인의 이해: GATK Best Practice

[계속 Update 예정] 자주 쓰는 linux 명령어 및 프로그램 관련 자료

1.jpg

III. Germline short variant discovery : Bam to Vcf

1. HaplotypeCaller

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_haplotypecaller_HaplotypeCaller.php

해당 위치의 변이와 read 갯수등을 바탕으로 하여, haplotype을 생성해주는 과정입니다. 일반적으로 우리가 heterozygote인지 homozygote인지에 해당하는 변이의 zygosity를 이 과정에서 생성한다고 볼 수 있습니다. 앞선 과정에서 BQSR을 통해 최종적으로 생성된 Bam 파일을 이용하여, haplotype call을 진행합니다.

gatk -- java-options 'Xmx16g' 'Xms8g' HaplotypeCaller -R [hg19_reference.fa] -I [sample01.final.bam] -ERC GVCF [-L targets.interval_list] \ -O [sample01.vcf] --genotyping-mode DISCOVERY --standard-min-confidence-threshold-for-calling 30 -ploidy 2 --output-mode EMIT_VARIANTS_ONLY

-L option에는 타겟 영역에 해당하는 BED 파일을 넣어주면 해당 위치에서만 Call이 일어나게 됩니다. Whole genome이 아닌 Target panel이나 Exome panel의 경우, 타겟 유전자들로 구성된 BED 파일을 지정해줍니다.

2. GenomicsDBImport

https://github.com/GenomicsDB/GenomicsDB/wiki

이전의 GATK와 달라진 부분입니다. 예전에 여러개의 샘플로 구성된 VCF를 Genotype GVCF로 합쳤는데, 추후의 편의성 및 연산 속도를 고려하여 GenomicDB를 구축하는 과정이 추가되었습니다. 이 단계를 건너 뛰어도 좋지만, 일반적으로 이렇게 DB를 구축하고 나면 여러모로 좋다고 소개하고 있네요. 여러개의 샘플로 구성된 환자 코호트의 경우는 이렇게 DB를 구축하고 하나로 합쳐주는 과정이 있는게 좋을 것 같고, Single sample이라면 이 과정은 건너뛰어도 좋을 것 같습니다.

gatk -- java-options 'Xmx16g' 'Xms8g' GenomicsDBImport -V [vcf file_list] \ [-L targets.interval_list] --genomicsdb-workspace-path \ [/genomicDB] --merge-input-intervals true --tmp-dir=/tmp 

[vcf file_list] 부분에는 -V sample01.vcf -V sample02.vcf -V sample03.vcf 와 같이 하나의 cohort로 구성하고자 하는 여러개의 샘플을 하나로 쭉 나열해주면 됩니다.

3. GenotypeGVCFs

위에서 GenomicDB를 구축하고, 여러 개의 샘플로 구성된 VCF가 있었다면, 이들을 Cohort 단위로 묶어서 하나의 VCF로 만들어 주는 과정입니다.

gatk -- java-options 'Xmx16g' 'Xms8g' GenotypeGVCFs -R [hg19_reference.fa] -V gendb://genomicDB --tmp-dir=/tmp -O [cohort.vcf]

4. VariantRecalibrator, ApplyRecalibration

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_VariantRecalibrator.php

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_ApplyVQSR.php

Bam을 구성할 때, 각 염기별로 recalibration을 진행했던 것과 비슷하게 call된 변이들에서도 recalibration을 진행해서 QC를 하는 부분입니다. 머신러닝을 이용하여, 기존 DB로 부터 변이들을 학습 시키고 걸러내는 방법이라고 하는데, 이를 위해서 hapmap, omni, 100G, dbsnp database를 다운 받아야합니다. 최근에는 CNN (convolutional neural network)에 기반한 모델을 수립하여, 테스트 중이라고 하는데 아직까지는 Beta 버젼에 머물고 있어 기존 코드를 이용합니다. 추후 이 부분은 CNN 기반 Recalibration code로 바뀔 가능성이 있습니다.

gatk -- java-options 'Xmx16g' 'Xms8g' VariantRecalibrator -R [hg19_reference.fa] -V [cohort.vcf] \
--resource:hapmap,known=false,training=true,truth=trueprior=15.0 hapmap_3.3.hg19.sites.vcf.gz \
--resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg19.sites.vcf.gz \
--resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.vcf.gz \
--resource:dbsnp,known=true,training=false,truth=false,prior=2.0 Homo_sapiens_assembly19.dbsnp138.vcf.gz \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode SNP -O [cohort.recal.vcf] --tranches-file [output.tranches] --rscript-file [output.plots.R]

gatk ApplyVQSR -R [hg19_reference.fa] -V [cohort.vcf] -O [cohort.final.vcf] \
--truth-sensitivity-filter-level 99.0 --tranches-file [output.tranches] --recal-file [cohort.recal.vcf] -mode SNP 

위의 과정을 거치면, 최종적으로 변이들이 call 되어 분석이 가능해집니다. 마지막 파트는 Annotation 과정인데, 해당 코드는 마지막 Part 3 포스팅에서 다루도록 하겠습니다.

[Reference]

GATK Best Practice for Germline short variant discovery (SNPs + Indels)

: https://software.broadinstitute.org/gatk/best-practices/workflow?id=11145

How to run VQSR from GATK

https://software.broadinstitute.org/gatk/documentation/article?id=2805

NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam

GATK4 (Genome Analysis Tool Kit)로 넘어오면서, NGS 시퀀싱 분석을 위한 파이프라인이 많이 개선 및 간소화된 것 같습니다. 덕분에 저도 최근 연구실에 구축되어 있던 파이프라인도 새롭게 뜯어고쳤는데, 이 참에 전반적인 분석을 위한 코드를 정리해볼까 합니다. 추후에 비슷한 파이프 라인을 구축하고자 하는 분들께 도움이 되었으면 하고 피드백도 환영합니다.

관련 포스팅 보기>

NGS 분석 파이프 라인의 이해: GATK Best Practice

[계속 Update 예정] 자주 쓰는 linux 명령어 및 프로그램 관련 자료

I. 들어가기에 앞서 (준비물)

우선 기본적으로 모든 작업은 리눅스 (linux) 환경에서 이루어집니다. 리눅스 환경에서 기본 명령어를 숙지하고 필요한 프로그램들을 설치합니다. 기본적으로 아래와 같은 프로그램들을 다운로드 및 설치합니다. 자세한 설치 방법의 경우, 해당 웹사이트에서 자세히 제공하고 있습니다. 저의 경우, 현재 최신버젼인 BWA 0.7.5, PICARD 2.20, GATK 4.1.3을 기준으로 설치 및 작업을 진행하였습니다.

BWA (Burrows-Wheeler Aligner): http://bio-bwa.sourceforge.net/

Samtools: http://www.htslib.org

PICARD: https://broadinstitute.github.io/picard/

GATK4: https://software.broadinstitute.org/gatk/download/

위의 프로그램이 설치가 끝나면, path 등록을 해서 리눅스 상에서 명령어를 입력하면 잘 실행되는지 확인합니다. 마지막으로 Reference Genome 파일을 다운로드 합니다. hg19 또는 hg38 버젼에 맞게 파일을 모두 다운로드하여 위치시키면 기본적인 준비들이 끝납니다. (파일은 UCSC genome browser에서 다운 가능)

II. Data pre-processing: Fastq to Bam

아래 코드는 GATK best practice의 Data pre-processing step에 해당하는 내용을 정리한 것입니다. 전체 Work-flow는 아래와 같습니다.

GATK Best practice: Data pre-processing for variant discovery

Fastq to Bam

1. Map to Reference: Fastq to Bam

bwa index [hg19_reference.fa]
bwa mem -t 8 -M -R '@RG\tID:foo\tSM:bar\tLB:library1' [hg19_reference.fa] [sample01_R1.fq.gz] [sample01_R2.fq.gz] > [sample01.sam]
samtools view -bt [hg19_reference.fa.fai] [sample01.sam] -o [sample01.bam]
samtools sort [sample01.bam] -o [sample01.sorted.bam]
samtools index [sample01.sorted.sam]

Reference sequence를 다운로드 받았으면, 기본적으로 indexing을 한번 진행해 준다. 이후, 위의 명령어는 Paired-end Read의 2가지 fastq 파일을 1개의 Sam 파일로 mapping하는데 필요한 명령어이다. []안의 부분을 적당히 변경하여 이용할 수 있다. @RG에 해당하는 부분은 기본적인 read에 붙은 바코드 정보를 말한다. 실제 명령어에서는 대괄호는 [] 입력하지 않는다. Fastq 파일은 압축 형태를 바로 이용할 수 있으므로, 압축을 모두 풀 필요는 없다. -t 는 thread option으로 컴퓨터 성능에 따라 적절히 지정해준다. samtools를 이용하여, sam > bam 변환 및 동시에 sorting, indexing을 진행한다.

2. MarkDuplicate + SortSam

java -jar picard.jar MarkDuplicates I=[sample01.sorted.bam] O=[sample01.dedup.bam] M=[sample01.markdup.metrics.txt]
java -jar picard.jar SortSam I=[sample01.dedup.bam] O=[sample01.sortsam.bam]

MarkDuplicate와 SortSam을 이용하여 중복되는 Read를 걸러주고, Read를 정렬시켜서 이후에 Variant call을 위한 기본 작업을 진행한다. 이 부분은 과거 Samtools를 이용하여, 주로 이용하였는데, SortSam으로 대체해서 진행이 가능하다.

3. Base Quality Score Recalibration (BQSR)

gatk --java-options 'Xmx16g' 'Xms8g' BaseRecalibrator -I [sample01.sortsam.bam] -R [hg19_reference.fa] --known-sites [dbsnp_151.GRCh37.vcf] -O [sample01.recal_data.table]
gatk --java-options 'Xmx16g' 'Xms8g' ApplyBQSR -I [sample01.sortsam.bam] -R [hg19_reference.fa] --bqsr-recal-file [sample01.recal_data.table] -O [sample01.final.bam] 

BQSR는 시퀀싱을 통해 생성된 염기들의 Quality Score에 발생하는 bias를 보정해서, 잘못된 mapping을 보정해주는 QC 과정의 일종이다. 위의 명령어를 성공적으로 실행시키면, 최종적으로 분석 가능한 Bam 파일이 생성된다. 생성된 Bam 파일은 이후에 변이를 Call하여 VCF 파일을 생성하는데 이용된다.

[References]

GATK Best practice: Data pre-processing for variant discovery

NCI GDC DNA-Seq Analysis pipeline: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline/

Samtools (Mapping to Reference): http://www.htslib.org/workflow/#mapping_to_variant

 

NGS 분석 파이프 라인의 이해: GATK Best Practice

최근들어 제 블로그의 방문자들이 눈에 띄게 증가했습니다. 대부분 구글 검색을 통해서 유입되는 분들인데, 검색어 통계를 보니 NGS 관련 내용에 대해서 검색하다가 들어오는 사람들이 많았습니다. 그래서 저에게는 쉬운 내용들이라 하더라도, 일반인들이 궁금해하는 내용에 대해서도 정리해두어야 겠다는 생각이 들어서 이번 포스팅은 NGS 장비를 통해서 생산된 데이터를 어떻게 분석하는지, 전반적인 파이프 라인에 대해서 개념을 소개하는 포스팅을 올려볼까 합니다.

우선, 차세대 염기 서열 분석법 Next-generation sequencing (이하 NGS)은 다양한 이름으로 불리는데, 이미 널리 사용하는 기술이기 때문에 현재는 차세대 기술이라고 보기 어렵습니다. 따라서, NGS라고 부르는 건 misnomer라고 할 수 있죠. 좀 더 정확한 명칭으로는 Massively parallel sequencing, High-throughput sequencing 등이 있는데, 동일 기술을 가르킨다고 생각하면 됩니다.

관련 포스팅 보기 >

휴먼 게놈 프로젝트, 그 이후.. 그리고 정밀 의료시대까지

NGS 검사: Whole Genome & Exome, Targeted Sequencing 비교

시퀀싱 기술의 발전으로 현재는 NGS를 널리 사용하게 되었는데, 다양한 방식이 있지만 지금은 Illumina 사의 Flow cell 기반의 short-read 시퀀싱 방식이 대부분의 시장을 독점하고 있는 상태입니다. 따라서, 아래의 NGS 분석 방법은 Illumina 시퀀싱 방식으로 생산되는 read를 기준으로 설명하도록 하겠습니다.

[Illumina sequencing 과정 소개 You-tube 영상]

 

시퀀싱 데이터 분석 과정은 위의 과정을 통해서 생산된 매우 많은 짧은 가닥의 read (50~150 bp 염기)들을 적절한 유전체의 위치에 퍼즐로 끼워 맞추고, 기존에 알려진 표준 유전체 (Reference sequence)와 대조하여, 바뀐 염기나 변이가 있는지를 검출하는 전체 과정을 말합니다. 따라서, 크게 아래와 같은 과정을 거쳐야, 생산된 read들로 부터 변이를 검출할 수가 있습니다.

GATK에서는 Germline 또는 Somatic variant 에 따라, 그리고 타겟 변이의 특성에 따라, 서로 다른 Best practice를 제공하고 있습니다. 이 중에서 아래는 가장 보편적인 Germline short variant 발굴 과정을 나타내고 있습니다.

[GATK Best Practice 보기]

Germline short variant discovery

Germline copy number variant discovery

Somatic short variant discovery

Somatic copy number variant discovery

 

gatk
[GATK best practice] GATK에서는 생산된 read로부터 변이를 검출하는 전체 과정에 대한 표준화 지침을 제공하고 있는데, 이를 GATK best practice라고 하며 자세한 과정은 GATK forum에서 step-by-step으로 제공해주고 있습니다.

I. 표준 유전체 서열에 read를 정렬하기 (FASTQ to SAM): 보통 NGS를 통해서 생산되는 개별 Read의 개수는 백만개 이상의 단위가 됩니다. 이때, 개별 read의 정보는 FASTQ 파일로 저장되고 그 크기는 수십~수백 Gb 단위가 됩니다. 각각의 생산된 read는 이미 알려진 인간 표준 유전체 서열에 가장 잘 맞는 위치에 정렬시키게 됩니다. 예로, 100 bp read를 기준으로 하여, 100개의 서열이 모두 일치하는 경우는 거의 유일하게 되므로, 해당 위치에 잘 찾아가게 됩니다. 다만, 반복 서열이나 특이적이지 않은 서열의 read는 제대로 mapping이 되지 않을 수가 있는데, 이는 Illumina 방식이 가진 한계점입니다. 보통 이 과정은 BWA (Burrow-Wheeler Aligner)라고 하는 프로그램으로 진행하게 됩니다.

fastq
개별 Read는 위의 그림과 같은 정보를 포함한 FASTQ 파일 형식으로 생산됩니다.

 

II. 정렬된 정보를 binary format으로 변경하기 (SAM to BAM): 위의 과정을 거친 SAM 파일은 여전히 용량이 매우 큽니다. 따라서, 용량을 절약하기 위해서 컴퓨터가 이해하는 2진법의 binary 형식으로 변경하면서 용량을 줄이고 연산 속도를 올립니다. 이 과정을 거친 파일이 흔히 말하는 1차적인 BAM 파일이 됩니다.

bam
인간 표준 유전체 서열에 개별 Read들이 정렬된 모습. BAM 파일

 

III. 개별 위치의 Base quality 보정하기 (Quality Control 과정): 1차적으로 read들이 표준 유전체에 정렬되면, 이것이 제대로 찾아 들어간 것인지를 평가하기 위해, 유전체 개별 위치에 대해 각각의 read 정보를 토대로 제대로 정렬된 것인지를 평가하고, 보정해주는 과정을 거칩니다. 특히, Indel이 발생한 read들의 경우, bias가 크기 때문에 따로 Indel realignment 라는 과정을 거치고, 개별 염기 위치에 대해서도 다시 한번 보정을 해주는 Base Recalibration 과정이 존재하게 됩니다.

IV. Variant Calling (BAM to VCF): 마지막 과정은 BAM 파일에서 실제로 변이를 검출하여, 변이들만 추출하는 과정입니다. Germline인지 Somatic mutation 인지에 따라서 다양한 알고리즘을 이용하게 되는데, Germline의 경우 가장 대표적으로 Haplotypecaller를 이용하게 됩니다.

V. Variant Annotation: 4번 단계까지 거치면, 무수히 많은 변이 정보를 포함하는 VCF 파일이 얻어지게 됩니다. 이렇게 call된 변이들 중 일부는 error를 일부 포함하고 있기 때문에, 보정 및 QC 과정을 한번더 거치게 됩니다. (Variant Recalibrator) 이후에 QC 과정을 거쳐서 Filtering된 변이들 중에서 관심 있는 변이만 얻기 위해서는, 기존 데이터 베이스의 자료를 토대로 각각의 검출 변이에 대해서 신원을 식별하고, ID를 발급하는 일종의 annotation 하는 과정이 필요합니다. 보통 이 과정에서 다양한 툴들이 이용되는데, 가장 대표적으로 Annovar 프로그램을 이용하여, gnomAD DB 등의 자료를 이용하게 됩니다. 최근에는 GATK의 Funcotator에서 기본적인 annotation을 지원하고 있습니다.

 

개별 과정의 코드까지 전부 올리기에는 너무 양이 많아질 것 같아서, 이번 포스팅에서는 전체 흐름에 대해서만 간략히 다루도록 하겠습니다. 다음 포스팅에서는 개별 변이를 해석하는 방법에 대해서 조금 더 자세히 다루도록 하고, 이번 포스팅은 여기서 마치도록 하겠습니다.

 

[Reference]

GATK Best Practice Forum https://software.broadinstitute.org/gatk/

Annovar http://annovar.openbioinformatics.org/en/latest/

NGS 데이터를 이용한 CNV 분석

Copy Number Variation (CNV)는 Single Nucleotide Variation (SNV)과 더불어, 유전적 다양성을 나타내는 주요한 원인으로 생각되고 있습니다. 유전자 sequence의 염기 하나가 치환된 SNV와 달리 CNV는 유전자 전체 또는 일부의 copy가 중복되거나 줄어들수도 있기 때문에 훨씬 넓은 영역에서 나타나는 Structural variation의 일종입니다.

관련 포스팅 보기 -> 유전학 중요개념 정리: Structural variation 및 Copy-number variation

일반적으로 NGS는 SNV를 보기 위한 목적으로 검사를 시행하지만, 해당 데이터를 활용하면 CNV 분석도 할 수 있기 때문에, 오늘은 NGS 데이터를 활용한 CNV 분석 방법에 대해 포스팅하고자 합니다.

NGS CNV
[그림1. NGS 데이터를 이용하여 CNV를 검출하는 원리] CNV 검출을 위해서는 mapping 되는 read 간의 정보, 그리고 각 영역에 mapping된 read의 depth 정보를 활용하게 됩니다.

위의 그림은 CNV 분석을 위한 NGS 데이터의 5가지 활용 원리를 나타내주고 있습니다. 그러나 가장 핵심이 되는 원리는 Read depth입니다. Target sequencing과 같이 Read depth가 충분한 경우에, 다른 검체들에 비해 해당 영역의 depth가 월등히 떨어지거나, 또는 월등히 높은 경우에는 해당 영역의 deletion 또는 duplication을 의심할 수 있습니다.

target_coverage_nd_FGFR2_4
[그림2. FGFR2 유전자의 Coverage (위) 및 Reference의 depth로 normalized한 depth (아래)를 나타내는 도표] 다른 검체들보다 Normalized depth가 월등히 높은 검체 (P27)는 해당 영역의 duplication, 월등히 낮은 검체 (P33)는 해당 영역의 deletion이 존재하는 것으로 의심할 수 있다.

사실 NGS 데이터는 CNV를 목적으로 한 것이 아니라, SNV 검출 목적의 데이터를 부수적으로 활용하는 것이기 때문에 많은 제한점이 있습니다. 따라서, 임상적으로 CNV 검사 목적의 NGS는 권장되지 않으며 적절한 가이드라인도 존재하지 않기 때문에 다양한 Computational tool 들이 개발되어 서로의 장점을 홍보하는 상황입니다. 다음은 다양하게 개발된 대표적인 CNV 검출 tool 들을 정리한 표입니다. 많은 경우  BAM 파일을 활용하는 것을 볼 수 있으며 대부분 R package를 제공하고 있어, 사용이 용이합니다.

NGS CNV2
[그림 3. CNV 검출을 위한 다양한 컴퓨터 툴] 어떠한 툴이 우수한가에 대해서는 명확하게 정립된 결론이 없기 때문에, 적절한 상황에 맞게 툴들을 활용하는 것이 필요합니다.
위의 표와 같이 다양한 툴들이 존재하지만, 실제로 몇가지 툴들을 사용하여 봤을 때, 결과들이 제각각이었고, 서로 일치하는 정도도 높지 않았습니다.  다양한 알고리즘을 활용함에도 불구하고, 위양성으로 보고되어 믿기 어려운 경우가 많았습니다. 가장 정확한 방법은 직접 그림 2와 같이 해당 영역의 coverage plot과 normalized depth를 보고 종합적으로 판단하는 것이었습니다. 아직까지 컴퓨터 툴들에 개선의 여지가 많음에도 불구하고, NGS 데이터를 활용하면 CNV에 대한 정보도 일부 얻을 수 있기 때문에 NGS는 더 폭넓게 활용될 것으로 전망이 됩니다.

[Reference]

Zhao, Min, et al. “Computational tools for copy number variation (CNV) detection using next-generation sequencing data: features and perspectives.” BMC bioinformatics 14.11 (2013): S1.