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.bedrGADEM을 설치하면 위 예제 데이터가 패키지 안에 기본으로 들어있다. 이 예제 데이터를 통해 실습을 해보도록 한다.
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
'Domains > Bioinformatics' 카테고리의 다른 글
| Samtools, Bcftools 설치 방법 (1) | 2018.05.14 |
|---|---|
| GATK4 설치 (0) | 2018.02.20 |
| Variant calling 이란? (0) | 2017.11.20 |
| IGV Tutorial [BAM File] (0) | 2017.10.16 |
| Sam File (0) | 2017.10.16 |
