|||
利用R/Bioconductor 软件包“chipseq” 和“Gviz”对蛋白的结合位点进行分析和显示。
# loading ‘chipseq‘and ‘Gviz‘
> library(chipseq)
> library(Gviz)
# Load alignedchipseq reads sample ‘cstest’ to workspace.
>data(cstest)
#use namesfunction to see samples included in ‘cstest’
#two sample inthis dataset—‘ctcf’ and ‘gfp’
>names(cstest)
[1]"ctcf" "gfp"
# display thedata of cstest dataset
> cstest
GRangesList oflength 2:
$ctcf
GRanges with450096 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr10 [3012936,3012959] +
[2] chr10 [3012941,3012964] +
[3] chr10 [3012944,3012967] +
[4] chr10 [3012955,3012978] +
[5] chr10 [3012963,3012986] +
... ... ... ...
[450092] chr12 [121239376,121239399] -
[450093] chr12 [121245849,121245872] -
[450094] chr12 [121245895,121245918] -
[450095] chr12 [121246344,121246367] -
[450096] chr12 [121253499, 121253522] -
...
<1 moreelement>
---
seqlengths:
chr1 chr2 ... chrY_randomchrUn_random
197195432 181748087 ... 58682461 5900358
# Get ‘ctcf’data for further analysis
> ctcf<- cstest$ctcf
# display thedata of ctcf
> ctcf
GRanges with450096 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr10 [3012936,3012959] +
[2] chr10 [3012941,3012964] +
[3] chr10 [3012944,3012967] +
[4] chr10 [3012955,3012978] +
[5] chr10 [3012963,3012986] +
... ... ... ...
[450092] chr12 [121239376,121239399] -
[450093] chr12 [121245849,121245872] -
[450094] chr12 [121245895,121245918] -
[450095] chr12 [121246344,121246367] -
[450096] chr12 [121253499,121253522] -
---
seqlengths:
chr1 chr2 ... chrY_randomchrUn_random
197195432 181748087 ... 58682461 5900358
#calculate coverage of ctcf
> ctcf.cov<- coverage(ctcf)
#get theisland of ctcf at chromosome 10
> island<- slice (ctcf.cov$chr10, lower=1)
#display theisland range
>range(width(island))
[1] 24371
# get thelargest island in binding region
>largeisland <- island[width(island) == 371,]
#display thelargest island information
>largeisland
Views on a129993255-length Rle subject
views:
start end width
[1] 6747584067476210 371 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 ...]
# get thetrack information of largest island in chromosome 10 including 5000 bp upstreamof largest island and 5000 bp downstream of largest island
> chiptrack<- window (ctcf.cov$chr10, start=start(largeisland) - 5000,end=end(largeisland) + 5000)
#display theinformation of track
> chiptrack
integer-Rle oflength 371 with 99 runs
Lengths: 17 7 5 6 6 12 5 1 11 5 9 7 ... 7 5 23 1 6 5 10 2 1 1110
Values: 1 2 1 2 3 2 1 2 1 2 3 4 ... 3 2 1 2 3 2 3 4 3 2 1
#transform theRLE format data to numeric format data
> chiptrack<- as.numeric (chiptrack)
#prepare DataTrack of this island for track ploting
> genome <- "hg19"
> dtrack<- DataTrack (data=chiptrack, start=(start(largeisland)-5000):(end(largeisland) +5000), end=(start(largeisland)-5000):(end(largeisland) +5000), chromosome="chr10", genome=genome,name="Density", background.title="brown", cex.title=1.2, background.title="white",fontcolor.title="darkblue", col.axis="darkblue",cex.axis=1)
#prepare ideogramTrack of chromosome 10, human hg19 genome
> itrack<- IdeogramTrack (genome=genome, chromosome="chr10")
#prepare GenomeAxisTrack
> atrack<- GenomeAxisTrack ()
# plot the tracks of largest binding domain of ctcf at chromosome 10
> plotTracks(list(itrack,atrack, dtrack), type="histogram",col.histogram="red", fill.histogram="red")
plotting result:
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 12:14
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社