EM 알고리즘을 통한 haplotype frequency 추정

우선 EM 알고리즘의 정의부터 알아보자.

EM 알고리즘은 관측되지 않는 잠재변수(unobserved latent variables)에 의존하는 확률 모델에서 최대우도(maximum likelihood)나 최대사후확률(Maximum A Posteriori)을 갖는 매개변수를 찾는 반복적인 알고리즘이다.

그렇다면 왜 haplotype frequency 추정할 때 EM 알고리즘 쓸 수 있는가? 관측되지 않은 변수가 있기 때문이다. 모수를 추정하고자 하는데 모수 추정에 필요한 값이 관측할 수 없는 잠재변수를 포함하고 있을 때 쓸 수 있다.

아래와 같은 세팅을 생각해보자.

Locus 1 has alleles A; a
Locus 2 has alleles B; b

Genotype counts are observed at the two loci for N individuals in the population.

homozygous for both loci: nAABB; naaBB ; nAAbb; naabb
homozygous for one loci: nAABb; naaBb; nAaBB; nAabb
heterozygous for both loci: nAaBb

근데 문제는 여기서 nAaBb는 2개의 가능한 하플로타입 nAB/ab or nAb/aB 이 있다. (즉, 두 loci 모두 genotype이 heterozygous인 경우 haplotype을 결정할 수 없다는 사실을 알 수 있다.)

우리가 알고 싶은건 pAB; pAb; paB; pab. 바로 이 haplotype frequency들이다.

여기서 pAB를 예로 들어 EM 알고리즘을 적용하는 과정을 살펴보자. 우선 pAB의 MLE는 아래와 같다.

p^AB = nAB/2N (N=관측된 사람 수. 하플로타입은 개인당 2개이므로 2N)

또한 nAB = 2nAB/AB + nAB/Ab + naB/AB + nAB/ab

이 때 문제는 nAB/ab를 모른다. 이를 잠재변수라고 하고, 이는 추정하려는 모수들 (pAB; pAb; paB; pab)과 관련되어 있다. 그래서 모수가 있다면 구할 수 있는 값이다. EM 알고리즘에서는 이를 모수를 적절한 값으로 초기화한 후, 반복적으로 nAB/ab를 구한다.

EM 알고리즘에서는 추정하려는 모수 pAB; pAb; paB; pab 를 적절한 값(ex. 0.25)으로 초기화한다. 그리고 이를 p0AB 처럼 표기한다.

관측한 값

Y = (nAABB; naaBB; nAAbb; naabb; nAABb; naaBb; nAaBB; nAabb; nAaBb)

와 초기화된 모수 p가 주어졌을 때, nAB의 기댓값은 아래와 같이 구할 수 있다.

여기서 E[nAB/ab]는 무엇인가?

왜 이렇게 되냐하면, nAB/ab는 nAaBb는 nAB/ab 또는 nAb/aB 이므로, 두 하플로타입(AB/ab, Ab/aB)이 등장하는 확률 값의 비율을 통해 nAB/ab의 기댓값을 구할 수 있다.

이를 통해 n0AB를 구했으면 pAB의 MLE를 통해 이를 업데이트한다.

p^AB = nAB/2N

이 과정을 반복하면 p^AB가 수렴해가고 이를 통해 추정값을 얻을 수 있다. 다른 모수들에 대해서도 이와 같이 한다.

Chip-seq

Chip-seq은 특정 단백질과 결합하는 시퀀스를 알아내기 위한 방법이다. 유전자 발현을 조절하는 전사인자의 결합위치를 알아내는데 많이 활용된다. 실제로 Chip-seq 데이터를 통해 연구를 해야되는 연구자들이나 분석가들을 위해 정리된 튜토리얼이 필요해서 찾아보다 발견한 것이 아래 링크이다.

http://www.biologie.ens.fr/~mthomas/other/chip-seq-training/

위 링크에서 Downloading 파트에 가면 Chip-seq 퍼블릭 데이터를 다운 받는 방법이 나와있다. 우선, 이 튜토리얼은 FNR이라는 단백질에 관한 Chip-seq 데이터를 대상으로 하고있다. 튜토리얼에서 사용할 데이터는 GSE41195이다. GSE41195란 GEO(Gene expression omnibus)에서 사용되는 식별자(Identifier)이다. 이를 통해 해당 데이터에 접근할 수 있다. GEO홈페이지에서 다시 SRA 식별자 SRX189773를 알아낸다. GEO, SRA는 다 퍼블릭 DB이다. 같은 시퀀스인데 여러 DB에 저장되어있는 것이다. 근데 다시 이 SRA 식별자를 아래 링크 ENA 데이터베이스에서 검색해야지 데이터를 얻을 수 있다.

위 링크에 들어가서 SRX189773를 검색한다. 그러면 위 사진처럼 나오고 여기서 SRX189773을 클릭한다.

그러면 이와같은 화면이 나오는데 여기서 FASTQ files File1을 클릭하면 된다. 저 링크는 ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR576/SRR576933/SRR576933.fastq.gz 여기다. (하지만 이 링크를 클릭하면 아이디/패스워드 인증이 뜬다. 따라서 위 방법대로 차근차근 들어가야함) 결국 Chip-seq 데이터는 fastq 파일로 제공되는 것이었다. sequencing된 read들이 fastq 파일 형태로 나오며, 이를 align 한 후에, 많이 겹치는 부분을 단백질 결합 위치로 예상하는 것이다. 이를 peak라고 한다.

다운로드 받은 fastq 파일을 열어보면 이렇게 생겼다.

dbSNP 관련 정리

유전체 데이터 분석 (김주한 저) 269p 부터 쓰여져 있는 SNP 데이터 분석을 참고하여 개인적으로 실습하고 정리한 자료. SNP 데이터 분석을 위한 데이터베이스 중 dbSNP 데이터베이스 관련하여 정리하였다. 모르는 부분이 많아 개인적으로 정리하는 차원의 포스팅이다.

dbSNP

dbSNP 은 NHGRI와 NCBI가 함께 만들어 운영하는 사이트로써 SNP에 대한 정보를 담고 있다. dbSNP에는 두 가지 고유 아이디가 있는데 ssId와 rsID이다. ssId는 사용자들이 제출한 데이터이고 rsID는 관리자가 교정한 데이터이다.

실습 예

rs380390 이라는 아이디가 붙은 SNP을 검색해보고 관련 유전자를 알아보자.

https://www.ncbi.nlm.nih.gov/snp 에 접속하여 rs380390 검색. 그러면 아래의 정보를 볼 수 있다.

해당 SNP는 1번 염색체의 196731921 위치에 있으며, CFH Gene에 존재한다. GeneView를 통해 CFG Gene 정보를 볼 수 있다. 또한 해당 SNP은 단백질을 코딩하지 않는 유전자 상의 부분인 intro에 생긴 intron variant이다. 다소 이해 안 가는 부분은, rs380390을 클릭하여 들어가면 

MAF

G=0.2450/1227 (1000 Genomes)
G=0.3264/9505 (TOPMED)

라고 된 것을 볼 수 있다. 이는 minor allele frequency/minor allel count 이다. ancestral allele이 G인 것으로 보아 G가 minor allele이 아닌 것 같은데 왜 MAF를 G= 이라고 표시 한 것인지 헷갈린다. G가 아닌 다른 allele을 minor allele로 정의하고 그들의 MAF가 저렇다는 것인가? 즉, 저부분엔 ancestral allel인 G=을 표기해주고 G 이외의 다른 allele이 24.5%가 있다는 것일까? 또 궁금한 것은 HGVS Names부분이다. 저 부분에 저렇게 많은 variant에 대한 notation들이 있는데, 왜 저렇게 많이 있는지 궁금하다. 어떠한 variant에 대한 HGVS notation 은 1개만 있으면 되는 것 아닌가..?

Gene view를 선택하면 해당 SNP가 포함된 gene인 CFH Gene에 있는 SNP들의 정보를 알 수 있다. 해당 snp은 intron variant이므로 in gene region을 클릭하여 검색해야지 나온다. (SNP들이 엄청 많아 페이지 로드하는데 시간이 좀 걸릴수도 있다..) default로 되어있는 cDNA를 클릭하면, 이건 exon 부분만 있는거기 때문에 안나오는듯 하다. 아래처럼 rs380390을 확인할 수 있고, heterozygosity=0.37과 MAF=0.245을 알 수 있다. (근데 heterozygosity와 MAF는 왜 다른거지?..)

dbSNP는 NCBI의 Entrez 시스템에 통합되었다. 따라서 질의문을 이용하여 SNP을 검색할 수 있는데, 예를 들어 이런 경우이다. '5번 염색체와 6번 염색체에 존재하는 모든 frameshift SNP를 찾아라' 같은 경우에 질의문을 통하여 검색하면 이 조건에 해당하는 SNP들을 골라낼 수 있다. 이 질의문에는 물론 문법이 있는데, 이 부분은 포스팅에서는 다루지 않았다. 어렵지는 않아보여, 필요할 때 바로 검색해서 사용하면 될 듯하다!

Introduction to Next Generation Sequencing data analysis

아래 글을 참고하였습니다.

http://genestack-user-tutorials.readthedocs.io/guide/intro-to-ngs.html

NGS 기술은 차세대 시퀀싱 기술이라고 말하며, DNA의 시퀀스 정보를 빠르고 값 싸게 읽어내는 기술이다. NGS 기술이란 광범위한 의미인데, 예를 들어, WGS, RNA-Seq, WES, WGBS, ChIP-Seq 등을 NGS 기술이라고 말한다. WES(Whole Exon Sequencing) 데이터의 예를 들어서 NGS 데이터 분석하는 과정에 대해서 알아보자. WES이란 DNA 상의 Exon 부위의 염기서열을 읽어내는 기술을 말한다.

우선 큰 그림을 살펴보면, 일반적인 NGS 데이터 분석은 raw reads quality control, preprocessing, mapping, post-alignment processing, variant calling, annotation, and prioritization 으로 이루어져 있다.

이제 차례대로 WES 데이터를 처리할 때 어떤 스텝을 거쳐야 하는지 알아보자. 시퀀싱 데이터를 가지고 첫 번째로 해야할일은 바로 시퀀스의 quality를 체크하는 일이다. 이를 quality control(QC)이라고 한다. 일반적으로 read의 길이와 수를 체크하게 되고 contaminating 시퀀스나 낮은 quality의 시퀀스가 있는지 찾아야한다. 이후에는 preprocessing 과정을 거치는데 이는 시퀀스의 퀄리티를 증가시키기 위한 과정이다. QC와 preprocessing 과정은 매우 중요하며, 이 과정이 제대로 되어야지 이후의 분석 결과를 신뢰할 수 있게 된다.

다음 스텝은 mapping이다. 또는 aligning이라고도 불린다. reads를 reference genome이나 reference transcriptome에 정렬하는 것을 뜻한다.  (reference genome은 표준게놈 또는 참조게놈 또는 참조 서열 등으로 불리며, HG19, HG38 등을 예로 들 수 있다.) 이렇게 함으로써 데이터의 뉴클레오타이드 서열을 de novo assembly 없이 연구할 수 있게 된다. 예를 들어 WES 데이터의 경우 WES 데이터의 read를 reference genome에 mapping 한다면, reference와 WES 데이터의 시퀀스 간의 다른 부분(variant)을 알아낼 수 있게 되고, 이 variant의 정확도는 mapping accuracy에 의존적이다.

따라서 그 다음에 해야할 일은 바로 mapping quality를 체크하는 일이다. 데이터의 특정 종류의 bias는 mapping step 이후 나타난다. mapping quality가 만족스럽지 않으면 이를 processing하는 과정을 거쳐야하는데 이를 post-alignment processing 이라고 부른다. 예를 들어, 중복된 mapped read를 제거하는 것을 의미한다. (이는 PCR의 인공물이다.) 이는 매우 중요하며 앞으로의 분석에 큰 영향을 미친다.

시퀀스가 참조 서열에 align 되었다면, 데이터를 experiment-specific한 방법으로 분석을 해야할 필요가 있다. 우리는 여기서 WES 데이터를 참조 서열에 mapping한 후에, variant 분석을 할 것이다. 즉, variant calling과 그 variant가 gene에 미치는 영향 (단백질의 변화, frame shift 등)을 알아볼 것이다. 이 과정에서 시퀀스를 참조 서열과 비교하고, 차이를 확인하여 이 차이가 유전자에 얼마나 큰 영향을 미칠지에 대해 분석하여야한다. 예를 들어, 이것이 snynonymous variant(mRNA가 생성하는 아미노산에 변화가 없는 변이)인 경우에는 이 영향이 미미할 것이다. 하지만 그 variant가 사이즈가 큰 deletion 이라면, 해당 시퀀스를 포함하는 유전자에 큰 영향을 줄 것이라고 예측해볼 수 있다. 분석과는 별개로 데이터를 visualization 해볼 수도 있다. mapped read 데이터를 visualising 함에 있어 가장 표준적인 tool 중 하나는 Genome Browser이다.

ClinVar DB 를 통한 질병 연관 변이 찾기

심근증(cardiomyopathy)과 관련있는 gene은 무엇이고, 그 gene의 어떤 위치의 variant가 pathogenic(질병의 원인이 될 수 있는)한지 알아보고자 한다. 이를 알아보기 위해서는 genotype-phenotype 데이터베이스를 이용하면되는데 이번 포스팅에서는 NCBI의 ClinVar 이라는 데이터베이스를 이용해보도록 하자. 우선 NCBI 홈페이지에 들어가 cardiomyopathy를 검색해본다.

그러면 아래와 같은 화면이 나오는데 ClinVar 을 클릭하여 들어간다.

이곳에서 왼쪽 메뉴바에서 Pathogenic을 클릭한다. 그러면 cardiomyopathy에 pathogenic한 gene의 variant를 찾을 수 있다.

또 많은 연구자들에 의해 보고된 pathgenic한 variant를 알아보고 싶다면, 아래와 같이 Multiple submitters를 체크한 후에 살펴보면 된다.

또 그 variant의 인해 분자적인 특성이 어떻게 변하는지 궁금할 수 있다. 이 때는 Molecular consequence에서 관심있는 분자적 특성에 체크를 하면된다. 예를 들어, missense(variant로 인해 mRNA의 코돈이 변해서 정상적인 경우와 다른 아미노산을 생성하는 분자적 변이)를 클릭하면 아래와 같은 화면을 볼 수 있다. 아래와 같이 다양한 pathogenic한 variant와 해당 gene을 알 수 있다. 해당 gene의 variant로 인해 missense가 일어나고 이로 인해 cardiomyopathy 질병을 일으키는 것이다.

크게보면 아래와 같이 variant를 나타내고 있다. 이렇게 variant를 표기하는 방법을 HGVS notation이라고 한다.

이를 통해 검색할 수도 있는데 예를 들어,

NM_213653.3(HFE2):c.959G>T (p.Gly320Val)

를 ClinVar 검색창에 검색해보자.

그러면 위와 같은 화면을 볼 수 있다. Links를 보면 해당 Variant가 UniProtKB, OMIM, dbSNP, 100 Genome Browser에 등록되어 있다는 것을 알 수 있다. 이를 통해 이 variant가 Multiple Submitters에 submit 되었다는 것을 알 수 있다. (Number of Submissions에 4가 표시되어 있다.)

ClinVar 데이터 다운로드 받기

https://www.ncbi.nlm.nih.gov/clinvar/에 접속 후 FTP => Go to the FTP site

여기서 vcf_GRCh37, vcf_GRCh38이 variant를 vcf 파일 형태로 보여주는 것이고, xml은 xml 파일 형태로 보여주는 것이다.

들어가면 다양한 버전이 존재하는데 가장 위에 있는 것을 받으면 된다. 

CIGAR Format

출처 - https://genome.sph.umich.edu/wiki/SAM

/* 2017.8.13 by 3months. */

CIAGR Format은 어떠한 Sequence A가 reference genome에 align되어있을 때 이를 나타내 주는 포맷이다.

예를 들어, 아래와 같은 reference sequence와 Read가 있다고 하자. 이 Read를 reference sequence에 align 해보면

RefPos:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T  G  A  A  C  T  G  A  C  T  A  A  C
Read: ACTAGAATGGCT

이와 같이 align 된다.

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T

이 때 이것을 나타내는 CIGAR Format은

POS: 5
CIGAR: 3M1I3M1D5M

이렇게 되는데 reference의 5번 position부터 시작해서 3개의 Matching, 1개의 insertion, 다시 3개 matching, 1개 deletion, 5개 matching 이라는 뜻이다. ref의 14번 position에 read와 시퀀스가 다르지만, align이 이렇게 되기 때문에 matching으로 쳐준다.

Lastz를 통한 CIGAR Format 생성

lastz_test.fa

lastz_test2.fa

의도적으로 lastz_test fasta 파일의 sequence에다가 4개의 sequence를 insert하고 2개의 sequence를 delete해서 lastz_test2.fa 파일을 생성하였다.

그 다음

lastz lastz_test.fa lastz_test2.fa --format=general:cigar > cigar.txt

를 실행 후, cigar.txt 파일을 열어보면

#cigar
208M4I163M2D138M

위와 같은 결과가 나왔음을 볼 수 있다. 결과를 해석하면 208개 matching 4개 insertion, 163개 matching 2개 deletion, 138개 matchin 이다.

2017.10.8 추가

lastz 사용해서 cigar을 얻을 때 다른 옵션을 주면 substitution도 고려가 가능하다.

lastz_test.fa

lastz_test2.fa

Lastz 실행

/* 2017.8.13 by. 3months */

Lastz를 설치하였으면 간단하게 결과를 뽑아보자.

두 개의 fasta 파일을 폴더에 붙여넣고,

해당 폴더에서 아래 명령어를 입력해보자.

lastz lastz_test.fa lastz_test2.fa

결과

그러면 콘솔창에 이런 결과가 나온다.

또는 아래와 같이 옵션을 주어서 실행하면 maf 파일로 내보낼 수 있다.

lastz lastz_test.fa lastz_test2.fa --format=maf > lastz_result.maf

lastz_result.maf 파일을 열어보면 아래와 같이 되어있다.  (이것이 무엇을 뜻하는지는 다음에 알아보자..)

##maf version=1 scoring=lastz.v1.04.00
# lastz.v1.04.00 --format=maf
#
# hsp_threshold      = 3000
# gapped_threshold   = 3000
# x_drop             = 910
# y_drop             = 9400
# gap_open_penalty   = 400
# gap_extend_penalty = 30
#        A    C    G    T
#   A   91 -114  -31 -123
#   C -114  100 -125  -31
#   G  -31 -125  100 -114
#   T -123  -31 -114   91

a score=47093
s ref 0 511 + 511 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCT--------AGTGCTCCACCCTG----CCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGTAAATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

s ref 0 521 + 521 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCTAGAGAGAGAGTGCTCCACCCTGCCCCCCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGT--ATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

아래 링크에서 lastz 실행과 관련된 다양한 예제를 볼 수 있다.

http://www.bx.psu.edu/miller_lab/dist/README.lastz-1.02.00/README.lastz-1.02.00a.html#examples

DNA 결합 단백질

Regulatory Genomics 관련 논문을 읽던 중 DNA binding protein을 보고 몰라서 공부하였다.

DNA => RNA => Protein이 되는 과정에서 DNA=>RNA로 변하는 transcription 과정 중 DNA가 RNA로 변하는 것의 시작은 전사개시복합체(transcription initiation complex)가 Promoter 부분에 붙으면서 시작된다. 이 복합체 단백질 중 하나인 RNA중합효소2가 유전자를 전사헤 pre-mRNA를 만든다. 그 다음에 RNA Processing(5'캡 생성, 폴리A 꼬리 생성, 인트론 splicing)을 통해 mature mRNA가 생성된다.

이 과정에서 RNA 중합효소가 전사를 개시하려면 전사인자(transcription factor)라는 단백질이 필요하다. 이들은 보편전사인자(general transcription factor)라고도 불린다. 하지만 이것만 가지고는 전사개시속도가 매우 낮으며 특수전사인자라는 조절 단백질의 작용이 필요하다.

promoter 주변에는 DNA 시퀀스에는 근거리 조절요소가 있다. 그리고 유전자에 멀리 떨어져 있는 원거리 조절요소인 인핸서가 있다. (한 개의 인핸서는 하나의 유전자에만 작용하며, 한 유전자는 여러개의 인핸서를 갖을 수 있다. 인핸서에 붙는 단백질들은 특수전사인자라고 불리며 전사인자 종류중 하나라고 볼 수 있다 )

진핵생물에서 유전자 발현율은 인핸서의 조절요소에 활성자(activator) 또는 억제자(repressor)가 결합함에 따라 증가되거나 감소된다. 결합되는 부분은 DNA Binding site라고 불리는듯하다. 즉 인핸서의 조절 요소에 특정한 단백질이 붙어서 그 인핸서의 표적 유전자의 발현을 조절하는 것이다.  활성자 단백질은 RNA 중합효소를 promoter로 recruting한다. 그래서 전사속도가 빨라진다. 반면 억제자 단백질은 RNA 중합효소가 이동하지 못하게 하여 전사가 될수 없게 만든다.

이렇게 gene expression 조절을 수행하는 여러가지 형태의 단백질(RNA 중합효소, 활성자, 억제자 단백질)을 DNA 결합 단백질 (DNA Binding protein)이라고 하며 이들은 sequence specific하다. 즉, DNA Binding protein을 코딩하는 시퀀스에 문제가 있을 경우 해당 protein의 표적 유전자의 발현에 이상이 생길 수 있다. 그래서 DNA Binding protein의 sequence specificity를 안다면 표적 유전자를 알 수 있게 되고, DNA Binding protein 시퀀스의 disease causing variant를 파악하는데 도움이 된다.

이러한 Transcription factor 등과 같은 단백질에 결합하는 DNA 시퀀스를 알아내는 방법으로 CHIP-seq 방법이 있다.

참고 - 캠밸 생명과학

유전체학을 위한 신경망 모델 (Neural Networks for Genomics)

원본글

http://blog.qure.ai/notes/neural-networks-for-genomics

위 글을 번역한 포스팅입니다.

Introduction

딥러닝을 처음 유전체학에 적용하려고 한다면, 최신으로 이루어지고 있는 연구(state-of-the-art)가 무엇인지, 연구자들이 딥러닝을 통해 어떠한 문제를 해결하고자 하고, 이를 어떻게 접근하고 있는지 알아야할 필요가 있다. 이 포스팅은 유전체학에서 유명한 네트워크 구조를 소개하고, 이 네트워크를을 트레이닝할 때 쓰이는 데이터의 타입, 그리고 모델을 통해 최종적으로 예측한 결과가 어떻게 나오는지를 요약한다.

시퀀싱 기술의 발달, 그리고 1000-genome project, encode, geo와 같은 public dataset의 범람에도 불구하고 아직까지 genotype-phenotype 구분, 그리고 시퀀스를 통해 질병을 예측하는 것을 완벽하게 해내지 못하고 있다. Bredan Frey의 이 talk는 딥러닝과 유전체학의 컨텍스트상에서 왜 더 많은 genome을 시퀀싱 하는 것이 반드시 해답이 아니라는 것을 설명하고 있다. genome은 매우 복잡하고, 서로 상호작용하는 수많은 정보 레이어를 갖고 있다. 현재의 대부분의 접근법은 시퀀스를 통해 직접적으로 표현형을 추측하는 것이 아니라, 시퀀스의 일부를 해석하는 시스템을 구축하는 것이다. 아래에 딥러닝이 유전체학이나 전사체학에 적용된 예를 소개한다.

얕은 fully connected networks 를 이용한 초기 연구

초기 연구중 몇몇은 주성분분석을 통해 차원을 축소한 gene expression 데이터를 single layer fully connected network에 적용한 것이다. 이러한 초기 연구는 gene expression을 통해 tumor type을 구분하거나, tumor grade를 예측하거나, 환자의 생존 여부를 예측하는데에 사용되었다. 이러한 연구의 발달은 가장 예측력이 높은 gene의 subset을 찾거나, gene의 signature를 찾는 연구의 발달을 가져왔다. 이와 비슷하게 인공신경망은 microRNA의 발현 패턴을 통해 대장암의 등급을 예측하거나, 암 상태를 예측할 수 있는 microRNA를 찾아내는 연구에 사용되었다.

Autoencoder를 이용한 차원 축소와 특징값 추출

Autoencoder는 유전체학에서 gene expression 데이터로부터 feature space를 줄이거나 유용한 feature를 찾아내는데 사용되어왔다. 한 예는 autoencoder를 사용해 labling 되지 않은 gene expression data로부터 정교한 feature representation을 구현하고, 이를 통해 tumor를 예측하는 classifier를 구현한 이 페이퍼이다. 더 최근 연구는 denoising autoencoder를 이용해 유방암 gene expression data에서 feature extraction을 한 연구이다. 다음은 ADAGE(paper and repository)인데, 박테리아 gene expression 데이터에 비슷한 접근을 한 연구이다. autoencoder를 통해 추출한 유용한 feature들은 암 분류, 에스트로겐 수용체의 상태, 환자의 생존 예측에 활용될 수 있다. 그러면 하나의 autoencoder로 부터 밝혀진 gene expression profile들을 통해 어떠한 조직에도 일반적으로 적용할 수 있는 feature를 밝혀낼 수 있을까? gene expression 데이터 말고 DNA sequence에 autoencoder를 적용하여 유용한 feature를 추출할 수 있을까?

Deep learning을 통해 gene expression, transcript expression level을 예측하기

 유전자의 발현은 공동으로 이루어지는 경우가 많기 때문에, 서로 다른 유전자의 발현 정도는 상관관계가 매우 높다. 이는 특정 subset의 유전자 발현 데이터만 가지고도 나머지 gene의 발현 정도를 예측해볼 수 있다는 것을 뜻한다. 이는 잠재적으로 유전자 발현 profiling에서 비용과 복잡성을 줄일 수 있다. D-GEX라고 불리는 방법이 개발되었는데 이는 퍼블릭 데이터인 CMAP을 통해 트레이닝된 심층신경망 모델로써, 1000개의 gene expression 데이터를 통해 나머지 gene 의 expression 정도를 예측한다. 최근의 topcoder challenge도 이와 비슷한 문제가 올라와있다. 이와 비슷하지만 더 hard한 태스크는 exon이나 transcript의 expression level을 DNA sequence data로 부터 예측하는 것이다. expression level은 sequence로만 결정되는 것이아니라, cellular context 에도 의존한다. ‘Deep learning of the tissue-regulated splicing code’ 라는 타이틀의 연구는 exon 주변의 dna sequence를 통해 그 exon이 splice in 되는 확률 (PSI)를 예측한다. 직접 손으로 만든 genomics feature들(이 genomics feature들은 특정한 mouse의 cell type에서 splicing pattern을 예측할 수 있는 것들이다.) 이 모델을 트레이닝하는데 사용되었다. autoencoder를 통해 이 feature의 차원을 축소하고, cell type을 나타내는 추가적인 input을 통해 multilayer fully connected network를 트레이닝하였다. 이 방법을 통해 저자들은 tool을 개발하고 검증하였는데, 이것을 통해 SNV가 splicing애 미치는 effect를 scoring할 수 있다.

Epigenomics를 위한 Convolutional Networks

최근 연구의 많은 부분이 바로 전사인자 결합 위치나, 인핸서 지역, 시퀀스를 통한 크로마틴 accessibility 예측과 같은 epigenomic의 문제에 CNN을 이용하는 것이다. DeepBind는 주어진 sequence를 통해 DNA나 RNA 결한 proteins의 specificity를 예측할 수 있는 방법이다. 이것을 하기 위해 많은 수의 protein binding micro-array chip-seq 데이터를 통해 CNN이 training 되었다. convolutional stage에는 one-hot encoding된 시퀀스에서 'motif detector' matiireces의 set을 scan한다. learning된 filter는 유전쳏학에서 DNA seuqnece motif를 설명하고 이해하는데 사용된다. 네트워크는 알려진 motif와 알려지지 않은 motif까지도 learning할 수 있다. DeepMotif는 비슷한 계열의 더 깊은 모델로써, 시퀀스에 따라 transcription factor binding(yes or no)를 분류하고, motif를 추출하는데 더 중점을 둔 모델이다. DeepSea는 ENCODE의 epigenomics 데이터로부터 트레이닝되었고, 단일염기의 mutation이 transcription factor binding과 DNAse sensitivity에 미치는 영향을 예측한다. 이 모델의 주요 기능은 다양한 규모의 시퀀스 기능을 학습하는 계층적 아키텍처의 사용, 광범위한 시퀀스 컨텍스트 스캔 기능 및 예측 기능을 공유하는 다양한 염색질 요소의 멀티 태스크 공동 학습이다. 다른 epigenomics에 사용되는 CNN의 예는 Basset인데, DNA 시퀀스로부터 chromatin accessibility code를 예측한다. DeepCpG는 DNA metylation을 모델링하였고, DEEP은 enhancer나 transcription factor binding으로 인해 gene의 transcription이 증가하는 DNA의 region을 예측하는 ensemble framework이다. 또한 chip-seq 실험에서 noise를 감소시키거나 signal을 증폭하기 위한 CNN-based 방법도 있다.

언어처리로부터 영감을 얻은 모델들

언어처리와 genome해석은 어느정도 비슷한 면이있다. 이는 언어처리에 사용되는 방법이 genomics를 이해하는데 유용할 수 있다는 것을 시사한다.

Recurent Neural Networks

RNN은 긴 범위를 갖는 DNA 시퀀스의 상호작용을 발견하기 위해 사용되어왔다. DanQ는 CNN과 bi-directional LSTM(long-short term memory)의 하이브리드로써, convolution layer는 regulatory motif를 발견하고, recurrent layer는 regulatory 'grammer'를 이해하기 위해 motif들 간의 long-term dependency를 발견하는데 사용된다. 다른 예는 DNA-level splice junction prediction인데 RNN이 splice junction이나 intro과 exon의 경계를 발견하기위해 training이 된다. 또한 바이러스 gemoe에서 proten coding region을 발견하는 것도 있다.

Word embedding이나 word2vec 같은 모델

‘Gene2vec: Neural word embeddings of genetic data’ 라는 타이틀을 갖고있는 한 프로젝트는 구글의 오리지날 word2vec을 genome 시퀀스로 구현한 것이다. (genome을 27bp의 'word'로 splitting한다.) 실제로 이것이 인간 genom에 실용적으로 적용될 수 있는지는 흥미로운 주제일 것이다. 비슷한 word2vec 스타일의 모델은 gene expression data로도 트레이닝이 가능하다.

기타 링크

dna binding protein motif

http://2013.igem.org/wiki/index.php?title=Team:XMU_Software/Project/promoter&oldid=359507

Regulatory Genomics 개념

regulatory genomics 관련 자료를 찾아보다 아래 링크를 찾아 번역하여 정리하였습니다.

https://biology.stackexchange.com/questions/17810/what-does-regulatory-genomics-mean

Genome은 크게 두 파트로 나뉘어진다.

1. 단백질을 코딩하는 파트 (coding region)

2. 단백질을 코딩하지 않는 파트 (non-coding region)

이 때 단백질을 코딩하지 않는 두 번째 파트는 또 다시 두 개의 클래스로 나뉘어진다.

2-1. transcribe가 되어 생물학적인 활동을 하는 파트 (long noncoding RNA, miRNA, competing endogenous RNA 등)

2-2. transcribe가 안되는 부분 (unscribed non-coding region)

이 중에서 2-2번 unscribed noncoding region은 아마도 1번파트와 2-1번 파트에 regulation 작용을 할 것으로 생각된다. (transcription factor나 transcriptional coactivators/corepressors에 binding 함으로써)

※ transcriptional coactivators/corepressors는 gene의 시작부분 (promoter) 또는 멀리 떨어진 부분(enhancer)에서 해당 gene이 expression이 될지 안될지를 조절한다. enhancer은 gene과 멀리 떨어져있지만 3차원 상에서는 붙어있을 수 있기 때문에 gene expression을 조절할 수 있다. (genome은 세포 안에서 접혀져 있기 때문이다.)

Regulatory Genomics는 바로 genomics "features"라고도 불리는 이 unscribed non-coding region에 대한 연구이다. 어떻게 그들이 gene regulation을 하는지 알아내는 것이다. 이 분야의 연구의 대표적인 예는 ENCODE project의 파트로써 publish된 페이퍼들이다.

ENCODE 프로젝트 링크

http://www.nature.com/encode/#/threads

ENCODE 프로젝트 소개 (홈페이지 내용 번역)

ENCODE(Encyclopedia of DNA Elements) 프로젝트는 National Human Genome Research Institute가 후원하는 프로젝트로써 genome 상에서의 transcription, transcription factor association, chromatin structure, histone modification을 밝혀내기 위한 프로젝트이다. 이러한 genome 상에서의 기능적인 요소들을 식별함으로써 현재 인간 genome의 80%의 부분이 최소 1개 이상의 biochemical function 을 한다는 것 알아내었다. 이러한 functional annotation에 관한 광활한 자원들은 genome과 유전자의 regulation과 oragnization에 있어 새로운 통찰력을 제공하고 있다.

ENCODE 관련 한글로 정리된 블로그

http://blog.daum.net/kimuks/7535222

https://madscientist.wordpress.com/2012/09/18/%EC%98%A4%EB%8A%98%EC%9D%98-%EB%85%BC%EC%9D%BD%EB%82%A8-%EC%97%94%EC%BD%94%EB%93%9C%EB%A5%BC-%EB%94%94%EC%BD%94%EB%93%9C%ED%95%98%EA%B8%B0/

Regulatory Genomics 및 Epigenomics 관련 강의 자료

https://simons.berkeley.edu/talks/regulatory-genomics-epigenomics