[medical] 유전역학-개론 6주차, 차세대 시퀀싱(NGS) 데이터 분석

지난 시간에 NGS에 대해 배웠다.

NGS는 단지 시퀀싱하는 기술 그 자체이고, 이번에 알아야할 것은 시퀀싱한 그 정보를 어떻게 분석을 할것이냐 하는 것이다.

chip array를 가지고 하던지, 생어시퀀서로 하던지, NGS로 하던지 간에 우리가 받는 데이터는 diplotype의 genotype data이다.

무엇을 하던간에 질환과의 연관분석을 하기 위해서는 이 genotype data를 가지고 살펴보면 된다.

 

그럼 NGS 분석이라는 것은 무엇이냐?

연관분석을 하기 전의 '데이터 전처리' 단계라고 생각하면 된다.

DNA조각들을 가지고 시퀀싱을 하게 되면, 우리가 알고있는건 단지 그 DNA 조각들(read)의 정보 뿐이다.

이것들을 통합하고 분석하는 과정이 필요하다.

 

chip array 데이터는 이렇게 복잡하지 않았다. 왜냐하면 단지 이미 genetic marker를 가지고 있고 그 마커에 대해 맵핑만 하면 되기 때문이다.

 

 

지난 시간에 DNA sample을 가지고 library(DNA 조각에 어뎁터를 붙힌 상태)를 만들어 시퀀싱을 하면 DNA 조각들의 염기서열을 읽을 수 있고,

우리가 이제 해야할 것은 이 DNA 조각들을 어떻게 처리할 것이냐 이다.

 

시퀀서를 통해서 우리는 genome과 chromatin, transcriptome 등 다양한 정보를 알 수 있지만, 여기서는 genome에 대한 내용만 다룬다. (빨간 점섬 부분)

 

 

시퀀싱을 목적은 크게 두가지, De novo와 리시퀀싱으로 나누어 볼 수 있다.

목적은 다르지만 30억의 염기서열인 아웃풋은 동일하다.

 

먼저 De novo 시퀀싱이다. 이는 어떤 유기체에 대한 최초의 시퀀싱을 의미한다.

이게 왜 중요하냐면, 보통 시퀀싱을 하면 유기체끼리는 대부분의 유전체정보를 공유하기때문에 그 유기체의 유전체 레퍼런스 정보를 가지고 맵핑을 하는데, 이같은 경우는 레퍼런스할 정보가 없기 때문이다.

그렇기 때문에 1kb이상의 긴 리드가 필요하다. 모나리자의 퍼즐 조각이 크면 클수록 맞추기 쉬운 것과 같은 개념.

휴먼 지놈 프로젝트에서 우리는 드노보 시퀀싱을 했다.

 

반대로 Resquencing은 레퍼런스를 통해 맵핑하는 과정을 통해 보다 수월하게 시퀀싱을 할 수 있다.

보통 대부분의 시퀀싱이라고 하면 이 리시퀀싱을 의미한다.

 

 

아웃풋은 동일하지만 목적이 달라, 그에 따른 프로세싱도 각각 De nono는 어셈블리, 리시퀀싱을 맵핑이라고 한다.

어셈블리는 리드간의 유사성을 통해 서로 맞춰나가면서 전체 지도를 완성하는 것이다.  

맵핑은 보이는것처럼 레퍼런스 지놈을 통해 나의 리드를 맞춰나간다.

 

 

사람의 유전정보 뿐만 아니라 어셈블리는 최근에 코로나 바이러스를 이 드노보 시퀀싱을 통해 완성할 수 있었다.

 

어셈블리 방식은 위 그림처럼 각각의 리드의 유사성을 통해 contig를 완성한다.

contig는 리드보다는 조금 더 긴 시퀀스를 의미한다.

이 contig가 제대로 만들어질려면 리드간의 겹치는 정도가 많아야 그 신뢰도가높아진다.

 

 

위 그림에서처럼 중간 노란색 염기서열 T는 리드간의 7개 겹치기 때문에 7뎁스, coverage 또는 7x라고 한다.

염기서열 각각의 coverage를 구해서 평균을 내는 것으로 그 신뢰도를 평가한다.

high 뎁스를 위해서는 high throuput 시퀀싱을 해야한다.

 

 

어셈블리 프로세스는 어셈블러라는 알고리즘을 통해 수천만개의 리드를 몇백개의 contig로 만든다.

이후 contig 역시도 여러개를 묶어 몇십개의 scaffold라는 단위로 묶을 수 있다.

그럼에도 불구하고 이어지지 않는 부분들이 나온다.

이런 부분은 repeated 시퀀스가 많은 부분이라고 한다.

이런 부분은 'NNNNN'과 같은 것으로 일단은 묶어두고, 휴먼 지놈 프로젝트에서도 이런 'NNNN' 부분들이 엄청 많았다.

최근에 들어서야 이런 부분들을 겨우 완성할 수 있었다고..

 

 

어셈블리가 어려운 이유는 바로 이 반복되거나 비슷한 시퀀스가 너무 많기 때문이다.

그래서 다른 위치에 존재하는 시퀀스를 어셈블러가 동일한 것으로 인식하여 에러가 발생한다.

또한 롱리드로 뽑아내는 것도 어렵고(에러가 많아서), high coverage로 뽑으려면 돈도 많이 든다.

 

이처럼 휴먼 지놈 레퍼런스도 계속 업데이트가 된다. 이 레퍼런스를 만드는 연구그룹은 GRC라고 한다.

지놈 레퍼전스가 계속 업데이트되어 가는걸 지놈 빌드라고 한다.

2003년에 완성하고도 계속 버전업이 되고 있다. 

 

이 지놈 빌드가 중요한 이유는 레퍼런스가 업데이트 된다는 뜻이기 때문이다.

몇 bp인지도 모르는 레퍼런스의 공백 부분이 버전업이 되면 매꿔진다. 이렇게 되면 SNP 같은 포지션 정보들이 뒤로 밀려나기 때문에, 시기에 따른 연구들을 통합하거나 공동연구에 어려움이 생길 수 있다.

툴(리프트오버)이 잘 나와있기때문에 활용을 하면 되고, 이러한 것들이 있구나 하는 것을 아는게 중요하다.

 

 

레퍼런스 지놈은 레퍼런스 어셈블리라고도 불리는데, 이 레퍼런스 지놈은 여러명의 지놈을 합성하여 diploid가 아닌 haploid 형태(DNA 한 줄, 부모 중의 한 allele 정보)로 만들었다.

그러는 이유는 오직 한 종류의 allele 정보만을 제공하여 variant 판정의 기준으로 사용하기 위함이다.

이는 단지 레퍼런스와 나의 DNA가 같은지 다른지만을 구분하기 위함.

 

현재 GRCh의 버전은 38버전의 패치 14번이다.

또 하나 구분해야하는게, reference standard나 reference panels은 population의 비교를 위해 사용한다.

단순히 human이라는 유기체를 비교하는 GRCh 레퍼런스가 아닌, 인종과 관련된 레퍼런스인 것이다.

 

 

1965년부터 어셈블리 지놈의 역사가 시작되었고, 사람뿐만 아니라 강아지, 원숭이 등 과 같은 다양한 유기체에 대한 어셈블리 시퀀싱이 이루어지고 있다.

 

 

대부분의 시퀀싱은 이 리시퀀싱을 의미하고, 이것의 목적은 SNP/INDEL 등과 같은 variants를 찾는 것이다.

 

그리하여 이 NGS 데이터 프로세스라는 것 자체가 리시퀀싱, 맵핑과 관련된 설명이다.

앞선 시간에 시퀀싱 기술에 대해 소개하였고, 그로부터 나온 결과물이 3번의 시퀀싱 로우 데이터이다.

이제 맵핑을 해야하는데, 맵핑할때 필요한 것은 레퍼런스 시퀀스인데 이는 NCBI나 GRC 등 홈페이지에서 다운로드 받을 수 있다. 

 

뒤에서 설명하겠지만, 레퍼런스와 똑같은 것은 굳이 콜을 할 필요가 없어서 넘어갈 수 도 있다.

스킵하는 것의 장점은 데이터의 양이 확연히 줄어든다는것!

레퍼런스와 비교하면서 다른 것은 SNP, 레퍼런스에는 없는데 나한텐 있는 것은 insertion, 나한테 없는데 레퍼런스에 없는 것은 deletion 즉 INDEL 등 으로 처리할 수 있다.

 

 

컴퓨터 용어로, 앞단의 아웃풋이 뒷단의 인풋이 되어 쭉 이어져 있는 어떤 워크플로우를 파이프라인이라고 하는데,

시퀀싱을 하는것은 이와 같다.

그리하여 앞단의 것이 꼬이면 뒷단의 것도 꼬이면서 연쇄작용이 일어난다.

과정과정 중에 log 파일을 남기고 QC를 하여 관리를 한다.

 

단계는 크게 3단계로 나눠볼수있다.

Data processing은 맵핑단계이다. 여기서는 처음 나온 리드고 맵핑을 아직 하지 않았기 때문에 unmapped reads라고 한다.

이 unmapped read가 들어있는 파일이 FASTQ 파일이다. 이게 시작점.

이후 이 파일을 가지고 레퍼런스에 맵핑을 한다.

맵핑이 되면 mapped read가 되고 이 파일은 SAM 또는 BAM 파일.

 

이 BAM 파일을 가지고 QC를 또 해야한다.

이렇게 되면 분석 준비가 된 리드가 된다.

 

이후 variant discovery에서 variant 콜링을 한다.

위 그림을 보면 판이 여러장 겹쳐져 있는데, 한 판당 한 사람이라고 보면 된다.

그래서 레퍼런스와 비교를 하는데, 이런 경우가 있다.

첫번째 사람은 특정 포지션 a에서 레퍼런스와 달라 variant로 되었고, 두번째 사람은 포지션 a에서 레퍼런스와 다르지 않아 콜링을 할 필요가 없지만, variant calling을 할때는 전체 표본들을 기준으로 보기때문에 콜링을 해야한다.

(내가 가진 population들과 비교를 해야하기 때문에)

따라서 콜링의 아웃풋인 GVCF는 전체 표본들의 variant의 합집합이라고 보면 된다.

 

VCF는 Variant Calling Format인데 GVCF는 Genome를 붙힌것이다.

 

한사람에서 variant는 400만개 정도가 나오는데 모든 사람들의 variant를 합집합으로 하면 1억개가 나올 수 도 있다.

이 모든 변이들을 다 분석할 것은 아니고, 시퀀싱의 목적은 rare vaiant를 찾는 것이 목적이기 때문에, common한 것은 버리고 많이 알려진 변이도 버리고 하는 등 variant 필터링을 거친다.

이런 과정을 거치면 effect가 크고 rare한 변이들만 남는다.

 

머리나 혈액 등 어디에서 DNA를 뽑아도 똑같은 정보인 germline DNA에 연구는 이렇게 한다.

 

그럼 cancer와 같은 somatic은 inherient 연구랑 조금 다르다.

맵핑하는 단계까지는 동일하지만, variant 콜링하는 단계부터는 tumor와 normal 페어끼리 비교하는 등 아예 방식이 다르다. 이 정도만 알아두면 된다.

 

 

NGS 연구를 하는 랩이라고 하면 시퀀싱만 하는 랩이 있는데 이를 wet lab이라고 하고, 데이터만 가지고 분석을 하는 랩은 드라이 랩이라고 한다.

 

실제로 데이터를 받는 형식은 FASTQ라는 데이터 이다.

이는 숏리드들의 시퀀스가 나열되어 있는 파일들이다.

원래는 FAST 파일인데 여기서 Quality 정보가 추가되었다고 해서 FASTQ 파일이 되었다.

 

일반적으로 paired-end 시퀀싱을 하기때문에 한 사람에 대해서 파일을 두개를 받는다.

포워드, 리버스 시퀀싱 각각 두개를 하기 때문이다.

paired-end 시퀀싱은 하나의 리드에 대해서 가운데에 일정한 간격을 두고, 포워드 한번, 리버스 한번 이렇게 시퀀싱을 수행하는 것이다. 이것의 장점은 나중에 통합을 할때 두 시퀀스 사이의 간격을 미리 알고있기 때문에 서로 짝꿍 시퀀스가 되어 맞추기 수월해진다. 예를 들면 500bp의 리드가 있을때 포워드로 150bp, 리버스로 150bp하고 300bp는 간격을 두는 것이다.

아무튼 그래서 R1, R2 파일 두개를 받고, 위 그림은 두명에 사람에 대한 FASTQ 파일인 것이다.

 

보통 1x로 시퀀싱을 하면 3Gb 용량의 데이터가 나오지만, 100x로 시퀀싱을 하면 300Gb로 용량이 매우 커진다.

하지만 퀄리티 정보 등 FASTQ 파일에는 염기서열 정보만 있는 것은 아니기때문에 보통 x3해서 한 사람당 1테라정도의 용량이 된다.

WES과 같은 엑손 시퀀싱은 100x로 하면 염기서열만 5기가바이트 정도 나온다.

 

로우 데이터를 열어보면, 한 리드에 대한 정보가 위 그림처럼 생겼다.

복잡하게 생겼지만, 단순하게 보자면,

첫번째 줄은 해당 리드의 이름이다.

두번째 줄은 리드의 length만큼 염기서열 정보가 있고,

세번째 줄은 단순히 한칸 띄운 것.

네번째 줄은 각 염기서열의 퀄리티 정보이다. 한 염기서열에 대응하는 한 칸으로 표현하기 위해서 아스키(ascii) 코드로 되어있다.

 

Q score는 로그스케일로 계산이 되기 때문에, Q10이라면 에러율이 10%라는 의미이고, Q30이면 에러율이 0.1%라는 의미이다. Q score의 스탠다드는 Q30이다. 이 Q score를 Phred score라고도 한다.

 

실제 리드를 열어보면 위와 같이 되어있다. FASTQ 파일은 리드들이 그림의 좌측상단처럼 무질서하게 엎어져있는 상태이다.

 

시퀀스에 영향을 미치는 에러 원인들은 시퀀서 자체에서 발생할 수 도 있고, 단계단계마다 매우 원인은 많을 수 있지만, 우리가 알고 있는 정보는 Q스코어밖에 없기 때문에, 이를 기준으로 QC를 한다.

 

Q score가 안좋게 나오는 이유는 지난 시간에 설명한 것 처럼 아웃오브싱크가 일어났을때 이다.

염기서열이 합성될때 같은 클러스터 안에서는 동시다발적으로 순차적으로 일어나야하는데, 몇몇의 리드들이 싱크가 안 맞게되면 Q score가 안 좋게 나온다.

 

데이터 QC할때 많이 보는 그림이 위와 같은 그림이다.

x축은 read length이고, y축은 Q score이다.

위 그림에서 왼쪽은 R1, 오른쪽은 R2를 나타낸다고 했을때, 초반에는 Q score가 좋다가 나중에 에러가 나는 모습을 볼 수 있다. 리버스는 포워드보다 더 좋지 않다. R1이 R2보다 대게 퀄리티가 좋다.

 

만약 시퀀싱 회사에 시퀀싱을 맡겨서 받아서 쓰는 입장이라면 QC를 반드시 해야한다.

주는대로 그대로 받아서 쓰게되면 퀄리티에 대한 보장을 할 수 없다.

 

 

이제 FASTQ 파일의 프로세싱은 끝났다. Q score를 가지고 QC를 하는 것이 FASTQ의 과정!

이제 레퍼런스와 맵핑을 하는 과정이 필요하다.

준비물은 레퍼런스 시퀀스와 QC가 끝난 나의 FASTQ 파일, 그리고 맵핑을 도와주는 툴. 이렇게 3가지만 있으면 된다.

 

맵핑을 완료하게 되면 SAM 또는 BAM 파일을 얻을 수 있고, 이를 도와주는 SW는 BWA 또는 samtools가 있다.

 

위 그림처럼 레퍼런스 시퀀스와 나의 리드를 비교하는데, 이게 레퍼런스와 100% 맞을 수 도 있지만 위 그림처럼 몇개는 안 맞을수도 있는데, 이를 어떻게 처리하는지는 맵핑해주는 툴들이 복잡한 알고리즘을 가지고 알아서 맞춰준다.

alignment 또는 mapping을 완료한 파일은 SAM 파일이라고 한다.

SAM 파일을 열어보면 위 그림처럼 생겼다.

FASTQ파일과 비슷하게 생겼지만 조금더 간단하다. FASTQ처럼 두번째줄과 세번째줄은 동일하다.

하지만 레퍼런스 포지션이 들어가있다. FASTQ에는 위치정보가 없다.

위 그림에서는 해당 염기서열이 19번 염색체의 8882171 위치라는 것이다.

CIGAR 정보에서 76M은 76의 리드 길이를 가지고 있는 시퀀스가 M이면 전부다 레퍼런스와 매칭이 되었다는 의미이다.

위 그림에서 우측 상단의 CIGAR strings를 보면 example: "M10 I4 M4 D3 M12" 라는 뜻은 시퀀스의 10개는 매칭이 되었고 그 다음 4개는 insertion이 된것이고(INDEL) D3은 deletion이 된것(INDEL)이라는 뜻이다.

 

SAM 파일에는 레퍼런스와의 비교한 정보도 있기때문에 FASTQ파일보다 용량이 크다.

그래서 압축을 해야하는데 그게 SAM의 바이너리 파일인 BAM 파일이다.

 

BAM 파일에는 각 리드에 대한 검색을 빠르게 처리하기 위해 index를 붙혀넣는다.

그리하여 BAM파일은 bam파일과 bai파일 이렇게 두개가 한 세트이다.

BAM 파일과 SAM 파일은 양방향 전환이 가능하기 때문에 너무 큰 SAM은 용량관계상 지우기도 한다.

 

SAM/BAM 파일을 시각화하기 위한 툴들은 위와 같다.

제일 많이 사용하는 것이 IGV(integrative genome viewer)이다.

 

위 그림이 IGV로 시각화한 결과이다. 왼쪽 그림을 확대하면 오른쪽 그림과 같이 보인다.

시퀀스마다 coverage가 다르다는 것을 볼 수 있다.

만약에 average coverage가 높다고 하더라도, 내가 찾은 variant의 coverage가 낮다면 신뢰할 수가 없게 된다.

저런 부위는 다시 시퀀싱을 해야한다.

한편, 만약 레퍼런스와 다른 variant라면 해당 부위는 컬러로 표시된다.

 

위 그림은 WGS과 WES의 차이를 그림으로 나타냈다.

WGS는 모든 염기서열에 대해 빽빽하지만, WES은 아래 파란 엑손 부위에서만 빽빽하게 보이는 걸 볼 수있다.

왜 엑손이 아닌 인트론 같은 부위에도 붙어있느냐 하면, 시퀀스가 유사한 것이 인트론에도 있기 때문에 몇몇은 거기에 붙어있는 것이다.

 

이게 2nd analysis의 read mapping 단계까지 한것이고, 이제 VCF를 만들어야한다. Variant calling format이라고 한다.

GATK이라는 툴은 germline할때 사용하고, muTect이나 VarScan2 같은 것은 somatic할때 사용한다.

 

위 그림은 레퍼런스와 비교했을때 다양한 상황의 예시이다.

하나의 염기서열이 다르다면 SNP이고, 비어있다면 INDEL이다.

 

cancer와 같은 경우는 DNA가 손상이 되어 다른 염색체에 갖다가 붙는 경우가 있다.

이것이 translacation인데, 이런 뮤테이션으로 cancer가 만들어지는 과정이다.

 

또 사람의 레퍼런스에 붙지않고 버려지는 애들을 박테리아와 같은 것에 붙히면 붙는 경우가 있다.

혈액에는 세균이나 바이러스가 있을 수 있기때문에 소량의 DNA들이 있는 것이다. 

이것이 pathogen이다.

 

그렇다면 variant calling은 어떻게 할까?

앞서 설명한것처럼, 레퍼런스 시퀀스는 haploid type으로 깔아놓는다.

하지만 내가 가지고 있는 시퀀스는 무조건 diploid로 나온다. (엄마꺼 하나, 아빠꺼 하나)

하지만 언페이즈 상황이다. (어떤게 엄마꺼고, 어떤게 아빠꺼고 알수는 없는 상황) (아는 상황은 페이즈)

 

우선 Ideal한 상황부터보자면,

homozygous로 나온다면 0/0으로 표기한다. 근데 여기서 페이즈 상황이라면 /가 아니가 |로 표기가 된다.

heteroztgous로 나온다면 0/1으로 표기한다.

alternative homozygous로 나온다는 1/1로 표기한다.

 

하지만 현실은 Readl한 상황과 같다.

앞서 본것처럼 많은 뎁스로 시퀀싱을 하기 때문에 SAM 파일에서는 REAL한 상황처럼 되어있다.

위 그림처럼 만약 리드 뎁스가 10이면, 확률적으로 접근하여 반반으로 쪼개서 하나는 엄마꺼, 하나는 아빠꺼로 생각한다.

A가 10개가 나왔으면 5/5로 해서 최종 'AA'로 콜을 해야한다.

다음으로 TC가 반복해서 5/5로 나왔다면, 'TC'로 콜을 한다.

 

근데 마지막 상황처럼 C가 8개, G가 2개가 나왔다면 어떻게 해야할까?

8/2로 나온다면 마이너한 G를 에러로 취급하여 'CC'로 콜링한다. 

이런 조건은 연구자가 정하는 것이고, 골드 스탠다드는 있다.

 

한편, germline연구가 아닌 cancer 연구에서는 8/2로 나오는것이 에러가 아니다.

각가의 리드들이 다 다른 셀에서 왔을 수 있기 때문에, 마이너한 것이 우리가 찾고자 하는 리드일 수 있기 때문이다.

 

아무튼 우리가 그래서 얻게되는 파일은 위와같은 VCF 파일이다.

그림이 조금 작지만, 우측하단에 있는 그림에서 ##은 주석이고 실제 데이터는 #부터 시작이다.

##을 살펴보면 각 컬럼이 어떤 의미를 가지고 있고, 컬럼 내의 정보를 어떻게 표기를 했는지 알수있다.

 

간단하게만 살펴보면,

1행은 VCFv4.2 포멧이고, 4행은 레퍼런스 DB는 NCBI36을 사용했고, 5행은 'homo sapiens' 즉, 휴먼이고,

7행부터 나오는 INFO에서 NS는 number of sample을 나타내고, TD는 total depth를 나타내고 등등 확인해보면 된다.

QUAL는 Q score를 나타내는데 옆에 FILTER를 보면 PASS는 기준을 통과했다는 뜻이고, q10은 기준이 Q score 10인데 통과하지 못했다는 뜻이다.

 

NA00001부터 NA00003은 사람(샘플)인데, 1번 샘플은 G|G, 2번은 A|G, 3번은 A/A 인 것이다.

GWAS 돌릴때 이 genoype 정보만 뽑아내면 genotype table로 만들 수 있는 것이다.

 

포멧을 좀더 살펴보면 위와 같다.

 

 

그럼 이제 우리는 레퍼런스와 다른 variant들을 찾아내어 VCF 파일을 가지고 있는 상태이다.

이제 해야할 것은 그 variant에 대한 주석을 달아주고, 필터링을 하는 과정이 필요하다.

 

모든 variant를 다 할수는 없으니, 위와 같은 질문을 던져 좀더 중요한 variant를 따져보아야한다.

 

이러한 것을 도와주는 SW는 위와같이 많다.

내가 저런 SW에 넣은 VCF 파일의 INFO 열에 우측 그림과 같은 정보들을 추가해준다.

조건을 설정할때, Synonymous만 하겠다 또는 intronic만 하겠다 등등의 조건을 설정해주면 다 필터링이 된다.

 

진단의학과에서는 이런 variant정보를 환자에게 보고를 해주어야하는데, 규칙이 필요하다.

Benign부터 pathogenic까지 variant의 영향력을 따져서 리포트한다.

 

보통은 멘델리안 disease(rare disease)에 기반하여 variant들의 영향력을 측정하고, GWAS 결과로 보고된 variant들도 함께 보고하기도 한다.

 

variant filtering은 위처럼, 퀄리티가 안좋으면 버리고, common하면 버리고, 이미 알려져있는 variant면 버리고 등등 

 

어떤 질환하고의 association을 봤을때, 첫째로 봐야할 것은 common한 것이냐 아니냐이다.

이후, population에서는 얼마나 존재하는 것인지. disease cohort에서는 얼마나 관찰되는지.

이러한 다양한 DB들이 잘 구축이 되어있어서 활용하면 된다.

 

이렇게 시퀀스된 나의 데이터를 어떻게 활용하고 의미를 찾을것인가 하는것은 이제 연구자의 역량이다.

단순히 유의미하게 나왔느냐 이런 문제가 아니다.

 

 

 

 

 

끝!!