[NGS DNA-SEQ] Functional Equivalence pipeline: CROMWELL, WDL

gnomAD, TOPMed 등 대규모 유전체 코호트들이 만들어지면서, 여기서 생산된 데이터를 이용하는데 중요한 문제가 부각되었는데, 바로 분석 결과 간의 재현성호환성에 있었습니다. 즉, 연구자가 GATK Best practice를 이용하여 Exome 또는 Genome 시퀀싱 분석을 진행하더라도 어떤 설정과 파라미터 값을 넣느냐에 따라, 최종 검출 변이의 결과가 달라지고, 이것은 연구 결과 간의 재현성의 측면에서 매우 중요한 문제가 되었던 것이지요. 참고논문에서 진행한 실험 결과를 보면, 동일한 샘플로 생산된 FASTQ 파일을 서로 다른 5개의 기관에 보내 각각의 파이프 라인으로 분석한 결과, Call된 변이들 간에 많은 차이가 있었다고 보고하고 있습니다.

[관련 포스팅 보기]

따라서, 점점 늘어나는 유전체 데이터만큼 유전체 분석 파이프 라인을 하나의 표준화된 파이프 라인으로 통합하는 것이 매우 중요해졌고, 그 결과 개발된 것이 “Functional Equivalence (FE)” Pipeline (기능적으로 동등한 파이프라인) 입니다. 사실 NGS 분석을 하는 사용자의 입장에서는 Input만 넣고, Output만 나오면 편한데, 그동안 개발된 툴들은 이를 모두 아우르는 것이 아니라, 그때 그때마다 필요한 부분들을 개발했기 때문에, 분석 파이프 라인도 이제야 어느 정도 성숙 단계에 이르렀다고 할 수 있습니다. 따라서 최근의 대규모 유전체 컨소시엄들은 모두 “Functional Equivalence ” Pipeline 을 통해 생산된 유전체 데이터를 생산하는 것으로 채택하고 있습니다. (그래서 저도 functional equivalent한 결과를 얻기 위해서 최근에 새롭게 공부를 하게 되었습니다.)

[ Functional Equivalence Pipeline Overview]

FE 파이프라인을 제공하기 위해서, Broad Institute의 개발진들은 WDL (Workflow Description Langauge)과 Cromwell이라고 하는 프로그래밍 언어를 개발하는데, 하나의 파이프라인을 패키지로 묶은 WDL 파일을 만들고, 이를 Cromwell이라는 프로그램으로 구동시켜주는 원리라고 합니다. 사실 사용자의 입장에서는 과거에 개별 프로그램을 설치하고, 개별 스텝을 따로 돌려야했다면, 지금은 이러한 것들이 모두 하나의 패키지 형태로 제공되기에 더욱 편해졌다고 할 수 있습니다 (?).

[Cromwell 페이지 바로 가기] https://cromwell.readthedocs.io/en/stable/

[WARP 페이지 바로가기] https://broadinstitute.github.io/warp/docs/get-started/

Cromwell의 로고: 꼬마돼지 베이브와 스타트렉에 출연한 배우 James Cromwell을 오마주한 로고라고 합니다 🙂

병원 검사실을 운영하는데, 검사 장비와 보고 방법을 표준화하는 것은 매우 중요한데, 유전체 분석 파이프 라인에도 이제야 이러한 개념이 들어왔다는 점이 반갑습니다. (분석 파이프 라인 하나도 이렇게 표준화하기가 어렵습니다.) Genome의 경우에는 처리해야할 데이터의 크기가 워낙 방대하기때문에 더욱 어려운 점이 있는 것 같습니다. 현재 이쪽 분야도 많은 Computational Scienctist들이 뛰어들어서 개발을 진행하고 있는 중이기 때문에, 추후에 더 User-friendly하고 간편한 파이프 라인이 개발되어 제공되기를 기대해 봅니다. (점차 대세는 클라우드로 옮겨가지 않을까 합니다?)

[References]

Regier, Allison A., et al. “Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects.” Nature communications 9.1 (2018): 1-8.

NGS DNA-seq pipeline: GATK Best Practice Code – Part3. Vcf manipulation

앞선 포스팅의 두가지 과정을 거치고 나서 생성된 VCF 파일을 이용하면 드디어 분석 가능한 변이들을 확인할 수 있습니다. 그러나 실제로 이 데이터를 열어보면, 지저분하고 활용하기 위해서는 어느 정도 가공이 필요합니다. 그래서 이번 포스팅은 VCF Filter를 적용하여 분석을 위한 변이들을 정제하고, 분석에 참조하기 위한 Annotation 작업을 위한 Code까지 정리해보도록 하겠습니다.

관련 포스팅 보기>

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

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

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

vcf_3

[일반적인 VCF 파일의 구조]

IV. Variant Filtration and Annotation : Genotype Refinement Workflow

VCF call을 하고 나면, 많은 변이 데이터가 생산되는데, 실제로 데이터의 퀄리티와 오류를 보정하는 작업이 필요합니다. 이 부분은 아직까지 확립된 최적의 Work Flow는 없지만, GATK에서는 Genotype Refinement Workflow를 제공하고 있습니다. 여기서는 Call된 genotype이 인구집단이나 가계도 정보에 근거하여 적절한 Call인지를 판단하게 됩니다. 가계도의 경우에는 양쪽 부모 정보가 모두 (Trio) 필요합니다. 아래에서는 이 Workflow를 따라 Code를 작성해보도록 하겠습니다.

Genotype refinement

1. CalculateGenotypePosteriors

VCF call을 통해 생성된 genotype을 기반으로 이번에는 거꾸로 Genotype의 posterior를 계산해줍니다. 이 결과를 토대로 변이 call을 판정해서, Genotype의 quality를 계산합니다. 이 과정에는 인구 집단이나 가계도의 Prior 정보를 이용하는데, 가계도 정보를 이용하기 위해서는 PED 파일 생성이 필요합니다.

PED 파일 생성 정보는 아래 링크를 참고하시기 바랍니다.

http://zzz.bwh.harvard.edu/plink/data.shtml

gatk --java-options 'Xmx16g' 'Xms8g' CalculateGenotypePosteriors -V [cohort.recal.vcf] \
-ped family.ped --skip-population-priors -O [cohort.GP.vcf]

2. Filter Low Quality Genotypes

위의 과정에서 계산된 GQ 값을 기준으로 Filter 조건을 걸어서 지저분한 변이 call들을 제외 시킬 수 있습니다. 아직까지 표준화된 Filter 조건은 없기 때문에, 사용자가 적절하게 필터 조건을 설정해주어야 합니다. 여기서 선언해주는 조건은 흔히 JEXL expression이라고 부르는 방법으로 작성을 해야한다고 하는데, 여기서는 널리 쓰이는 Filter option인 GQ < 20 을 걸도록 하겠습니다.  필요에 따라 수치를 변경시키거나, 조건을 바꾸면 됩니다. 자세한 내용은 아래에서 찾아보시기 바랍니다.

article about using JEXL expressions

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

마지막으로 filter조건에 부합하는 변이만 남기기 위해, GATK의 SelectVariants 이용하도록 하겠습니다.

gatk --java-options 'Xmx16g' 'Xms8g' VariantFiltration -R [hg19_reference.fa] -V [cohort.GP.vcf] \
--genotype-filter-expression "GQ<20" --genotype-filter-name "lowGQ" -O [cohort.filter.vcf]
gatk --java-options 'Xmx16g' 'Xms8g' SelectVariants -V [cohort.filter.vcf] --exclude-filtered true --max-nocoll-fraction 0.1 -O [cohort.filtered.vcf]

3. VCF annotation using Annovar

이제 어느 정도 정제가 된 변이 데이터셋이 생성되었습니다. 이제 불러온 변이들에 맞는 annotation 작업을 진행하면 분석을 진행할 수가 있습니다. 최근에 GATK에서 Funcotator를 개발하여, annotation이 가능해졌는데 아직까지는 Annovar를 이용하는게 대세이므로, Annotation 작업은 Annovar를 이용하도록 하겠습니다.

관련 포스팅 보기>

Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용

Annovar 홈페이지를 통해 다운로드 신청을 하면,  등록된 메일로 Kai Wang에게서 자동으로 다운로드 url을 받을 수 있습니다. Annotation을 원하는 Database의 크기가 크기 때문에, 처음에는 DB 다운로드에 시간이 많이 소요가 될 수 있습니다.

Annovar를 설치하고, 아래 부위에 원하는 DB의 이름을 찾아서 넣으면 다운로드가 진행됩니다.

annotate_variation.pl -buildver hg19 -downdb -webfrom annovar [refGene] humandb

Annotation은 아래 code를 통해서 진행합니다. 여기서는 refGene, cytoBand, avsnp150, dbnsfp35a, exac03, gnomad_exome, clinvar_20190305 DB를 이용하여 annotation 작업을 진행하도록 하겠습니다. 각각의 DB는 Gene-based, Region-based, Filter-based annotation에 해당하는 카테고리가 있고, 이를 operation 옵션에 g, r, f를 이용하여 나타내줍니다.

table_annovar.pl [cohort.filtered.vcf] humandb/ -buildver hg19 -out [myanno] -protocol refGene,cytoBand,avsnp150,dbnsfp35a,exac03,gnomad_exome,clinvar_20190305 -operation g,r,f,f,f,f,f -nastring . --vcfinput --remove --thread 8

[References]

VCFtools를 이용한 VCF filtering: 위에서는 GATK를 이용하였지만, 과거에 많이 사용하였던 VCFtools를 이용하는 방법도 있습니다. 자세한 내용은 아래를 참조하시기 바랍니다.

http://www.ddocent.com/filtering/

GriffithLab에서 소개하고 있는 VCF filter 과정: 해당 내용은 GATK4로 넘어오면서 코드가 약간 달라졌습니다. 다만 전반적인 과정은 유사하므로, 참고하시기 바랍니다.

https://pmbio.org/module-04-germline/0004/02/02/Germline_SnvIndel_FilteringAnnotationReview/

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&nbsp;-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&nbsp;-V [vcf file_list] \ [-L targets.interval_list] --genomicsdb-workspace-path \ [/genomicDB] --merge-input-intervals true --tmp-dir=/tmp&nbsp;

[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

 

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

[Linux 자주 쓰는 명령어]

  • ls: 현재 위치 display
  • cd /folder : 이동
  • rm -rf /folder : 하위 폴더 포함 전체 삭제
  • mkdir /abc: abc 디렉토리 생성
  • chmod -R 777 ‘file1’: file 1 권한 변경 / 읽기: 4, 쓰기: 2, 실행: 1
  • cat, grep: cat /BiO/project/example.txt | grep “test”  내용 중 특정 문자열을 포함하는 줄만을 불러옴.
  • screen -ls : screen list 표시 / -S ‘screen1’ : screen1 생성 / -r ‘screen1’ : screen1 재입장 / Ctrl + A + D: screen 나가기
  • pkill screen : screen 모두 강제 종료
  • gzip ‘file1’ : file1 압축하기 / gzip -d ‘file1’ : file1 압축 풀기
  • tar 압축하기 : tar -cvf [파일명.tar] [폴더명]
  • tar 압축풀기: tar -xvf [파일명.tar]
  • tar.gz 압축하기 및 풀기: tar -zcvf / tar -zxvf : option에 z 추가
  • 여러개 gz 파일 합치기: cat file1.gz file2.gz file3.gz > allfiles.gz
  • cat /proc/cpuinfo | grep name: 서버 CPU 모델명 확인
  • nproc: 서버 코어 갯수 출력
  • df -h: 시스템 디스크 공간 확인
  • free -h: 현재 메모리 상황 확인
  • sar -r 10: 10초 간격으로 메모리 사용 현황 표시
  • rsync [file1] [path] -r –progress : 하위 폴더 포함 복사 & 복사 진행 상황 함께 보기
  • sudo: super user do – root 권한으로 실행
  • apt-get upgrade / apt-get update / apt-get install : 우분투 관련 패키지 업데이트 및 설치
  • 파일명 변경: mv or rename 이용
  • mv [파일명1] [파일명2]
  • rename [변경전문자] [변경후문자] [대상파일]: rename은 파일패턴문자(?, * 등)를 인식하기 때문에 여러 개의 파일에 동시 사용이 가능. Regex 정규 표현식을 이용해 사용 (링크 참조)

ex> file1~file9.png > file1~file9.jpg 확장자 변경
>rename png jpg file?.png

https://qgp9.github.io/blog/2016/03/24/rename-batch-rename-files

 

[Linux PATH 관련 설정]

  • echo $PATH: 현재 걸려있는 PATH 모두 확인
  • export PATH=”$PATH:/path/to/dir”: 경로 디렉토리를 일시적으로 추가

경로 디렉토리를 영구적으로 추가하기 위해서는, vi 에디터를 이용하여 ~/.profile에 위 경로를 추가하여 수정해주어야한다. (우분투 16의 경우, /etc/environment 에 등록하여도 됨.) 저장하고 나가기”:wq”
>sudo vi /etc/environment
export PATH="$PATH:/path/to/dir

이후에 해당 내용을 적용시키려면, source code를 실행시킨다.
>source /etc/environment

[java 관련 설정]

GATK 4를 돌리기 위해서는 JAVA 1.8 이상의 버젼이 필요한데, 최근에 오라클에서 유료로 전환되면서 업데이트가 힘들어졌다. 다행히 GATK는 무료로 배포되고 있는 Open JDK를 지원하기 때문에, Open JDK로 업데이트를 할 경우에 무리 없이 실행이 가능하다. 다음은 Ubuntu에서 Open JDK의 하나인 Zulu-8으로 업데이트를 위한 실행 명령어이다.

>sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 0x219BD9C9
>sudo apt-add-repository 'deb http://repos.azulsystems.com/ubuntu stable main'
>apt-get update
>apt-get install zulu-8

위와 같이 zulu를 설치 후에, 기본 자바 홈을 vi 에디터를 이용하여, 설치 폴더로 변경해주어야 한다.

> vi /etc/profile 실행, 기존 자바 홈 설정 문구를 지우고, 아래 문구를 추가해준다. “:wq”
export JAVA_HOME=/usr/lib/jvm/zulu-8-amd64

>apt-get remove oracle*
>apt-get autoremove --purge
>apt-get autoclean

마지막으로 Java version을 확인하여, 설치가 제대로 되었는지 본다.
>> java -version
openjdk version "1.8.0_181"
OpenJDK Runtime Environment (Zulu 8.31.0.1-linux64) (build 1.8.0_181-b02)
OpenJDK 64-Bit Server VM (Zulu 8.31.0.1-linux64) (build 25.181-b02, mixed mode)

[Java option에 따른 의미]

>C:\java –X
–Xms : set initial Java heap size
–Xmx : set maximum Java heap size
–Xss : set java thread stack size

[conda, pip SSL error 관련]

> conda config --set ssl_verify False
> pip --trusted-host pypi.org --trusted-host files.pythonhosted.org install pip

회사 네트워크 보안 등의 설정으로 방화벽이 걸려있는 경우, SSL 인증서 문제가 발생하면서 업데이트가 안된다. 통신 에러 또는 방화벽 차단에 의해 conda update가 안되는 경우에는,  위의 명령어를 이용, ssl_verify 옵션을 꺼준다. pip의 경우는 –trusted-host 부분에 다운 받을 저장소의 주소를 등록하여 준다. 매번 이런 명령어를 쓰는 것이 귀찮으면, config 파일에 아래와 같이 등록하면 된다.

>sudo vi /etc/pip.conf
[global]
trusted-host = pypi.org
files.pythonhosted.org

[Samtools 이용 bam file header 변경]

for file in *.bam; do filename=`echo $file | cut -d "." -f 1`; samtools view -H $file | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - $file > ${filename}_chr.bam; done

Last updated 2020-08-13

Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용

오늘 포스팅은 Annovar에 관한 내용으로 작성해볼까 합니다. 언젠가 한번은 정리할 필요가 생각하고 있던 내용인데, 이제서야 차례가 왔네요. 시퀀싱 기술이 보편화되면서, 대부분 표준화된 파이프라인을 이용하여 시퀀싱 raw data로부터 변이를 검출하는 Up-stream analysis는 대부분의 업체에서 대동소이한 결과를  주고 있습니다. 연구자의 입장에서 사실 더 중요한 것은 이 정보를 어떻게 효과적으로 이용할까 하는 부분에 있고, 그래서 Bam file 또는 VCF file에서 시작하는 Down-stream analysis가 더 중요하다고 볼 수 있습니다. 그 첫 단계로 필요한 적절한 정보를 활용하여 주석을 다는 과정이 있고, 이를 우리는 Annotation (주석 달기) 과정이라고 합니다. 그리고 이 과정에서 우리는 대부분 Annovar를 사용하게 됩니다. 물론 업체에 의뢰하면 대부분 기본적인 Annotation이 끝난 파일도 전달을 받게 되는데, 사실 필요 없는 내용이 잔뜩 달려서 파일의 용량만 무지막지하게 커진다거나, 정작 필요한 내용이 빠진 경우도 종종 발생합니다. 그래서 이번에 다룰 내용은 주석 달기의 각 항목에 대한 의미와 주요 활용 항목에 대해서 정리해보겠습니다.

관련 포스팅 보기>

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

Annovar는 Perl script 기반으로 짜여 있어서, linux에서 명령어를 주면 자동으로 annotation을 달아주게 되어있습니다. 이때 몇가지 option에 따라서 원하는 내용대로 다양한 Database에서 관련 내용을 annotation 할 수가 있습니다. 더 많은 내용을 추가함에 따라서 연산 시간과 최종 파일 용량은 커지고, 가끔은 이것이 오히려 독이 되기도 합니다. 따라서 원하는 내용을 입맛에 따라 최적화하여 활용하는 것이 필요하게 됩니다. 아래 그림은 Annovar를 이용하여, 주석을 다는 과정의 전반적인 흐름을 보여주고 있습니다. 간단히, “Database 파일을 통째로 다운 받아서, Perl script 기반의 annotate_variation.pl을 실행하면, 그에 맞는 내용을 찾아서 주석으로 달아준다.” 정도로 말할 수 있겠습니다.

ANNOVAR_main_package_workflows.svg.png

Annovar에서 제공하는 Database 목록 보기

위 페이지에서는 다운로드 가능한 Database의 목록들을 보여주고 있습니다. 필요한 항목의 이름과 업데이트 날짜 등을 참고하여, DB를 다운로드하고 적절하게 활용하기 바랍니다.

annotate_variation.pl -buildver hg19 -downdb -webfrom annovar [위 목록의 Database 이름] humandb/

 


I. Gene-based Annotation

VCF 파일에서 가장 핵심 정보는 사실 몇개 없습니다. 정확하게는 5개 정보만 있어도 충분한데, “몇번째 염색체 (Chr)의 몇번째 염기 위치 (Position)가 원래 무엇인데 (Ref) 무엇 (Alt)으로 바뀌었다.” (흔히 이러한 형태의 파일을 MAF format이라고 합니다.) 이 말을 하기 위해서, 부수적인 내용들이 잔뜩 달려 있는 셈이죠. 그러나 이 정보만을 이용해서는 사람들이 알아들을 수가 없습니다. 우리는 “어떤 유전자의 몇번째 엑손 영역의 몇번째 아미노산이 무엇으로 바뀌었다“와 같은 정보가 필요하죠. 그런데 사실 아직도 유전자의 정확한 정의와 위치에 대해서는 100% 밝혀지지 않았습니다. 따라서 이 유전자라고 하는 부분도 다양한 데이터 베이스가 존재하게 됩니다. 대표적으로 RefGene, UCSC/Ensemble Gene, Known Gene, CCDS 등등의 데이터 베이스가 있습니다. 따라서 여기서 어떤 데이터 베이스를 이용하여, annotation을 하냐에 따라서 출력이 달라집니다. (하지만 사실 큰 차이는 없습니다. 대부분의 업체에서는 2~3개의 데이터 베이스를 이용하여 annotation을 해주는데, 대부분 내용이 중복되어 용량만 커짐… 그래서 개인적으로는 그냥 RefGene만으로도 충분합니다.)

<가장 핵심적인 VCF 파일의 정보> = MAF format

Chromosome : Position (Start_End) : Reference sequence > Alternative sequence

위의 database를 이용하면, 위의 정보가 어떤 유전자에 속하고, 해당 유전자에서 어떤 기능을 하는 어떤 부위의 변이인지, 기능적으로 변화가 있는지 없는지 등에 대한 기본적인 정보를 제공해주게 됩니다. 어떻게 보면 가장 핵심적인 정보를 추가하는 부분이라고 할 수 있습니다.

 

II. Filter-based Annotation

사실 변이를 Genome Browser에서 찾는 가장 빠른 방법은 rsID를 이용하는 것입니다. 그런 점에서 dbSNP 또는 avSNP의 rsID를 주석으로 달아놓는 것은 활용도가 높습니다. rsID는 변이 보고가 점점 늘어남에 따라서 계속 갱신되고 있는데, 가장 최근 database는 avSNP 151 버젼이지만, 보편적으로 아직까지는 avSNP 147 버젼을 사용하고 있는 것 같습니다.

관련 포스팅 보기>

바이오 연구자를 위한 Genome Browser 비교 및 활용

아미노산 치환의 효과 예측: In silico tool의 원리와 종류

변이빈도와 효과 크기

사실 변이의 생물학적 의미를 해석하는 과정에서 정말로 중요한 내용은 일반 인구 집단에서 얼마나 흔하게 존재하는가?에 있다고 할 수 있습니다. 그러한 의미에서 Population Frequency 정보는 변이의 의미를 파악하는데 매우 중요한 단서를 제공해 줍니다. 이와 같은 Population Frequency는 다양한 집단에서 시퀀싱을 시행하여, 그 빈도를 계산한 다양한 데이터 베이스에 기반하게 되는데, 대표적으로 EXAC, gnomAD, 1000Genome, ESP database 등이 있습니다. 이때 중요한 점은 인구 집단이 얼마나 대표성을 띄는가?에 있다고 볼 수 있는데, 빈도를 계산한 인구 집단의 크기 (n수)와 인종적 배경 (Ethnicity) 등이 특히 중요합니다. 더불어, Rare variant의 경우는 더욱더 인종에 따른 차이가 크기 때문에 일반적으로는 한국인의 경우는 EAS (East Asian population) 정보를 이용하지만, 해당 데이터 베이스가 충분히 한국인을 대표하지 못하는 경우가 많습니다. (위의 데이터 베이스에서는 gnomAD가 가장 n 수가 가장 크기 때문에 저는 주로 EXAC과 gnomAD의 EAS 인구 집단 정보를 활용하고 있습니다.)

사실 데이터 용량이 가장 뻥튀기 되는 부분이 in-silico prediction tool 부분에 있습니다. 아미노산 치환 효과를 예측해주는 tool은 100가지가 넘게 있는데, 각 tool 들이 예측해주는 정보를 얼마나 포함할 것인가에 따라서 추가되는 정보도 달라집니다. 크게 얼마나 진화적으로 보존된 지역인지 또는 아미노산 구조적으로 변화를 유발하는지 등에 기반한 알고리즘으로 개발되어, SIFT, PolyPhen과 같은 고전적 tool부터, GERP, SiPhy, MutationTaster, FATHMM, MetaSVM, CADD, DANN과 같은 다양한 tool이 존재합니다. (역시나 이쪽도 절대 지존은 없기 때문에 경우에 따라 적절하게 활용합니다. 보통 Voting Method, 즉 10개 또는 선택된 갯수의 서로 다른 알고리즘 중에서 몇개가 Deleterious로 예측하는가? 와 같은 방식으로 많은 경우 활용하게 됩니다. 저는 개인적으로 Ensemble 기반의 tool을 주로 활용하고 있습니다.)

마지막으로 임상 정보들을 annotation 하는 database가 존재하는데, 대표적으로 OMIM, HGMD, ClinVar, COSMIC 등등이 있습니다. 희귀 유전 질환에 대해서 연구를 한다면, OMIM이나 HGMD, Cancer 관련 연구를 한다면, COSMIC database 정보를 annotation 하는 것이 도움이 될 수 있습니다. 그러나 사실 이렇게 annotation을 덕지덕지 붙이다보면 파일 크기가 엄청나게 불어나게 됩니다. 현재까지 대부분의 연구자들은 Coding region의 Functional variant에만 관심이 있기 때문에, 1차적으로 Gene-based annotation 후 exon 영역의 functional variant만 filtering하고나서, 해당 변이들에 대해서 annotation 하는 방법이 시간과 데이터를 절약하는 방법이 될 수 있습니다.

 

III. 기타 annotation 방법

관련 포스팅 보기>

암유전체 분석: Driver mutation prediction tools

위의 annovar를 이용하는 방법은 linux 기반의 서버를 통해서 대용량으로 실행하는 방법입니다. 그러나 서버를 구축하지 못하거나, linux를 친숙하게 이용하지 못하는 경우에는 그러면 어떻게 annotation을 하는가? 에 대한 문제가 발생합니다. 이를 위해서 다양한 Web 기반의 annotation tool 들이 존재하게 됩니다. 가장 대표적인 것이 wANNOVAR입니다. 기타 cancer를 다룬다면, Oncotator 또는 Cancer Genome Interpreter도 대안이 될 수 있습니다. 그러나 역시 이러한 tool들은 Annovar에 비해서 자유도는 떨어지기 때문에 기능에 제약이 있다는 단점이 있습니다. 마지막으로 R을 활용하여, annotation이 가능한 몇가지 package들이 개발되어 있습니다. 대표적인 package로는 MAFtools, VariantAnnotation 등이 있으나, 역시 기능이 AnnoVar에 비하면 제한적입니다. 그러나, 소수의 변이에 대해서 빠르게 annotation이 필요한 경우라면 이러한 도구들도 적절하게 활용하는게 도움이 될 수 있습니다. 더 자세한 정보는 아래 github를 활용하시기 바랍니다.

 


 

References>

Yang, Hui, and Kai Wang. “Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR.” Nature protocols 10.10 (2015): 1556.

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

wANNOVAR: https://github.com/WGLab/doc-ANNOVAR

Oncotator: https://github.com/broadinstitute/oncotator

VariantAnnotation: https://github.com/Bioconductor/VariantAnnotation

MAFtools: https://github.com/PoisonAlien/maftools

Web resources for Bioinformatics database: https://netbiolab.org/w/Web_Resources

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/