Chip-seq 의 기초 이해 

Genome 상의 단백질을 코딩하는 부분이 아닌 지역 (non-coding region) 에서 기능적인 부분을 찾기 위한 노력들이 이루어지고 있습니다. genome 연구 초창기에는 non-coding region 이 junk DNA 즉, 아무런 기능을 하지 않는다고 잘못 알려져 있었던 적도 있지만, non-coding region 에 위치한 다양한 기능 부위 (functional elements)들은 유전자 발현에 중대한 영향을 미치며, 이것이 인간의 복잡성을 결정하는 것으로 이해되고 있습니다. 이런 의미에서 non-coding region 을 이해하는 것은 중요하며, 이를 위해 다양한 실험 데이터가 모이고 있습니다. 그 중 하나가 바로 chip-seq 데이터입니다. 

Chip-seq의 핵심

• Chip-Seq = Chip + Next Generation Sequencing 
• Chip = Chromatin + Immunoprecipitation 

Immunoprecipitation (면역침강반응) 어떠한 sample에서 '특정 물질' (target) 을 찾기 위해 그 특정물질의 항체 (antibody) 를 이용하는 것입니다. Chip은 Chromatin (염색질) 에 Immunoprecipitation 을 하는 것입니다. 이것이 Chromatin 에 적용되는 경우, antibody 를 이용해 DNA 에 특정 위치에 결합한 단백질 (ex. transcription factor) 을 찾을 수 있습니다. 이 부분을 침강시킨 후, NGS 기술을 이용해 시퀀싱하는 것이 Chip-seq 입니다. 

Chip-seq 의 종류 

ChIP–seq and beyond: new and improved methodologies to detect and characterize protein–DNA interactions (Nature review genetics, 2012)

• Transcription factor chip-seq 
• Histone modification chip-seq 

Chip-seq은 크게 transcription binding site (tfbs) 를 찾기 위한 tf chip-seq 과 histone modification site 를 찾기위한 hm chip-seq 으로 나뉩니다. 두 경우 실험적으로 약간의 차이가 납니다 (위 그림 참). 또한 sequencing 이후에, 관찰되는 결과도 상당히 다른데, tf chip-seq 의 결과 나타나는 peak는 수십에서 최대 수백 서열 정도 (narrow)이지만, hm chip-seq 의 경우, 많게는 100만 서열 단위까지 길게 나타납니다 (broad peak).

Chip-seq control 

Chip-seq의 목적 중 하나는 genome 상의 특정 기능을 수행하는 부분: tfbs 또는 hm 을 찾기 위한 것입니다. 이를 chippeak finding 라 하는데, 이를 위해 control 데이터가 필요합니다. control 데이터가 필요한 이유는 간단히 말해 Chip-seq 에서 발생하는 noise 때문인데, noise 는 실험적 noise 와, 실제 생물학적 noise 로 나누어볼 수 있을 것 같습니다. 실험적 noise 의 경우, 해당 target 이 아닌 genome  상의 다른 부분이 sequencing 이 된 것입니다. 약 80~90 % 정도가 target 이 아닌데도 sequencing 이 된 부분 입니다. 이를 background noise 라고도 부릅니다. 생물학적인 noise 의 경우, genome 상의 특정 부분 예를 들어 repetitive sequence 가 있는 부분은 실험적인 이유에서가 아니라 그냥 reference genome 에 alignment 가 잘 되기 때문에 peak 처럼 보이기도 합니다. 따라서 이러한 noise 들을 보정하기 위해 control 데이터가 필요합니다. 

control 데이터는 크게 "input"과 "mock" 으로 나뉩니다. 

• input: cross-link와 sonication 은 됐는데, immunoprecipitation 안된 것 
• mock: IgG 라고하는 특별한 항체를 이용해 genome 상의 random 한 부분이 immunoprecipitation 되도록한 것 

이 중에서 일반적으로 input 이 control 로 많이 쓰입니다. 

Peak Finding

Chip-seq 도 NGS 이기 때문에, 최종 결과로 fastq 파일이 나옵니다. 이를 reference genome 에 alignment 를 해서 bam 파일이 나오게 되는데, 이 bam 파일에서 chip-seq 의 한가지 목적 = functional region 찾기를  하기 위한 것 peak finding 입니다. peak 가 있는 부분은 더 많이 sequencing 이 되었다는 것이고, chip-seq 에서 immunoprecipitation 이 많이 된 부분이기 때문에 알고자 하던 부분 (tfbs, hm)일 가능성이 크기 때문입니다. peak finding 을 할 때 여러가지 이슈가 있기 때문에 여러 복잡한 알고리즘들이 많이 쓰입니다. 이러한 알고리즘을 종합해서 구현해 놓은 tool 중 가장 많이 쓰이는 것이 MACS2 라는 Tool 입니다. Chip-seq 의 경우, single-end sequencing 이기 때문에 strand-dependant bimodality 의 문제가 생기는데, MACS 에서는 이를 보정하기 위한 shifting 모델을 사용하고, local dynamic 파라미터를 활용한 포아송 분포 모델을 통해 특정한 크기를 갖는 window로 genome 을 훑으며 통계적으로 유의한 지역 (peak) 을 찾습다. 

참고

http://epigenie.com/wp-content/uploads/2013/02/Getting-Started-with-ChIP-Seq.pdf

SAMtools BCFtools 1.8버전 설치

다운로드 링크

http://www.htslib.org/download/

설치 절차

1. 위 링크에서 우분투 운영체제 wget을 통해 samtools, bcftools tar.gz2 파일을 다운받는다.

2. tar -xvf [파일명] 을 통해 압축을 해제한다.

3. 이 때, 압축 해제는 아무데나 해도된다.

cd samtools-1.x    # and similarly for bcftools and htslib
./configure --prefix=/where/to/install
make
make install

4. 이후 /where/to/install/bin에 samtools binary 파일이 생성된다.

5. /usr/bin 에 심볼릭 링크를 만들어준다.

sudo ln -s /where/to/install/bin /usr/bin/samtools

6. samtools --version  명령어를 통해 잘 설치 되었는지 확인한다.

BCFTools도 이와 마찬가지 방법이다. 폴더 명만 바꿔주면된다.

GATK4 설치

이곳에서 원하는 버전 zip파일 다운로드

https://github.com/broadinstitute/gatk/releases

예를 들어, 4.0.1.1 버전이면 다운로드를 받은 후,

unzip gatk-4.0.1.1 

명령어를 통해 압축을 풀면 설치가 끝난다.

만약에 binary 파일을 다른 곳에 위치시키고 싶다면(예를 들어 /usr/bin), gatk 파일을 카피해서 원하는 폴더 (/usr/bin)으로 옮긴다.

cp gatk /usr/bin

하지만 binary 파일이 다른 파일들과 떨어져있을 경우, 기존 폴더의 경로를 명시적으로 지정해주어야한다. 이는 환경변수를 지정함으로써 할 수 있는데 환경변수 $ GATK_LOCAL_JAR 을 gatk-package-4.0.1.1-local.jar 파일의 위치로 지정하면 된다.

예를 들어, export GATK_LOCAL_JAR=/home/gatk-package-4.0.1.1-local.jar 을 통해 환경변수를 지정하면, binary 파일이 떨어져있어도 실행 가능하다.

Chip-seq 데이터를 통한 binding motif의 분석

Chip-seq 을 통해 대략적인 protein binding site (혹은 histon modification)의 시퀀스를 알고, 이 시퀀스들 중에서 비슷한 시퀀스를 찾아내어 해당 protein에 specific한 시퀀스가 무엇인지를 알아내는 분석이다. 이를 수행하는 bioinformatics Tool은 매우 많으며, [Evaluating tools for transcription factor binding site prediction] 논문에서 좋은 평가를 받았던 rGADEM을 이용하여 binding motif를 찾는 분석을 해보도록 한다. rGADEM은 GADEM 방법을 r로 구현한 것이며, GADEM은 motif를 찾는데 Genetic 알고리즘과 EM 알고리즘을 결합하여 사용한다. 자세한 방법론은 본 포스팅에서 다루지 않는다. 본 포스팅은 이 튜토리얼을 참고하였다.

# rGADEM Tutorial
# reference : https://bioconductor.org/packages/devel/bioc/vignettes/rGADEM/inst/doc/rGADEM.pdf

source("https://bioconductor.org/biocLite.R")
# biocLite("rGADEM", lib="C:/Users/Documents/R/win-library/3.4")
# biocLite("BSgenome.Hsapiens.UCSC.hg19", lib="C:/Users/Documents/R/win-library/3.4")
# biocLite("DNAStringSet", lib="C:/Users/R/win-library/3.4")

# install.packages("RCurl")

library(RCurl)
library(rGADEM)
library(BSgenome.Hsapiens.UCSC.hg19)

rGADEM은 bioLite 패키지에 있으며, 위와 같이 설치할 수 있다. rGADEM, BSgenome..., DNAStringSet 3개의 bioLite 패키지를 설치한 후, 이와 의존성이 있는 RCurl 패키지까지 설치하면 라이브러리 로드가 완료된다.

데이터

Test_100.bed

Test_100.fasta

rGADEM을 설치하면 위 예제 데이터가 패키지 안에 기본으로 들어있다. 이 예제 데이터를 통해 실습을 해보도록 한다.

chipseq 데이터로 binding motif를 분석하기까지의 pipeline은 QC, peak calling 등의 과정을 거쳐야한다. 하지만 이번 분석에서는 peak calling까지 끝난 데이터로 분석을 한다고 가정한다. peak calling이 완료되면 protein binding이 많이 일어난 부분의 시퀀스를 얻어낼 수 있다. 이 데이터는 주로 BED 혹은 FASTA 파일로 주어진다.

1. 데이터가 BED 포맷인 경우

BED <- read.table("C:\\Test_100.bed",header=FALSE,sep="\t")
BED <- data.frame(chr=as.factor(BED[,1]),start=as.numeric(BED[,2]),end=as.numeric(BED[,3]))

BED

rgBED <-IRanges(start=BED[,2],end=BED[,3])
Sequences <- RangedData(rgBED,space=BED[,1])

rgBED
Sequences

> BED

chr     start       end
1  chr1 145550911 145551112
2  chr1 112321064 112321265
3  chr1 120124183 120124384
4  chr1 114560780 114560981
5  chr1   8044002   8044203
6  chr1 154708057 154708258
7  chr1  21537833  21538034
8  chr1 117197889 117198090
9  chr1  22547577  22547778
10 chr1 170982215 170982416
11 chr1 203938925 203939126
12 chr1 198607893 198608094
13 chr1 244838696 244838897
14 chr1  22613354  22613555
15 chr1  36725231  36725432

50 chr1  35024860  35025061

> rgBED
IRanges object with 50 ranges and 0 metadata columns:
           start       end     width
       <integer> <integer> <integer>
   [1] 145550911 145551112       202
   [2] 112321064 112321265       202
   [3] 120124183 120124384       202
   [4] 114560780 114560981       202
   [5]   8044002   8044203       202
   ...       ...       ...       ...
  [46]  24574361  24574562       202
  [47]  46110458  46110659       202
  [48] 114377432 114377633       202
  [49] 152657429 152657630       202
  [50]  35024860  35025061       202

> Sequences
RangedData with 50 rows and 0 value columns across 1 space
       space                 ranges   |
    <factor>              <IRanges>   |
1       chr1 [145550911, 145551112]   |
2       chr1 [112321064, 112321265]   |
3       chr1 [120124183, 120124384]   |
4       chr1 [114560780, 114560981]   |
5       chr1 [  8044002,   8044203]   |
6       chr1 [154708057, 154708258]   |
7       chr1 [ 21537833,  21538034]   |
8       chr1 [117197889, 117198090]   |
9       chr1 [ 22547577,  22547778]   |
...      ...                    ... ...
42      chr1 [221975293, 221975494]   |
43      chr1 [ 59123669,  59123870]   |
44      chr1 [ 76297843,  76298044]   |
45      chr1 [188010599, 188010800]   |
46      chr1 [ 24574361,  24574562]   |
47      chr1 [ 46110458,  46110659]   |
48      chr1 [114377432, 114377633]   |
49      chr1 [152657429, 152657630]   |
50      chr1 [ 35024860,  35025061]   |

2. GADEM 실행

gadem <- rGADEM::GADEM(Sequences,verbose=1,genome=Hsapiens)

interation을 돌며 GADEM 알고리즘이 실행된다.

총 2개의 motif 가 생성되었다.

> length(gadem@motifList)
[1] 2

이 중에서 첫 번째 motif는 아래와 같이 생겼다.

plot(gadem@motifList[[1]])

> gadem@motifList[[1]]@alignList
[[1]]
An object of class "align"
Slot "seq":
[1] "GCCCCGACCTCTTATCTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 100

Slot "pval":
[1] 2.209965e-09

Slot "fastaHeader":
[1] 38


[[2]]
An object of class "align"
Slot "seq":
[1] "ACCCCAACCTCTTATCTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 58

Slot "pval":
[1] 2.365196e-09

Slot "fastaHeader":
[1] 38


[[3]]
An object of class "align"
Slot "seq":
[1] "ACCCCAACCCCTTATTTCTG"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 43552181

Slot "end":
[1] 43552382

Slot "strand":
[1] "+"

Slot "seqID":
[1] 38

Slot "pos":
[1] 121

Slot "pval":
[1] 5.472386e-09

Slot "fastaHeader":
[1] 38

GADEM 알고리즘은 이 motif를 설명할 수 있는 총 26개의 시퀀스를 찾았다. 이 중 3개의 시퀀스가

GCCCCGACCTCTTATCTCTG
ACCCCAACCTCTTATCTCTG
ACCCCAACCCCTTATTTCTG

이다. 이러한 비슷한 시퀀스를 tfbs 라고 하며, tf protein이 이 곳에 달라붙어 regulatory role을 수행하게 된다.

[[26]]
An object of class "align"
Slot "seq":
[1] "CCTTCCTTCCTTTCTTCCTT"

Slot "chr":
[1] "chr1"

Slot "start":
[1] 200254533

Slot "end":
[1] 200254734

Slot "strand":
[1] "-"

Slot "seqID":
[1] 36

Slot "pos":
[1] 193

Slot "pval":
[1] 0.0001878598

Slot "fastaHeader":
[1] 36

위 결과처럼 - strand의 결과도 볼 수 있는데, 이 경우에는 BED파일에서 reverse strand를 보면 된다. 하지만 위 결과는 자동으로 reverse strand를 보여준다. http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:200254533,200254734를 확인하면 위 시퀀스의 reverse인 aaggaagaaaggaaggaagg 를 확인할 수 있다.

3. INPUT이 FASTA 파일인 경우

# 2. Input이 FASTA 파일일 경우
path <- system.file("extdata/Test_100.fasta", package="rGADEM")
FastaFile <- paste(path,sep="")
Sequences <- readDNAStringSet(FastaFile, "fasta")
Sequences

Fasta파일의 경우 위와 같이 전처리 하면 되며, 이후의 절차는 BED 파일과 동일하다.

> Sequences
  A DNAStringSet instance of length 49
     width seq                                                                                                      names               
 [1]   202 CTACAGCTGTTCCTTGTCATCAGCCTGGGGGTGGGTAGTATTTTGATCTTA...GTCCCCATGTGACTGTAGGGTTCCTGAAACCTGGCAGGCCACTCTGCTTG FOXA1 _ 1
 [2]   202 TTATTCTGATGTGGTTTTGCGGTTATACAGTAAGCAGCACTGCTTATGTGG...ACACCAGAGGCCACCAGGACCAGAATGTTTACCAATGTAGGCAGTCACTA FOXA1 _ 10
 [3]   202 AAAGGAGAAACACAGCCAAATAATAAAACAATATCTTCTGTAAGTAAAGAG...CAAAGGTGCAAAGCCCTCTTTCAGATCCATCTCCACCATTTCCCTTCAGG FOXA1 _ 11
 [4]   202 TGTACCCCCCCAATATTTCATGATATTTACATGTTTGCATAGTACTTTCCT...CACAGGCAGAGTCAACATTGGAACTCGGAACACTGAGCCTGGCATTCCAA FOXA1 _ 12
 [5]   202 TTTAAGACTGCCACCTGAAATCAAGTCCAGTGGCTGTTCCTGTCTTCCCGT...TTATTTAACCTTTTGCCATTAGTTTACAAGAAGACGGGTTGAGCGAGTCA FOXA1 _ 13
 ...   ... ...
[45]   202 ATTAGTTCATGCCAGGCAGGGTTTGACCAAAACCTTGGACTCACTCCTGTC...TGCTGCTGCTAAGCCAGTGAGATAGATGGTGGTTTGGAAACACCCTCATG FOXA1 _ 5
[46]   202 CCAGAGCCACCACAGCCAGGCCTCTGCGGCCAGTGTCAACAAGGGGCCAGG...TGGGGTTTTGTTTACACACCTAGGGTCACCTGAAAACACCTGCTGTTTAA FOXA1 _ 6
[47]   202 CAGATCAGAGCCTGGGAGCGGGCCATGTGCAGCTGGATGGGCAGCTGGAGG...GCCCACCCTCCTGCCTTCCCAGCCAAGGAGGGAAGGGAGCTGTGGGAGGA FOXA1 _ 7
[48]   202 GGAATGTGATTTACCCAGATAAATCATCAGCTCAAGGGACTGCTTGGAGAG...TCCAGAGCTAGACGCAGAGGAAGGGGTCAAACATGTCCACATGGAACCTG FOXA1 _ 8
[49]   202 GGCTTTTCCAAGAACAATAGTGTTCTCCTAACAAACATTCGTTCCATCAGC...GGTGTTAGGGACAGAGTCCCAGGTGGCATAAAGTCGGGTTGGTCCTTGGC FOXA1 _ 9

예제 Fasta 파일의 경우 아래와 같은 motif를 찾았고, 총 74개의 sequence가 이 consensus에 부합했다. 즉, FOXA1라는 transcription factor의 binding motif (혹은 transcription factor binding site)의 consensus를 위의 Fasta 파일의 데이터를 통해 예측해볼 수 있다.

plot(gadem@motifList[[1]])

[[72]]
An object of class "align"
Slot "seq":
[1] "CTATTGACTTT"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "-"

Slot "seqID":
[1] 4

Slot "pos":
[1] 93

Slot "pval":
[1] 0.0001517677

Slot "fastaHeader":
[1] 4


[[73]]
An object of class "align"
Slot "seq":
[1] "ATATTTACTTG"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "+"

Slot "seqID":
[1] 42

Slot "pos":
[1] 135

Slot "pval":
[1] 0.0001719873

Slot "fastaHeader":
[1] 42


[[74]]
An object of class "align"
Slot "seq":
[1] "CTATTTGCTGG"

Slot "chr":
[1] "chr"

Slot "start":
[1] 0

Slot "end":
[1] 0

Slot "strand":
[1] "+"

Slot "seqID":
[1] 41

Slot "pos":
[1] 53

Slot "pval":
[1] 0.0001754087

Slot "fastaHeader":
[1] 41

Variant Calling

정의

variant call 은 어떤 개인의 genome혹은 transcriptome에서 reference와 뉴클레오타이드 서열의 차이가 있는지를 결정하는 것이다. 이 과정은 variant frequency와 confidence measure가 고려된다.

유즈케이스

DNA-seq : variants

• 질병과의 유전적 연관성
• 암에서의 돌연변이 알아내기
• 이형접한 인구 빈도 알아내기

RNA-seq : allele-speicifc expression

• Allelic imbalance
• Association with isoform usage (splicing QTLs)
• RNA editing (allele absent from genome)

Chip-seq : allele-specific binding

• protein binding이 allele-specific하게 일어나는 현상

핵심은 variant call을 통해 어떤 locus에 'variant가 있다' 라는 정보를 알 수 있기 때문에 이를 통해 위와 같은 연구를 할 수 있음

Variant Call은 Genotype 보다 더 일반적인 접근이다.

genotype은 어떤 locus의 allele의 집합을 나타낸다. 근데 배수성을 보통 2n으로 가정하고 한다.  또한 genotype은 주로 SNP에 대해서만 고려된다. 하지만 RNA-seq 같은 경우 2n이 아니며 genotype이 없다. 또 암유전체학 등에서는 copy number variation이 중요하기 때문에 variantl call이 genotype보다 일반적인 접근법이라 할 수 있다.

참고

https://www.bioconductor.org/help/course-materials/2014/CSAMA2014/3_Wednesday/lectures/VariantCallingLecture.pdf

파일첨부

VariantCallingLecture.pdf

참고할 사이트들

Best Practices workflow

https://software.broadinstitute.org/gatk/best-practices/bp_3step.php?case=GermShortWGS

IGV Tutorial

IGV(Intergrative Genomics View)는 게놈 데이터를 시각화하고 인터랙티브하게 볼 수 있는 생명정보학 툴이다. 또한 다양한 포맷의 데이터를 로드할 수 있어 매우 편리한 것 같다. array-based 데이터, ngs 데이터, annotation 데이터 모두 로드할 수 있다고 한다.

1. 다운로드

http://www.broadinstitute.org/software/igv/download

실행은 다운받은 폴더에서 igv.bat 을 실행하면 된다.

2. reference fasta 파일 로드

Genomes - Load Genome From File을 통해 reference file을 로드한다. 필자는 hg19 13번 염색체 fasta 파일을 로드하였다.

3. bam 파일 로드

File - Load from File에서 bam 파일을 로드한다. 이 때 bam 파일 인덱스 파일도 필요하다. IGV에서 로드하기 전에 samtools를 통해 bam 파일의 인덱스 파일을 만든다. bam 파일은 reference 에 mapping된 시퀀스 데이터를 나타내주기 때문에 이처럼 reference 파일과 같이 로드해야한다. bam 파일은 ENCODE 데이터 중에서 아무거나 가져와서 로드했다. 용량은 91.6Mb 정도 됐던 것 같다.

둘 다 로드를 완료하면 위 화면과 같이 나온다. 매핑된 부분이 옅은 회색으로 나오는 것을 볼 수 있다. 짙은 회색은 coverage를 나타낸다. bam 파일에는 어떤 reference 데이터에 매핑되어있는지를 이름으로 나타 내는데 만약 이 이름이 다르다면 bam 파일을 로드할 수 없다. 예를 들어, 앞서 로드한 reference fasta 파일의 >chr13 부분을 다르게 변경하면 bam 파일을 로드할 수 없다. 왜냐하면 bam 파일에는 chr13이라는 이름으로 align 되어있기 때문이다. 이 데이터는 ENCODE에서 받았는데 시퀀스가 이어져있지않고 조각조각 흩어져있다. "특정 기능" 을 하는 부분만 시퀀스 된 정보로 볼 수 있다. 저 회색 부분에 마우스를 올려보면 아래와 같은게 뜬다.

이는 그 리드의 기본적인 정보(이름, 길이)와 reference sequence에 매핑된 위치와 CIGAR 정보를 나타내준다.

참고

TreeGenes_IGV_Tutorial.pdf

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)에 나와 있다.

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를 이용한다.

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들을 골라낼 수 있다. 이 질의문에는 물론 문법이 있는데, 이 부분은 포스팅에서는 다루지 않았다. 어렵지는 않아보여, 필요할 때 바로 검색해서 사용하면 될 듯하다!