Sam File

sam file은 시퀀스 정보를 저장하는 파일이다. 근데 fasta 포맷처럼 시퀀스가 일렬로 쭉 나열되 있는 것이 아니라 reference 시퀀스에 align된 시퀀스를 저장한다. 그래서 시퀀스와 함께, 그 시퀀스가 referecne 시퀀스에 매핑된 정보를 CIGAR 포맷으로 나타낸다. 또 bam file로 변환할 수 있는데 bam file은 sam file과 같은 정보를 갖고 있지만 binary 형태로 변환된 데이터이다. sam file은 텍스트 데이터인데 이를 binary로 변환함으로써 데이터를 압축하는 효과를 나타낸다.

이것이 sam 파일의 기본적인 포맷. (https://genome.sph.umich.edu/wiki/SAM)

그리고 각각의 필드에 대한 설명은 https://en.wikipedia.org/wiki/SAM_(file_format)에 나와 있다.

적절한 조합에 대한 다양한 시도가 stacking의 성능을 높일 수 있다.

http://blog.kaggle.com/2016/12/27/a-kagglers-guide-to-model-stacking-in-practice/

k-mer 이란

genomics에서 k-mers란 어떤 시퀀스가 주어졌을 때, 길이가 k 가능한 모든 substring의 집합이다.

위키피디아에 따르면 다음과 같다.

The term k-mer typically refers to all the possible substrings of length k that are contained in a string.

즉, 아래 시퀀스의 경우

ATCGAAGGTCGT

k=4이면 4-kmers는 아래와 같다.

ATCG  TCGA  CGAA  GAAG  AAGG  AGGT  GGTC  GTCG  TCGT

이를 이처럼 Sequence Assembly에도 활용할 수 있다.

Bioinformatics에서는 주로 k-mers가 어떤 시퀀스의 "시그니쳐" 를 나타낸다. 즉, 어떤 시퀀스에서 feature를 뽑을 때, 이 k-mers를 이용한다.

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