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

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

 

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에서 제공하는 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

Exome sequencing을 위해 고려할 요소들: capture kit와 target coverage 선택

Exome (엑솜)이란 유전자의 exon 영역을 모두 포함하는 집합체를 말합니다. 최근 가격적으로 시퀀싱 비용이 많이 떨어지면서, 연구 목적의 엑솜 시퀀싱이 광범위하게 사용되고 있는 것 같습니다. 사실 시퀀싱 이후의 데이터 분석도 중요한 과정이긴 하지만, 많은 사람들이 간과하고 있는 것이 그보다 상위 단계에 있는 시퀀싱 데이터 생산 계획 과정입니다. 여기서 말하는 시퀀싱 데이터 생산 계획이라 함은, 목적에 맞는 적절한 시퀀싱 장비 선정, 엑솜 키트 선택, 적절한 target coverage depth 선정 등을 모두 포함합니다.

Whole exome 은 일반적으로 모든 exon 영역을 지칭하지만, 사실 아직까지도 전체 유전자의 정체(?)를 정확히 모르고 있는 시점에서 논란이 많이 남아있는 부분이 있습니다. 일반적으로 시퀀싱 장비는 Illumina 사의 장비가 가장 보편적으로 이용되고 있기 때문에 차치하고라도, 유전체 (genome)에서 Exome 부분만 capture 하는데 사용되는 키트도 아래와 같이 다양한 제품이 존재합니다. 아래의 표에서 가장 눈여겨 볼 부분은 Target Region의 크기인데, 일반적으로 Exome이라고 말하는 부분의 크기도 39 ~ 64 Mb로 차이가 나는 것을 볼 수 있습니다. 이는 여러가지 기술적인 이유로 타겟 영역을 서로 다르게 디자인한 부분과 엑손 영역의 타겟 유전자의 수도 차이가 있기 때문입니다.

관련 포스팅 보기>

NGS Target enrichment method: Hybridization vs. Amplicon capture

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

시장의 선두 주자는 A사 였으나, 최근에 많은 회사에서 경쟁적으로 capture 효율과 coverage 를 개선시킨 저렴한 가격의 kit를 개발하여 공급하고 있습니다. I사의 경우에는 UK Biobank의 Exome sequencing에 이용되어, 호환성에서 장점이 있습니다. 사실 서로의 제품이 더 좋다고 광고하는 상황에서 Exome capture kit 선택은 쉽지 않지만, 최소 DNA 요구량, 관심이 있는 주요 유전자에서의 Coverage 정도, 그리고 가격 등의 요소를 종합적으로 고려하여 최적의 키트를 선택하는 것이 좋습니다. 제품마다 유전자별로 cover 되는 효율에 차이가 있는데, 특히나 관심이 있는 유전자들이 잘 cover 되지 않는 제품이라면, 검체에서 해당 변이를 검출하는 민감도에 큰 차이를 보일 수 있습니다. 특히나 낮은 수준으로 존재하는 종양의 변이들을 검출하는 연구 목적의 검사에서는 변이 검출 유무의 중요한 요소로 작용할 수도 있기 때문입니다.

일반적으로 시퀀싱 비용은 생산되는 데이터의 크기에 비례하여 증가하게 됩니다. 이때, 타겟 영역의 크기와 캡쳐 키트의 효율, 그리고 원하는 Coverage depth를 이용하면, 대략적으로 필요한 시퀀싱 데이터의 크기를 계산할 수 있습니다. 예를 들어, 위 표의 Agilent SureSelect V6 Exome kit를 이용하여, 100×의 depth로 시퀀싱을 하고 싶다고 가정하고, 일반적인 target capture 효율 (0.6 ~ 0.7)을 적용해봅시다. 아래 계산에 의해 총 10Gb의 데이터를 생산해야 원하는 coverage를 달성함을 계산할 수 있습니다.

(시퀀싱 데이터 크기) = (타겟 영역의 크기) × (Depth) / (On-target ratio)

(시퀀싱 데이터 크기) = 60 Mb × 100 / 0.6 = 10 Gb

 

시퀀싱을 위한 총 비용은 위 표의 capture kit 가격과 생산되는 데이터의 크기, 그리고 이용되는 시퀀싱 장비 및 검체 처리에 사용되는 시약의 가격 등이 합쳐져서 결정됩니다. 이 중에서 이용자가 결정할 수 있는 부분은 capture kit의 종류와 전체 시퀀싱 데이터의 크기 이기 때문에, 위의 내용들을 잘 숙지하여 필요한 만큼의 데이터를 효율적으로 생산하는 것이 중요합니다. 위의 내용은 Genohub의 Whole Exome Sequencing Guide의 내용을 정리한 것입니다. 마지막으로, 위에서 언급한 내용들을 정리하면서, 포스팅을 마치도록 하겠습니다.

Considerations for Whole Exome Sequencing

1. What sequencing instrument and read length should I choose for exome-seq?

2. How much sequencing coverage do I need for exome sequencing?

3. How do I calculate the sequencing coverage or depth required for my whole exome sequencing study?

4. Which exome sequencing capture kit should I use for my study?

5. How can I compare the annotation and exome capture design between each kit?

 

[Reference]

Genohub: Whole Exome Sequencing Guide

Goh, Gerald, and Murim Choi. “Application of whole exome sequencing to identify disease-causing variants in inherited human diseases.” Genomics & informatics 10.4 (2012): 214.

암유전체 분석: Driver mutation prediction tools

이전 포스팅에서 암에서 발생하는 mutation을 driver와 passenger로 구분하는 개념에 대해서 언급했습니다. 이번에는 실제로 시퀀싱을 진행했을 때 검출되는 많은 변이들을 실제 암 발생의 driver와 passenger를 구분하는 방법과 다양한 툴들에 대해서 정리해 보고자합니다.

관련 포스팅 보기>

[유전학 중요개념 정리] Driver vs. Passenger mutation in cancer

[유전학 중요개념 정리] Mutational signature

사실 Somatic mutation이나 Germline mutation이나 질병 발생의 원인 유전자와 변이를 찾는 방법이라는 데에서 큰 틀의 접근 방법은 동일합니다.  Germline 변이를 판독하는 큰 틀을 제시하는 가이드라인이 ACMG guideline이라고 한다면, Cancer 변이 판독의 기준으로는 흔히 AMP/ASCO/CAP guideline의 Tier system이 사용되고 있습니다. 즉, 개별 변이들을 아래와 같이 증거 수준과 임상적 중요도 등에 따라 Tier 1~4로 구분을 하여, 중요도가 높은 변이들을 임상적으로 활용하는 것이지요. 하지만, 이 역시도 한계가 많고 구분도 모호하기 때문에, 실질적으로 검출된 변이들의 driver mutation을 예측할 수 있는 다양한 툴들이 개발되고 있습니다.

[Somatic variant 변이 판독을 위한 AMP/ASCO/CAP guideline에 따른 Tier classification]

NGS를 시행하게 되면, 다양한 변이들이 쏟아져 나오게 됩니다. 이때 해당 변이의 판독은 크게 아래와 같은 접근법을 이용하게 됩니다.

1. 기존의 암에서 자주 보고된 알려진 변이인가? Database에 이미 널리 알려진 변이 (매우 소수)
2. Database에 등록 되어 있지는 않지만, 정상 인구 집단에서는 관찰되지 않는 매우 드문 변이인가? (Population genetics 관점에서 allele frequency)
3. 여러가지 in-silico prediction tool이 해당 변이의 deleterious effect를 예측하고 있고, 해당 변이가 단백질의 매우 중요한 3차원적 위치에 있는 경우 (Mutational hot-spot, Functional genetics 관점에서 protein의 기능 및 domain)

 

이러한 접근법에 근거하여 다양한 tool들이 개발되고 있으며, 대표적으로 널리 쓰이는 몇가지 tool들을 소개하며, 이번 포스팅은 마치고자 합니다. NGS 검사를 통해 검출된 변이에 아래의 DB에서 제공하는 다양한 정보를 annotation하고, 이에 근거하여 driver mutation을 예측하게 됩니다.

I. COSMIC (https://cancer.sanger.ac.uk/cosmic)

Wellcome Sanger Institute에서 제공하고 있는 암 유전체 관련 DB입니다. 보통 개별 변이마다 DB에 등록되면서 COSMIC ID가 부여되는데, 가장 방대한 DB를 구축하고 있어서 새로운 변이들을 확인할 때 가장 먼저 살펴보게 되는 DB입니다.

 

II. cBioPortal (http://www.cbioportal.org/)

다양한 암종과 TCGA cancer genome 데이터를 기반으로 하여, 보고된 다양한 mutation에 대한 정보들을 제공하고 있습니다. 대표적이고 유명한 paper들에 사용된 cancer genome DB를 포함하고 있고, 실제 유전자들의 functional domain과 hot-spot 정보들을 함께 제공하고 있어서 유용하게 이용할 수 있는 DB입니다.

 

III. OncoKB (https://oncokb.org/)

Memorial Sloan Kattering Cancer Center에서 구축한 DB로 조금 더 임상적으로 중요한 변이들이 명확하게 curation 되어 있습니다. 임상적으로 중요하고 근거 수준이 명확한 변이들을 Level에 따라서 잘 정리한 장점이 있으나, 변이 데이터는 상대적으로 조금 빈약한 편입니다.

 

IV. Cancer Genome Interpreter (https://www.cancergenomeinterpreter.org)

이미 구축된 다양한 DB와 기존 논문 보고 데이터들을 통합하여, 변이들의 driver mutation 여부를 종합적으로 잘 판독해주는 툴로 유용하게 사용할 수 있습니다. 다만, 프로그램이 공개되어 있지 않고 서버에 직접 본인의 데이터를 업로드하여야 하고 한번에 업로드할 수 있는 변이의 수가 5,000개로 제한되어 있는 점은 단점이라고 할 수 있습니다.

 

V. MutaGene (https://www.ncbi.nlm.nih.gov/research/mutagene/)

가장 최근에 개발된 tool로 Python package도 제공되어 있어, 따로 서버에 자료를 올리지 않고 설치해서 바로 사용할 수 있는 장점이 있습니다. Driver mutation 예측 외에도 Mutational signature 분석도 함께 할 수 있어서, 유용한 정보를 제공하고 있습니다.

 

위의 내용을 살펴보면, 아직까지 완벽한 변이 판독 방법은 없구나 하는 것을 느끼게 됩니다. 사실 이전에 약물 유전자와 관련된 연구에 대해 포스팅 했었는데, 비슷한 연구가 암 관련 유전자에 대해서도 함께 진행 중입니다. 따라서 최근의 연구 추세는 다양한 변이의 임상적 판독을 위한 충분한 정보를 제공할 수 있는 대용량 변이 판독 방법에 집중되고 있으며, Functional genomics 분야의 큰 부분을 차지하며 연구비가 몰리고 있는 상황입니다.

관련 포스팅 보기>

약물유전체 정밀의료의 실현, F-CAP 프로젝트

유전자 변이의 해석: 대용량 기능 검사의 필요성

 

[References]

Li, Marilyn M., et al. “Standards and guidelines for the interpretation and reporting of sequence variants in cancer: a joint consensus recommendation of the Association for Molecular Pathology, American Society of Clinical Oncology, and College of American Pathologists.” The Journal of molecular diagnostics 19.1 (2017): 4-23.

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 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)라고 하는 프로그램으로 진행하게 됩니다.

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

 

II. 정렬된 정보를 binary format으로 변경하기 (SAM to BAM): 위의 과정을 거친 SAM 파일은 여전히 용량이 매우 큽니다. 따라서, 용량을 절약하기 위해서 컴퓨터가 이해하는 2진법의 binary 형식으로 변경하면서 용량을 줄이고 연산 속도를 올립니다. 이 과정을 거친 파일이 흔히 말하는 1차적인 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 검사: Whole Genome & Exome, Targeted Sequencing 비교

임상 검사 및 연구 목적의 검사로 시행하는 염기서열 분석 방법은 기존의 생거 시퀀싱에서 차세대 염기서열 분석법 (Next-generation sequencing; 이하 NGS)로 빠르게 바뀌어서, 이제는 대부분 NGS로 생산된 엄청난 유전체 데이터가 쏟아지고 있습니다. 하지만 이러한 데이터의 생산보다 더 중요한 것은 어떠한 목적과 목표를 가지고 생산된 데이터를 활용할 것인가에 있다고 하겠습니다. 그래서 이번 포스팅에서는 NGS 기술을 이용하여 DNA 시퀀싱을 하는 대표적인 3가지 검사법에 대해서 비교 및 정리하고자 합니다.

관련 포스팅 보기>

NGS 타깃 시퀀싱 패널 검사의 분석 및 해석시 고려할 사항

NGS Targeted Gene Panel 디자인을 위한 고려 사항

Exome sequencing을 위해 고려할 요소들: capture kit와 target coverage 선택

SNP array와 array CGH의 원리 및 UK Biobank Array, Korean Chip

 

I. Whole Genome Sequencing (WGS)

유전체 전체의 시퀀스를 모두 검사하는 방법입니다. 유전체 전부를 보기 때문에, 개별 시퀀스의 깊이 (depth)는 낮아지지만, 이론적으로 유전체 모든 영역의 SNP (Single Nucleotide Polymorphism), INDEL (Insertion & Deletion), SV (Splicing variant), CNV (Copy number variant) 등을 확인할 수 있습니다 (물론 short-read에서는 repetitive 영역을 모두 보는 것은 불가능합니다). 단, 검사 단가가 비싸고 생산된 유전체의 크기가 매우 커서 분석 및 저장에도 많은 비용이 들어갑니다. 그러나 non-coding 영역의 regulatory variant를 발굴하고, 전체 유전체 영역을 동일하게 가정하는 Genome-wide Null 분석이 가능한 장점이 있습니다.

 

II. Whole Exome Sequencing (WES)

유전체 중에서 단백질을 직접 코딩하는 엑손 영역 (Exome; 엑솜)의 유전체만을 분석하는 방법입니다. 사람의 경우, Exome은 전체 유전체의 2% 미만이기 때문에 WGS보다 생산되는 데이터의 크기가 작아, 저장 및 연산 용량도 줄어들고 가격도 저렴합니다. WGS을 분석하더라도 많은 경우, non-coding 영역의 변이는 해석이 어렵고 크게 의미 있는 결과를 얻는 경우가 많지 않기 때문에 WES이 더 경제적이고 더 높은 depth를 얻을 수 있는 장점이 있습니다. 하지만, 타겟 영역을 제외한 부위 (Intergenic 또는 Intron 영역의 regulatory variant)는 검출이 어렵습니다.

 

III. Targeted Gene Sequencing (TGS)

특정 질병 또는 목적에 부합하는 유전자들로만 구성된 Customized 패널을 구성하여 검사하는 방법입니다. 상대적으로 좁은 영역의 유전체만을 대상으로 하기 때문에, 적은 용량의 데이터를 생산해도 높은 시퀀싱 Depth를 얻을 수 있고, 가격적으로도 가장 저렴합니다. 효율적으로 원하는 유전자들로 입맛에 맞게 유전자 패널을 구성하여, 데이터의 연산 및 저장 용량도 줄일 수 있습니다. 대표적으로 약물 대사에 관련된 유전자로 구성된 약물 유전체 패널, 암 발생과 관련된 유전자들로 구성된 고형암 패널 등이 있습니다. 타겟으로 하는 유전자가 명확하고 보고자 하는 영역이 명확한 경우에 가장 적절하여, 임상 검사실에서 가장 많이 사용하는 검사 방법입니다.

[유전체 검사의 범위에 따른 가격, 용도 및 도구] 더 폭넓은 유전체를 검사할수록 가격은 일반적으로 더 비싸지고, Depth는 낮아지게 됩니다. 검사 목적에 따라 적절한 modality를 활용하는 것이 중요합니다.

위의 3가지 검사법을 결정해야할 때 고려해야할 요소는 크게 다음과 같습니다.

1. 검사 목적: 진단용인가? 연구용인가? 연구용이라면, 연구 대상의 유전자가 제한된 연구인가 또는 새롭게 발굴하는 것이 목적인가?
2. 분석 대상: 특정 유전자에 한정할 것인가? 해석이 용이한 엑손 영역에서 새로운 유전자를 발굴할 것인가? 전체 유전체에서 탐색적으로 연구할 것인가?
3. 검체 이용 및 시퀀싱 깊이: 어떠한 검체를 이용할 것이며, 시퀀싱 깊이는 어느 정도가 적절한가?
4. 검사 비용 및 분석 능력: 시퀀싱 결과 생산된 데이터의 크기가 매우 크기 때문에 이것을 분석하는 컴퓨팅 파워 및 용량, 전체적인 비용도 고려해야 합니다.

올해부터 한시적으로 임상 검사의 목적으로는 타겟 시퀀싱 패널 (TGS)에 대해서 보험 급여가 인정되어 임상 검사실에서 시행되고 있고, WES 및 WGS의 경우에는 진단의 목적보다는 연구 목적으로 새로운 후보 유전자 또는 영역을 발굴하는 목적으로 많이 사용하고 있습니다. 위의 3가지 외에도 최근에는 Clinical Exome Sequencing이라고 하여, 전체 엑손 영역 중에서 임상적으로 질병과 연관된 유전자들로만 구성된 일종의 광범위 타겟 시퀀싱 패널과 같은 검사도 WES보다 좀 더 저렴하게 검사가 가능합니다.

또한, 연구자의 입장에서 비용 만을 생각하면 시퀀싱보다는 micro-array 기반의 검사가 더 저렴하고 간편한 경우도 있기 때문에 (e.g. 한국인칩 ㅠㅠ), 목적에 적합한 검사를 선택하는 것이 중요하다고 할 수 있겠습니다. Array 검사법과 시퀀싱 검사법의 비교는 나중 포스팅에서 다루도록 하겠습니다.

 

[References]

https://blog.genohub.com/2015/02/21/whole-genome-sequencing-wgs-vs-whole-exome-sequencing-wes/

https://blog.genohub.com/2016/10/24/targeted-gene-panels-vs-whole-exome-sequencing/

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

The landscape of Human Genome Project (National Human Genome Research Institute 에서 발췌)

인간 게놈 프로젝트 (Human Genome Project; HGP)는 인간의 모든 유전자 염기 서열을 밝히는 것을 목표로 1990년 처음 시작되었습니다. 그리고 당초 목표보다 2년 빠른 2003년도에 목표를 완수하게 되었죠. 이 당시만 해도, 인간의 유전체 암호가 모두 해독되어서 마치 모든 유전병을 정복가능하게 될 거라고 기대가 매우 컸습니다. 그렇지만 당시의 기대와 흥분도 잠시, 지금까지 무려 14년이 흘렀지만 아직도 게놈 프로젝트는 현재 진행형입니다. 이번 글에서는 게놈 프로젝트의 역사와 의의, 그리고 앞으로 유전체 연구 진행 방향에 대해서 논의해 보겠습니다.

사실 게놈 프로젝트를 말하는데 있어서, 염기서열 분석 기술의 발전을 떼놓고는 말하기가 어렵습니다. 하지만 각 세부적인 분석 기술 소개하는 것 만으로도 엄청난 분량이 되어버리므로, 이번 글에서는 간단히만 언급하겠습니다.

1958년 노벨 화학상 수상자, Frederick Sanger (1918 – 2013). 그의 이름을 딴 생거 시퀀싱 방법은 아직까지도 염기 서열 분석 방법의 고전적 분석 표준(Gold standard)으로 인정받고 있습니다.

유전자의 염기 서열을 분석하는 가장 고전적인 방법은 생거 시퀀싱(Sanger sequencing)법 입니다. 영국의 생화학자였던 프레데릭 생거에 개발된 방법은 아직까지도 염기 서열 분석의 표준 방법(Gold standard)로 인정되고 있습니다. 그러나 기술의 발전으로 점점 더 빠르고 정확하게, 대용량의 유전 정보를 시퀀싱하는 것이 가능하게 되었고, 최근에는 이를 두고 차세대 염기 서열 분석법 (Next-generation sequencing, NGS) 라고 부르고 있습니다. 사실 우리는 이미 NGS의 시대에 살고 있기 때문에 엄밀한 의미에서 차세대 라는 말은 맞지 않습니다. 따라서 일부에서는 3세대 염기서열 분석 (Third generation sequencing, TGS) 또는 대용량 염기 서열 분석 (Massively parallel sequencing, MPS)라고 부르기도 합니다. 하지만, 편의상 널리 NGS로 통용되고 있습니다. 예전에는 염기 서열을 밝히는데 수년이 걸렸던 인간 게놈 프로젝트도, NGS 기술의 발달로 이 글을 작성하는 지금도 매우 저렴한 가격에 며칠이면 가능해졌으며, 소요 시간과 가격은 점점 감소하고 있습니다. 그렇다면 이렇게 저렴해진 유전체 서열 분석 기술이 의미하는 바는 무엇일까요?

2000년대 초, 백만불이었던 전체 유전체 분석 가격은 2017년 현재 이미 1,000불 아래로 떨어졌습니다. 즉, 개인이 100만원의 비용이면 본인의 모든 염기 서열 정보를 알 수 있는 시대에 도달했다는 의미가 됩니다.

사실 인간 게놈 프로젝트 이후에 진행된 프로젝트는 1000 게놈 프로젝트 (1000 Genome Project)가 있습니다. 1000 게놈 프로젝트는 1000명의 사람들의 전체 유전체를 분석해서 서로 차이를 보이는 염기 서열에 무엇이 있고, 이것이 개인마다 어떻게 다름으로써 개인별 특징을 나타내는지를 찾고자하는 첫 시도였다고 볼 수 있죠. 그리고 현재 알게된 사실은 사람들은 무수히 다른 단일 유전자 변이 (Single Nucleotide Variant, SNV)를 가지고 있으며, 이러한 SNV가 각기 어떻게 작용을 해서 서로 다른 형질을 보내는지에 연구의 포커스가 맞춰진 상태입니다. 인간이 가지고 있는 형질과 질병은 무수히 많습니다. 그리고 인간의 유전체도 SNV를 포함해서 매우 다양하며, 무궁 무진한 상호 작용과 조절을 받고 있습니다. 즉, 인간 게놈 프로젝트 초기에 기대했던 목적을 완전히 이수하려면 이러한 모든 유전체의 발현과 조절, 그리고 개개인의 유전 정보와의 관계를 밝히는 것이 필수라는 것이죠. 1000명에서 시작된 프로젝트는 이제 나라와 인종별로 10만명, 100만명 수준으로 확대되고 있으며, 앞서 언급한 유전체와 형질간의 관계를 확인하기 위한 데이터 수집 과정입니다.

즉 앞으로의 유전체 프로젝트와 연구 방향은 이러한 형질 또는 질병과 유전 정보 간의 관계를 파헤치는데 집중될 것입니다. 유전자 지도 완성 뿐 아니라, 유전자 지도 안의 각 위치가 어떠한 역할을 하고, 어떤 형질과 질병을 일으키는지 완벽하게 이해를 해야한 엄밀한 의미의 정밀 의료가 실현 가능하다는 것이죠.  유전체 정보는 매우 방대하고, 이를 분석하는데는 매우 많은 시간과 노력이 필요합니다. 따라서 뒤이어 발전하게 된 것이 이러한 유전 정보를 분석하는 생물 정보학 (Bioinformatics)입니다. 더 나아가 최근에는 이러한 유전 정보를 빅데이터로 간주하여, 인공지능 방법론을 활발하게 적용하고 있는 상황입니다. 정밀 의료를 논의하는데 있어, 유전학, 생물 정보학, 그리고 인공 지능 등을 함께 이해할 필요가 있는 대목입니다.

유전 정보 데이터를 바탕으로 한 정밀 의료 를 실현하기 위해서는 생물 정보학과 같은 분석 도구가 필수적이며, 유전체 데이터가 매우 크기 때문에 이를 효과적으로 처리하기 위해 빅데이터와 인공 지능 기술 적용이 크게 각광 받고 있습니다.

궁극적으로 정밀 의료 시대에는 모든 개인이 자신의 유전 정보를 보유하게 될 것이고, 앞에서 알게된 지식을 바탕으로 사는 시대에 있게 될 것입니다. 그리고 지금은 연구자들이 그러한 지식의 틈을 메꿔가는 과정이라고 볼 수 있겠죠.

유전 정보 분석 및 생물 정보학, 그리고 빅 데이터와 인공지능. 각각의 세부적인 내용은 나중에 더 자세히 알아보도록 하고 이번 글은 여기서 마치겠습니다.