# packages ‘Gviz’ and ‘chipseq’
> library (Gviz)
> library (chipseq)
# Load aligned chipseq reads sample ‘cstest’ to workspace.
> data(cstest)
#use names function to see samples included in ‘cstest’
#two sample in this dataset—‘ctcf’ and ‘gfp’
> names(cstest)
[1] "ctcf" "gfp"
# display the data of cstest dataset
> cstest
GRangesList of length 2:
$ctcf
GRanges with 450096 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 more element>
---
seqlengths:
chr1 chr2 ... chrY_random chrUn_random
197195432 181748087 ... 58682461 5900358
# Get ‘ctcf’ data for further analysis
> ctcf <- cstest$ctcf
# display the data of ctcf
> ctcf
GRanges with 450096 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_random chrUn_random
197195432 181748087 ... 58682461 5900358
#calculate coverage of ctcf
> ctcf.cov <- coverage(ctcf)
#get the island of ctcf at chromosome 10
> island <- slice (ctcf.cov$chr10, lower=1)
#display the island range
> range(width(island))
[1] 24 371
# get the largest island in binding region
> largeisland <- island[width(island) == 371,]
#display the largest island information
> largeisland
Views on a 129993255-length Rle subject
views:
start end width
[1] 67475840 67476210 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 the track information of largest island in chromosome 10 including 5000 bp upstream of largest island and 5000 bp downstream of largest island
> chiptrack <- window (ctcf.cov$chr10, start=start(largeisland) - 5000, end=end(largeisland) + 5000)
#display the information of track
> chiptrack
integer-Rle of length 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 11 10
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 the RLE 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="darkblue", fill.histogram="darkblue")
plotting result:
帖子
日志
相册
分享
记录
活动
投票
悬赏
辩论
商品
收藏
勋章