利用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: