![[07.15(수)]_제14차_CBMC_희귀질환_센터_화상_컨퍼런스](https://2wordspm.wordpress.com/wp-content/uploads/2020/07/07.15ec8898_eca09c14ecb0a8_cbmc_ed9daceab780eca788ed9998_ec84bced84b0_ed9994ec8381_ecbba8ed8dbceb9fb0ec8aa4.png?w=563&h=813)
COVID-19는 학회와 강의의 풍경도 많이 바꾸어 놓은 것 같습니다. 개인적으로는 Zoom을 이용하여 화상 인터뷰 및 미팅 등은 한 적이 있는데, 온라인 강의는 이번이 처음입니다.
의국 선배님으로 부터 NGS에 대한 기초적인 개괄에 대한 강의를 부탁받아, 분당 차병원 임상과 선생님들을 대상으로 강의 자료를 준비하였습니다. 강의 자료가 필요하신 분들을 위해 블로그에도 업로드 합니다. (사용은 자유롭게 하시되, 출처만 명시해주시기 바랍니다.)
관련 포스팅 보기>
바이오 연구자를 위한 Genome Browser 비교 및 활용
NGS 결과의 임상 적용: Genotype-phenotype correlation
NGS 분석 파이프 라인의 이해: GATK Best Practice
Annovar: Population frequency, in-silico prediction tool 및 기타 database 활용
강의 자료 다운로드 > CBMC conference
지난 글에 이어서, 생성된 Bam 파일로부터 변이들을 읽어 들이고, haplotype call을 하는 두번째 파트의 code를 정리해보겠습니다. 아래 코드는 GATK 4.1.3 버젼을 기반으로 작성되었습니다. GATK 버젼에 따라서 조금씩 Tool과 명령어에 차이가 있습니다. 전반적인 흐름은 아래와 같습니다.
관련 포스팅 보기>
NGS DNA-seq pipeline: GATK Best Practice Code – Part1. Fastq to Bam
해당 위치의 변이와 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
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
GATK4 (Genome Analysis Tool Kit)로 넘어오면서, NGS 시퀀싱 분석을 위한 파이프라인이 많이 개선 및 간소화된 것 같습니다. 덕분에 저도 최근 연구실에 구축되어 있던 파이프라인도 새롭게 뜯어고쳤는데, 이 참에 전반적인 분석을 위한 코드를 정리해볼까 합니다. 추후에 비슷한 파이프 라인을 구축하고자 하는 분들께 도움이 되었으면 하고 피드백도 환영합니다.
관련 포스팅 보기>
우선 기본적으로 모든 작업은 리눅스 (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
위의 프로그램이 설치가 끝나면, path 등록을 해서 리눅스 상에서 명령어를 입력하면 잘 실행되는지 확인합니다. 마지막으로 Reference Genome 파일을 다운로드 합니다. hg19 또는 hg38 버젼에 맞게 파일을 모두 다운로드하여 위치시키면 기본적인 준비들이 끝납니다. (파일은 UCSC genome browser에서 다운 가능)
아래 코드는 GATK best practice의 Data pre-processing step에 해당하는 내용을 정리한 것입니다. 전체 Work-flow는 아래와 같습니다.
GATK Best practice: Data pre-processing for variant discovery
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을 진행한다.
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으로 대체해서 진행이 가능하다.
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 파일을 생성하는데 이용된다.
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
5월 29, 30일간 K 호텔에서 진행된 대한진단유전학회 학술대회에 참석하였습니다. 기억이 잊혀지기 전에 몇가지 인상 깊었던 내용들에 대해서 정리하고 가고자 합니다.
우선 많은 임상과에서 연자 선생님들이 초청되어서, 실제 진료에 NGS 검사 결과를 적용한 사례들을 발표해주셨는데, 아직까지는 갈 길이 멀어보이지만, 그래도 검사를 통해서 더 나은 진료와 치료를 받을 수 있었던 케이스들을 통해서 정밀 의료라는 큰 흐름은 계속 발전해 나갈 것으로 기대가 됩니다. 특히, 어느 정도 경험이 쌓인 분들은 변이 판독과 Genotype-Phenotype correlation에 대해서 많은 고민을 하고 계신 것 같습니다. 그러나 여전히 검출되는 수많은 변이의 임상적 의미에 대한 판독 부분은 시간이 더 필요할 것으로 보입니다.
관련 포스팅 보기>
특히, 분자 유전학적 메커니즘에 따라서 치료제가 개발된 경우에 치료 효과를 기대할 수 있는 환자들을 선별하고 해당 환자들이 눈에 띄는 임상적 호전을 보이는 몇몇 대표적인 케이스들은 상당히 인상 깊었습니다. 다만, 아직까지는 검사 대비 효용성이나 검사를 위한 가이드 라인 수립, 검사 후 치료에 대한 보험 적용 문제 등 넘어야할 산은 매우 많아 보였습니다.
더불어 암 환자들의 혈액에 존재하는 매우 미량의 DNA, 일명 순환 종양 DNA (Circulating tumor DNA; CtDNA),를 검출하여 조기 진단 및 치료 반응 추적 등에 활용하고자 하는 연구가 많이 진행되고 있는데, 이번 학회에서는 CtDNA에 대한 많은 관심과 진행되고 있는 연구에 대해서 살펴볼 수 있었습니다. 특히나 암 환자의 경우 조직을 이용하여 진행되었던 많은 검사들의 무게 중심이 궁극적으로는 혈액 검사를 이용하는 쪽으로 옮겨가고 있다는 느낌이었습니다. CtDNA의 경우는 워낙 핫한 분야이기 때문에 다음에 기회가 되면, 관련 내용을 리뷰해 보도록 하겠습니다.
이외에도 Human Microbiome 분야와 post-GWAS 분야 (SNP 발굴을 넘어서, meta-analysis 및 complex trait에 대한 polygenic risk score 수립을 통한 예측) 등도 주로 다루어졌습니다. 개인적으로는 전반적으로 최근에 많은 관심이 쏟아지는 주제에 대해서 고르게 구성된 유익했던 학회였습니다.
최근들어 제 블로그의 방문자들이 눈에 띄게 증가했습니다. 대부분 구글 검색을 통해서 유입되는 분들인데, 검색어 통계를 보니 NGS 관련 내용에 대해서 검색하다가 들어오는 사람들이 많았습니다. 그래서 저에게는 쉬운 내용들이라 하더라도, 일반인들이 궁금해하는 내용에 대해서도 정리해두어야 겠다는 생각이 들어서 이번 포스팅은 NGS 장비를 통해서 생산된 데이터를 어떻게 분석하는지, 전반적인 파이프 라인에 대해서 개념을 소개하는 포스팅을 올려볼까 합니다.
우선, 차세대 염기 서열 분석법 Next-generation sequencing (이하 NGS)은 다양한 이름으로 불리는데, 이미 널리 사용하는 기술이기 때문에 현재는 차세대 기술이라고 보기 어렵습니다. 따라서, NGS라고 부르는 건 misnomer라고 할 수 있죠. 좀 더 정확한 명칭으로는 Massively parallel sequencing, High-throughput 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
I. 표준 유전체 서열에 read를 정렬하기 (FASTQ to SAM): 보통 NGS를 통해서 생산되는 개별 Read의 개수는 백만개 이상의 단위가 됩니다. 이때, 개별 read의 정보는 FASTQ 파일로 저장되고 그 크기는 수십~수백 Gb 단위가 됩니다. 각각의 생산된 read는 이미 알려진 인간 표준 유전체 서열에 가장 잘 맞는 위치에 정렬시키게 됩니다. 예로, 100 bp read를 기준으로 하여, 100개의 서열이 모두 일치하는 경우는 거의 유일하게 되므로, 해당 위치에 잘 찾아가게 됩니다. 다만, 반복 서열이나 특이적이지 않은 서열의 read는 제대로 mapping이 되지 않을 수가 있는데, 이는 Illumina 방식이 가진 한계점입니다. 보통 이 과정은 BWA (Burrow-Wheeler Aligner)라고 하는 프로그램으로 진행하게 됩니다.
II. 정렬된 정보를 binary format으로 변경하기 (SAM to BAM): 위의 과정을 거친 SAM 파일은 여전히 용량이 매우 큽니다. 따라서, 용량을 절약하기 위해서 컴퓨터가 이해하는 2진법의 binary 형식으로 변경하면서 용량을 줄이고 연산 속도를 올립니다. 이 과정을 거친 파일이 흔히 말하는 1차적인 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/
암유전체 분석은 크게 SNV/INDEL 수준의 Mutation 분석과 Chromosome/CNV 수준의 Somatic Copy number Alteration (SCNA) 분석으로 나눌 수 있습니다. 과거에 SCNA는 주로 SNP array 또는 Array CGH과 같은 Microarray를 이용하여 시행하였지만, 최근에는 NGS 데이터를 활용하여 2가지 분석을 모두 시행할 수가 있습니다. 이번 포스팅은 NGS 데이터를 활용하여, SCNA를 분석하는 Genomic Identification of Significant Targets in Cancer (GISTIC) 분석 방법에 대해서 정리해보고자 합니다.
관련 포스팅 보기>
ftp://ftp.broadinstitute.org/pub/GISTIC2.0/GISTICDocumentation_standalone.htm
Segmentation File (-seg)
(1) Sample (sample name)
(2) Chromosome (chromosome number)
(3) Start Position (segment start position, in bases)
(4) End Position (segment end position, in bases)
(5) Num markers (number of markers in segment)
(6) Seg.CN (log2() -1 of copy number)
원래 GISTIC은 Array 기반으로 개발된 프로그램이기 때문에, Probe 정보를 받아들이게 되는데 NGS 데이터는 Probe 정보가 없습니다. 따라서, 타겟 영역의 엑손 하나 하나를 일종의 Probe로 간주하고 데이터를 변환하여 넣어주면 분석을 할 수가 있습니다.
NGS를 통해 생산된 Bam 파일을 이용하여, 타겟 영역의 Copy number 정보를 구하고 이 데이터를 활용해서 아래와 같이, 1차적으로 Segmentation을 해줍니다. 그리고 여러 샘플의 이러한 정보를 합쳐서, 통계적으로 유의미한 CN 수 변화가 발생한 곳을 검출하게 됩니다. 최근에는 유용한 R package가 많이 있는데, DNA copy를 이용한 아래 코드는 segmentation과 데이터 변환에 유용하여 함께 올립니다.
![[참석 후기] 2019 대한진단유전학회 학술대회](https://2wordspm.wordpress.com/wp-content/uploads/2019/05/sub_top.jpg?w=825&h=170&crop=1)
#Preparing CopywriteR output for GISTIC 2.0 analysis
library(DNAcopy)
load("/PATH/TO/segment.Rdata")
segmentation.values <- segment.CNA.object$output
colnames(segmentation.values) <- c("Sample", "Chromosome", "Start Position", "End Position",
"Num markers", "Seg.CN")
write.table(segmentation.values, file = "/PATH/TO/segmentation_values.tsv", quote = FALSE,
row.names = FALSE, sep = "\t")
markers <- data.frame(paste(segment.CNA.object$data$chrom, segment.CNA.object$data$maploc,
sep = ":"),
segment.CNA.object$data$chrom, segment.CNA.object$data$maploc)
colnames(markers) <- c("Marker Name", "Chromosome", "Marker Position")
write.table(markers, file = "/PATH/TO/markers.tsv", quote = FALSE, row.names = FALSE,
sep = "\t")
위와 같은 데이터 변환 후 Input 데이터 변환이 끝나면, GISTIC 분석을 위한 모든 준비가 끝나게 됩니다. GISTIC은 Matlab 기반 프로그램이지만, 다행히 cloud를 통해 GenePattern에서 웹기반으로도 이용할 수 있습니다. 마지막으로, GISTIC을 이용하여, 4,934개의 암 샘플을 분석한 Nature genetics의 논문을 소개하며, 포스팅을 마치도록 하겠습니다.
[GISTICS 분석 결과] 유의하게 amplification 또는 deletion이 발생한 위치에 존재하는 Tumor driver gene을 발굴함으로써, 암 발생에 대한 연구를 할 수 있습니다.
최근 많은 유전체 연구가 유전체의 3차원적인 구조와 직접 단백질을 coding 하지 않는 non-coding region의 역할과 질병 발생 메커니즘에 이뤄지면서, 유전체에 대한 이해의 폭이 넓어지고 있습니다. 최근에는 과거의 단순하게 A,T,G,C로 이루어진 염기 서열의 2차원적인 시각에서 벗어나, 유전자의 발현 기작은 훨씬 복잡하게 이뤄진다는 것을 알게 되었습니다. 사실 이러한 이유 때문에 GWAS 또는 NGS를 이용하여 유전자의 염기 서열을 확인하는 것으로는 유전형과 실제 표현형 간의 괴리가 크지 않았나 하는 생각입니다. 이러한 이유 때문에, 저는 DTC 검사로 단순한 몇개의 SNP을 검사해서 개인의 다양한 표현형을 예측해준다는 회사들에 회의적입니다. 오늘은 유전체 영역 중 non-coding region 에 위치하면서, 실제로 유전자들의 발현을 조절하는 Enhancer와 Super-enhancer의 개념을 정리하는 포스팅을 남기고자 합니다.
관련 포스팅 보기 > DTC 유전자 검사의 딜레마: 과학과 산업 사이
우선 들어가기에 앞서, Enhancer의 개념에 대해서 살펴보겠습니다. 유전자의 발현은 유전자의 위쪽에 존재하는 Promotor의 transcription binding site에 신호 전달 물질을 통해서 이뤄지는데, enhancer는 이러한 promotor와 작용함으로써 유전자 발현을 조절하는 부위입니다. Super-enhancer는 그 이름에서 유추할 수 있듯이, 이러한 enhancer의 작용이 특별히 더 강한 그룹을 의미하는데, 더 자세한 정의는 아래에서 살펴보겠습니다. 조직별로 유전자 발현 패턴은 서로 상이하게 되는데, 이러한 유전자의 발현 양상의 차이는 결국 다양한 enhancer들의 작용 기작과 관계되어 나타난다는 것이 현재까지의 관찰입니다.
ChIP-seq에 대한 포스팅 보기 > Non-coding 영역의 GWAS 신호 해석: 3C based method
Pott, Sebastian, and Jason D. Lieb. “What are super-enhancers?.” Nature Genetics 47.1 (2015): 8.
저는 작년 2월부터 1년 반정도의 기간을 약물유전체 연구를 하고 있습니다. 지도 교수님이신 이민구 교수님과 다양한 약물 반응에 대한 유전적 바이오 마커를 발굴하는 연구를 하고 있는데, 생각보다 쉽지가 않고, 좋은 결과가 나오지 않고 있습니다. 그에 비해 최근에 선천성 기형의 일종인 두개골 조기유합증이라는 희귀질환에 대해 성형외과 및 신경외과와 공동연구에도 참여하고 있는데, 많은 환자들의 유전적 원인들을 확인할 수 있었습니다. 유전적 소인과 형질 간에는 어떠한 연관이 있는 것일까요? 이번 글은 흔히 말하는 Common disease 와 Rare disease 의 차이와 더불어, 지난 1년반정도의 기간을 약물 유전체 연구를 하며 느낀 점들과 왜 약물 유전체 연구가 어려운지에 대해서 정리해보고자 합니다.
기본적으로 약물 유전체 연구는 크게 여러 사람들이 동일한 약물을 먹었을 때 혈중 유효 농도가 다양하게 나타나는 것에서, 어떤 유전적 차이가 이러한 약물 대사에 기인하는지와 부작용 발생 유무의 위험도를 예측할 수 있는 유전적 바이오마커가 있는지에 관심을 갖춰 연구되고 있습니다.
관련 포스팅 > 약물 유전학은 왜 정밀의료에서 중요한가?
I. 약물 반응은 복합 형질 (Complex trait)이다 : 기본적으로 약물의 대사 과정에는 다양한 약물 효소가 관련합니다. 또한 약물이 흡수되어 배출되기까지의 대사 과정 (ADME) 또는 약동학 (Pharmacokinetics) 과정에는 다양한 요소들이 관여하기 때문에, 한 두가지 유전적 소인이 형질에 결정적 차이를 나타내기 어렵습니다. 복합 형질로 가장 많이 연구되는 질병 중 하나가 2형 당뇨병 (Type 2 Diabetes mellitus; T2DM)인데, 당뇨병 발생의 원인과 그 유전적 요인에 대해서 많은 연구가 진행되었지만 여전히 속 시원한 유전적 원인에 대해서는 알지 못하고 있습니다. 특히 이러한 복합 형질에서 발굴된 유전적 마커들은 형질의 차이에 기여하는 정도가 매우 작아서, 대부분의 효과 크기 (Effect size)가 매우 작습니다. 그래서 그나마 연구가 잘되고 결과가 잘 나오는 것은 효과 크기가 매우 큰 한 두가지의 유전적 인자가 약물의 부작용 발생 유무에 영향을 미치는 경우입니다.
II. 약물 반응의 측정 자체가 어렵다 : 체내 약물 대사능에 영향을 주는 유전적 인자를 확인하고자 하는 연구의 경우, 일단 환자에서 해당 약물 농도 측정 자체가 매우 어렵습니다. 현실적으로 환자들에게는 의사들이 체중이나 대사능 등을 고려하여 약을 처방하기 때문에 복용한 약물의 양도 간격도 전부 달라지게 되며, 약물 농도라는 것도 매우 변동성이 심하기 때문에 언제 채혈하였는지, 다른 약과 함께 복용하였는지 (drug-drug interaction), 음주 & 흡연 여부, 성별 등 다양한 요소에 영향을 받게 됩니다. 기본적으로 이러한 요소들에 대한 명확한 통제가 어렵고, 보정을 한다고 하더라도 그 측정 약물 농도가 명확하게 그 사람의 약물 대사능을 대변하지도 못합니다. 즉, 처음부터 얻어지는 정보 자체에 매우 큰 변동성이 있기 때문에 해당 데이터와 유전적 정보 간의 연관성을 찾으려고 해도, 그 영향이 명확하게 큰 경우가 아니면 연관성을 찾기가 매우 어렵습니다.
III. 약물 대사 경로에는 다양한 대체자가 존재한다. : 이 세상에는 정말로 다양한 약물이 존재합니다. 기본적으로 약물은 간에서 대사되어 신장을 통해 배설된다고 알려져 있습니다만, 약물 개별로 보면 어떤 약물이 정확하게 어떠한 효소에 의해 대사되어 어떠한 형태로 배설되는지, 명확하게 알려져 있는 약물은 그리 많지 않습니다. 희귀 질환의 경우에는 생명에 필수적인 역할을 하는 어떠한 유전자에 문제가 생겨서 바로 질환으로 나타나는 경우가 많습니다. 이는 해당 유전자가 만들어내는 단백질이 중요한 역할을 하고, 다른 유전자가 대신 기능을 해주지 못하기 때문입니다. 반면에 약물 유전자가 만들어내는 약물 효소의 종류는 워낙 다양해서 한 두가지 효소에 문제가 생긴다고 하더라도, 비슷한 다른 효소가 이러한 역할을 대신해주게 됩니다. 그리고 대사 경로 자체가 한가지 방향으로만 정해져 있는 것이 아니라, 어떠한 길이 막히면 다른 길로 돌아갈 수 있는 대체 경로가 존재하게 됩니다. 즉, 약물 대사능은 한가지 유전자와의 1:1 대응이 아니라, 다수의 효소들이 관여하여 복합적으로 나타나기 때문에 동시에 고려해야할 요소들이 많아지게 됩니다. 이를 유전학적으로 나타내보면 다음과 같습니다.
IV. 연구 방법의 한계 : 유전적 바이오 마커 발굴의 연구 방법으로 많이 사용하고 있는 것이 SNP array chip 또는 NGS를 통한 시퀀싱입니다. SNP array는 주로 GWAS 연구에 사용하기 때문에 인구집단에 흔하게 존재하는 common variant 연구에 사용하고, NGS 시퀀싱은 유전자의 개별 변이까지 모두 확인하기 때문에 rare variant 발굴에 사용하게 됩니다. 그러나 두 연구 방법 모두 한계가 있습니다. 앞에서 언급한 것처럼 복합형질에서 common variant는 그 효과 크기에 대부분 매우 작기 때문에 GWAS 연구로는 새로운 마커의 발굴이 쉽지 않은 편입니다. 반면 Rare variant 발굴에 유리한 NGS 방법으로는 rare variant를 발굴하여도 그 변이의 해석이 쉽지 않고, 더불어 통계적으로 의미 있는 결과를 얻기 위해 필요한 n수가 매우 커서 현실적으로 연구가 어렵게 됩니다.
관련 포스팅 >
[유전자칩 비교] SNP array와 array CGH, 그리고 한국인칩
위에서 언급한 여러가지 이유들로 인해, 약물 유전체 연구는 정말 어려운 분야인 것 같습니다. 하지만 다른 한편으로는 정밀의료 분야의 발전으로 가장 많은 사람들이 혜택을 볼 수 있는 분야도 약물과 관련된 분야이기 때문에, 그만큼 의미가 크다고 할 수 있겠습니다. 이러한 여러가지 어려운 점에도 불구하고, 열심히 노력하고 있는 연구자들이 함께 좋은 연구가 나올 수 있었으면 하는 바람입니다.
NGS 데이터를 이용하여 환자의 변이를 판독하는데 있어 제일 어려운 부분은 Intron 영역과 splicing site 인 것 같습니다. Exon 영역의 경우는 직접적으로 아미노산 서열에 영향을 주는 부분이기 때문에 Null variant 또는 Missense variant에 따라 어느 정도 예측이 가능하지만, Non-coding 영역인 Intron 영역은 해석하기 어렵기 때문이지요. 그래서 종종 Intron 영역은 배제하고 Coding region만 변이 판독을 하는 경우가 많습니다. 그러나 이럴 경우, splicing site mutation을 종종 놓칠 수 있습니다. 판독에 포함하더라도 실제 검출된 변이가 어떤 영향을 주는지 알기 어려운 경우도 많습니다. 그래서 이번 포스팅에서는 splicing이 일어나는 과정 및 실제로 splicing site에서 mutation이 검출된 예를 통해서 splicing site 변이 해석 방법에 대해서 정리해 보고자 합니다.
위의 과정에서 transcription이 제대로 일어나지 않은 경우, 잘못된 mRNA가 생성될 수 있고 이러한 mRNA의 산물로 잘못된 단백질이 형성되어 환자의 형질이 나타날 수가 있습니다. 아래는 splicing 과정 중에서 기능적으로 중요하여 보존된 영역의 sequence를 보여주고 있습니다. 따라서 일반적으로 NGS 변이 판독 시에는 exon 영역 전후 10bp 또는 50bp 까지도 판독에 포함하기도 합니다. 그러나 많은 경우, 판독이 쉽지 않아서 실제로 mutation을 검출하는 경우는 드뭅니다.
최근에 두개골 조기 유합증 환자의 NGS 결과를 판독하다가 나온 예를 통해 Splicing site 변이를 판독하는 방법을 살펴 보겠습니다. 해당 환자는 TCF12 유전자의 c.1468-7A>G 변이가 heterozygote로 확인되었습니다. 아래 그림과 같이 원래 AA sequence이던 부분이 변이로 인해 AG로 바뀌면서 원래 splicing acceptor site로 작동해야할 부분의 앞쪽이 splicing 되면서 잘못된 transcription이 발생한 case 입니다. 위 그림 (c)의 Cryptic splice site usage에 해당합니다.
이렇게 되면 원래 exon 17 앞의 intron 영역의 CTTTAG sequence가 포함되어, 실제 mRNA에는 Leu(CUU)-Stop(UAG) codon이 포함되고, 결국 해당 mRNA는 inserted stop codon에 의해 exon 16번까지만 발현되는 Stop gain variant와 같은 결과를 보이게 됩니다.
마지막으로 이러한 splicing variant를 simulation 하는 in-silico tool을 소개하면서 포스팅을 마치고자 합니다. 아래 논문에서는 splicing에 영향을 주는 SNV의 효과를 예측하는 in-silico tool에 대해서 소개하고 있는데, scSNV score로 명명하여 여러 컴퓨터 알고리즘을 적용하고 있습니다. 위의 환자의 변이는 ADA score 0.9995 / RF score 0.9739 으로 ADA 또는 RF 알고리즘으로 예측한 결과 모두 영향을 받을 가능성이 아주 높음 (1에 가까울 수록) 을 보여주고 있습니다.
[Reference]
Singh, Ravi K., and Thomas A. Cooper. “Pre-mRNA splicing in disease and therapeutics.” Trends in molecular medicine 18.8 (2012): 472-482. https://doi.org/10.1016/j.molmed.2012.06.006
Xueqiu Jian, Eric Boerwinkle, Xiaoming Liu; In silico prediction of splice-altering single nucleotide variants in the human genome, Nucleic Acids Research, Volume 42, Issue 22, 16 December 2014, Pages 13534–13544, https://doi.org/10.1093/nar/gku1206
Copy Number Variation (CNV)는 Single Nucleotide Variation (SNV)과 더불어, 유전적 다양성을 나타내는 주요한 원인으로 생각되고 있습니다. 유전자 sequence의 염기 하나가 치환된 SNV와 달리 CNV는 유전자 전체 또는 일부의 copy가 중복되거나 줄어들수도 있기 때문에 훨씬 넓은 영역에서 나타나는 Structural variation의 일종입니다.
관련 포스팅 보기 -> 유전학 중요개념 정리: Structural variation 및 Copy-number variation
일반적으로 NGS는 SNV를 보기 위한 목적으로 검사를 시행하지만, 해당 데이터를 활용하면 CNV 분석도 할 수 있기 때문에, 오늘은 NGS 데이터를 활용한 CNV 분석 방법에 대해 포스팅하고자 합니다.
위의 그림은 CNV 분석을 위한 NGS 데이터의 5가지 활용 원리를 나타내주고 있습니다. 그러나 가장 핵심이 되는 원리는 Read depth입니다. Target sequencing과 같이 Read depth가 충분한 경우에, 다른 검체들에 비해 해당 영역의 depth가 월등히 떨어지거나, 또는 월등히 높은 경우에는 해당 영역의 deletion 또는 duplication을 의심할 수 있습니다.
사실 NGS 데이터는 CNV를 목적으로 한 것이 아니라, SNV 검출 목적의 데이터를 부수적으로 활용하는 것이기 때문에 많은 제한점이 있습니다. 따라서, 임상적으로 CNV 검사 목적의 NGS는 권장되지 않으며 적절한 가이드라인도 존재하지 않기 때문에 다양한 Computational tool 들이 개발되어 서로의 장점을 홍보하는 상황입니다. 다음은 다양하게 개발된 대표적인 CNV 검출 tool 들을 정리한 표입니다. 많은 경우 BAM 파일을 활용하는 것을 볼 수 있으며 대부분 R package를 제공하고 있어, 사용이 용이합니다.
![[유전학 중요개념 정리] Enhancer, Super-enhancer](https://2wordspm.wordpress.com/wp-content/uploads/2018/11/molce-40-3-169f1.gif?w=825&h=475&crop=1)
두마디 정밀의료 