NGS 데이터로부터 미지의 유전체 크기 (genome size)를 In silico적인 방법으로 예측한다.

Genome size의 다양성 #

(출처: Alberts et al. Molecular Biology of The Cell, 4th edition, 2002, Garland Science, New York.)

생물군에 따라 다양한 유전체사이즈를 보이고 있고 특정 군종에서는 상당히 넓은 범위에 걸쳐있다. 주목할 것은 흔히 생각하듯이 그 생물체의 organismal complexity와 유전체 크기는 큰 상관관계가 없을 수도 있다는 것이다. 예를 들어 human보다 수십 배, 수백 배 더 큰 유전체 크기를 자랑하는 Protist, Plant, Amphibian 종들도 존재하는 것을 확인할 수 있다. 이를 "C value paradox"라 한다. 생물종의 복잡도와 유전자 개수, 그리고 유전체 크기 등의 상관분석도 재미난 연구 topic이 될 수 있을 것 같다.

생물정보학적인 유전체 분석의 제일 첫 스텝은 아무래도 de novo genome assembly가 아닐까 한다.

만약에 어떤 생물종의 유전체 De novo assembly를 진행하고자 할 때 사전지식이 없어 분석하고자 하는 그 생물종의 유전체 크기를 알지 못한 상태에서 진행한다면 어떨까? 염색체가 몇 개인지 전체 사이즈가 얼마인지 모른다면 상당히 답답할 것이다. 당장 시퀀싱을 얼마나 어떻게 해야 할지도 막막할 것이다. 이는 현재 진행되고 있는 모델종이 아닌 수많은 생물종의 유전체 연구 프로젝트에서 맞닥뜨리고 있는 문제이기도 하다.

일반적으로는 실험적인 방법 (Flow cytometry나 PCR 방법 등)을 통해서든 아니면 유전체 크기가 알려진 유사성이 높은 다른 종으로부터 어느 정도 사이즈를 예상한 상태에서 NGS 분석을 진행한다. 하지만 De novo assembly를 통해 만들어낸 draft genome이 맞게 잘 만들어진 것인지 확인하기 위해서라도 대략적인 유전체 크기를 알면 많은 도움이 될 것이다.

다행히도 이런 고민을 먼저 했던 사람들이 있고 NGS 데이터만 가지고 In silico적인 방법으로 유전체 크기를 예측하는 방법이 고안되었다.

첫 번째 방법: WGS/K-mer 분석 #

BGI (Beijing Genome Institute)에서 소개한 방법으로 이미 여러 편의 논문에서 인용되었고 이제는 전 세계적으로 유전체 분석에서 통용되는 기본적인 분석 항목이 되었다.

원리는 간단하다.

예를 들어, 특정 종의 유전체의 60-fold에 해당하는 NGS reads를 가지고 있다면 임의의 한 부분 (가령 25-mer 영역)을 무작위로 추출했을 때 이론적으로 그 부위가 평균 60번 정도는 읽혔을 것이다.

가령 200MB 유전체를 가진 종의 NGS 데이터를 3.4GB 가지고 있다면 이는 17-fold에 해당하는 양이고 임의의 유전체 부위를 추출한다면 평균 17번 발견할 확률이 높다는 것이다.

만약 genomic sequencing library를 유전체 조각의 유실 없이 제대로 잘 만들었고 시퀀싱도 충분한 양 (15-fold 이상)이 됐다면 특정 시퀀스들을 대상으로 depth frequency를 plotting 했을 대부분의 유전체 조각들은 17-fold에서 peak를 그리는 정규분포의 모양을 할 것이라고 예상할 수 있다. 즉, 현재 시퀀싱을 통해 가지고 있는 서열이 전체 유전체의 17배에 해당하는 서열이므로 이를 17로 나눠 원래 유전체 크기를 예측할수 있다.

$$ \textrm{coverage depth} = {\textrm{genome size} \over \textrm{total base num} } = { \textrm{genome size} \over {\textrm{total read num} \times \textrm{avg read length} }} -- 수식 (1) $$

위 수식을 보면, total read num, avg read length는 이미 알고 있으므로 coverage depth를 알면 미지의 genome size를 계산할 수 있게 된다.

여기서 특정 사이즈의 k-mer시퀀스를 살펴보자. NGS data로부터 만들어질 수 있는 전체 k-mer 개수는 다음과 같다.

$$ \textrm{total kmer num} = \textrm{total read num} \times ( \textrm{avg read length} - \textrm{kmer size} + 1) $$

예를 들어 100-mer짜리 Illumina HiSeq read 하나에서 만들어질 수 있는 17-mer (임의의 k-mer를 17로 정했을경우)의 개수는, 100 - 17 + 1 (= 84)개가 될 것이다. 이러한 HiSeq reads가 300백만 개가 있다면, average read length는 100-mer, k-mer size는 17-mer로 상정할 경우,

$$ \textrm{total read num} = 3,000,000 $$

$$ \textrm{avg read length} = 100 $$

$$ \textrm{kmer size} = 17 $$

$$ \textrm{total kmer num} = 3,000,000 \times (100 - 17 + 1) = 252,000,000 $$

이 된다.

이 2억 5천2백만 개의 k-mer들 중에는 동일한 시퀀스들을 가진 것들이 있을 것이고 unique한 k-mer들의 리스트를 작성, 각 k-mer가 몇 번씩 관찰되었는지 그 frequency를 계산한다. 그렇게 해서 depth가 같은 k-mer들이 전체 k-mer pool에서 차지하는 비율을 plotting 한다. x축은 k-mer depth, y축은 전체에서 차지하는 frequency. k-mer frequency plot의 example은 아래 giant panda 데이터를 참조.

그리고 다음의 수식도 성립한다.

$$ \textrm{kmer coverage depth} = \textrm{coverage depth} \times \frac{ \textrm{avg read length} - \textrm{kmer size} + 1}{ \textrm{avg read length}} $$

위에서 작성한 k-mer frequency plot에서 peak을 보이는 depth point가 k-mer_coverage_depth가 되며 아래 수식에 의해 NGS 시퀀스 전체의 genome coverage depth를 구할 수 있게 된다.

$$ \textrm{coverage depth} = \textrm{kmer coverage depth} \times \frac{ \textrm{avg read length}}{ \textrm{avg read length} - \textrm{kmer size} + 1} $$

여기서 구해진 coverage_depth로부터 전체 genome_size를 위 수식 (1)로 부터 구할 수 있다.

$$ \textrm{genome size} = \frac{\textrm{total base num}}{\textrm{coverage depth}} $$

JELLYFISH: K-mer frequency distribution 측정 #

다음으로 NGS data로부터 k-mer frequency plot을 얻는 방법이다. 여러 개의 프로그램이 논문에 나와 있는데 그중에서 JELLYFISH를 소개한다.

(출처: Marcais and Kingsford. A fast lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics (2011) 27:764–770)

JELLYFISH는 free software로 해당웹사이트에서 (http://www.cbcb.umd.edu/software/jellyfish/)에서 무료로 다운로드 받을 수 있다.

다음 그림은 JELLYFISH를 이용해 얻은 k-mer frequency plot을 얻은 일례이다. BGI에서 분석한 Giant panda의 경우인데, 17-mer짜리를 이용해서 frequency plot을 얻은 결과이다.

(출처: Li et al, The sequence and de novo assembly of the giant panda genome. Nat. Genet. (2010) 473:311-317)

이 frequency plot으로부터 위 수식에 따라 계산한 결과 Giant panda의 genome size는 2.46GB로 예측되었고 실제 만들어진 draft genome은 2.25GB로 예측결과와 큰 차이가 없는 것으로 나타났다.

즉, k-mer method는 가지고 있는 NGS data가 반영할 수 있는 genome size를 예측해내는 것이라고 보면 될 것이다. 만약 genomic sequence library가 제대로 만들어지지 않아 genome의 일정 부분만 포함되었다면 그로부터 얻어지는 draft genome도 partial일 수밖에 없을 것이고 k-mer frequency plot을 이용한 방법도 마찬가지로 원래보다 작은 genome size를 예측결과로 내어줄 것이다. 무엇보다 시퀀스 라이브러리를 잘 만들고 NGS data를 제대로 얻는 것이 중요하다고 하겠다.

왜 17-mer인가? #

이론적으로는 어떤 k-mer를 사용하든 유사한 결과를 얻을 것이다. 하지만 제대로 된 결과를 얻기 위한 가정은 최소한 특정 k-mer가 해당 유전체에서 딱 한 번 unique하게 발견되어야 한다. 5-mer 짜리는 확률적으로 1KB당 한 번씩 발견될 수 있고, 10-mer는 1MB당 한 번씩 발견될 수 있습니다.

이렇게 사이즈가 작은 k-mer들은 시퀀스가 유전체에 여러 번 무작위적으로 발견될 수 있기 때문에 올바른 결과를 담보할 수 없다. 16-mer쯤 되면 4.3GB당 한 번 정도 무작위적으로 발견되고 이 정도면 나쁘진 않지만

맨 처음 그림에서 보시다시피 4.3GB 보다 더 큰 사이즈의 유전체도 많이 존재하기 때문에 안심할 수 없을 것이다. 17-mer 정도면 (17GB당 1번) 현재 진행되는 유전체 프로젝트에서는 안전한 결과를 보장할 수 있다고 하겠다. 하지만 17Gb 보다 큰 유전체를 가지고 있는 종들도 많이 있어서 그 경우에는 당연히 17-mer보다 큰 사이즈의 k-mer를 이용해야 할 것이다.

그래프의 해석 #

Jellyfish를 이용한 k-mer density plot이 아래와 같이 그려졌다면 몇가지 유전체의 특징을 예측해 볼수있다. 그래프의 X축은 K-mer로 구성된 서열들의 depth를 나타내며, Y축은 K-mer로 구성된 서열들의 수 (charater의 수)를 의미한다. 이중 몇가지 사례를 들어 살펴보자

유전체 서열내 heterozygosity가 많은 경우 #

유전체 시퀀싱은 haplotype 상태에서 시퀀싱을 진행해야하지만, 보통은 haplotype 샘플링이 매우 어렵기 때문에 diploid 혹은 다배체수의 유전체를 그대로 진행하기도 한다. 이때 dipoid 상태의 유전체상 heterzygosity가 높다면 main peak의 앞쪽으로 불룩한 둔턱이 생길수 있다. 즉, X축으로 보면 main peak보다는 낮은 depth로 일부 서열이 상당부분 존재한다는 의미로 일부 영역에서 heterozygosity가 있음을 알수 있다. 심지어 4배체 6배체의 샘플의 경우 main peak가 4개, 6개로 나타날수 있다.

유전체 서열내 repeat 영역이 많은 경우 #

식물이나 해양 생물의 경우 일반적인 포유류에 비해 유전체내 repeat의 비율이 매우 높다. Repeat 영역은 다른 영역과 다르게 반복되는 서열이 많다. 즉, K-mer로 서열화 했을때 동일한 charater를 갖는 서열이 많아지는 효과가 있다. 동일하게 17배 시퀀싱을 했다 할지라도 repeat영역의 서열은 k-mer로 조각되어졌을때 그의 2배, 4배 이상 나타날 것이다. 따라서 이들은 k-mer density plot에서 X축의 depth상 높은, 먼쪽에 peak로 나타난다. 따라서 그래프의 왼쪽 끝 부분이 볼록하거나 낮은 peak가 형성되어 있다면, 해당 유전체에는 repeat의 비율이 높을수 있음을 인지하여 시퀀싱 전략을 수립해야 한다.

두 번째 방법: WGS/RNA-seq/K-mer 분석 #

두 번째로 소개할 방법은 BMC Genomics (2011) 논문에 소개된 내용이다. 버블콘이라 불리는 바다달팽이 (Conus bullatus)의 전사체를 분석한 논문인데 (Ref. Hu et al. BMC Genomics 2011, 12:60) 기본적인 원리는 위에서 소개한 BGI 방법과 유사하다. 차이점은 Genome Shotgun NGS data만 이용하는 것이 아니고 RNA-seq data도 함께 이용한다는 것이다.

RNA-seq에서 얻은 sequence reads에 genome NGS reads를 매핑하는 것으로 Transcriptome을 레퍼런스로 삼아서 genome sequence reads를 붙이고 각 Transcript당 평균적으로 몇 fold로 붙었는지 계산해보면 그 genome NGS reads가 genome의 몇 fold를 반영하고 있는지를 판단할 수 있다는 원리이다.

이 방법의 장점은 genome NGS coverage가 1.5 fold 정도여도 좋은 예측결과를 보여준다는 것이다. 이 방법을 검증하기 위해 논문에서는 C. elegans와 D. melanogaster의 NGS 데이터로부터 가상의 1.5 fold 짜리 모의데이터를 구성해서 테스트를 해보았는데 다음과 같은 좋은 결과치를 얻을 수 있었다고 한다.

(출처: Hu et al, Characterization of the Conus bullatus genome and its venom-duct transcriptome. BMC Genomics (2011) 12:60)

WU-BLASTN으로 "M = 1 N = -3 Q =3 R = 1 wordmask seg lcmask" 옵션을 사용해서 NGS reads를 RNA-seq 어셈블리로 만들어진 transcript contigs에 매핑하고 각 transcript별 coverage depth를 length 기준으로 계산해서 얻은 다음, 위 그래프에서 보는 것처럼 coverage-depth vs transcript percentage plot을 얻는다. 이때 얻어진 maximum peak coverage가 전체 genomic NGS reads의 coverage가 되며, 수식 (2)를 이용해서 전체 genome size를 예측할 수 있다.

만약 소량의 genome NGS data만 있어서 첫 번째로 소개한 k-mer 방법을 사용하기 어려운 상황에서 RNA-seq data를 같이 가지고 있다면, 써볼 만한 방법이라고 생각된다.

References #

1. Marcais and Kingsford. A fast lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics (2011) 27:764–770
2. Li et al, The sequence and de novo assembly of the giant panda genome. Nat. Genet. (2010) 473:311-317
3. Hu et al, Characterization of the Conus bullatus genome and its venom-duct transcriptome. BMC Genomics (2011) 12:60

Incoming Links #

Related Bioinformaticses (Bioinformatics 0) #

• JELLYFISH

Suggested Pages #

• 0.025 Whole-Genome transcriptome map
• 0.025 FASTQ
• 0.025 Reference genome
• 0.025 microarray
• 0.025 Sequencing
• 0.025 FASTA
• 0.025 Roche 454 GS-FLX
• 0.025 염색체
• 0.025 DNA
• 0.025 시퀀싱장비
• More suggestions...