教你学会ChIP-seq分析 | 第七讲

#我这里直接从GEO里面下载了peaks结果，它们详情如下：wc -l *bed

6768 GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed

3660 GSM1278643_Xu_MUT_rep2_BAF155_MUT.peaks.bed

11022 GSM1278645_Xu_WT_rep1_BAF155.peaks.bed

5260 GSM1278647_Xu_WT_rep2_BAF155.peaks.bed

49458 GSM601398_Ini1HeLa-peaks.bed

24477 GSM601398_Ini1HeLa-peaks-stringent.bed

12725 GSM601399_Brg1HeLa-peaks.bed

12316 GSM601399_Brg1HeLa-peaks-stringent.bed

46412 GSM601400_BAF155HeLa-peaks.bed

37920 GSM601400_BAF155HeLa-peaks-stringent.bed

30136 GSM601401_BAF170HeLa-peaks.bed

25432 GSM601401_BAF170HeLa-peaks-stringent.bed

#PeakID chr start end strand Normalized Tag Count region size findPeaks Score Clonal Fold Change

chr20 52221388 52856380 chr20-8088 41141 +

chr20 45796362 46384917 chr20-5152 31612 +

chr17 59287502 59741943 chr17-2332 29994 +

chr17 59755459 59989069 chr17-667 19943 +

chr20 52993293 53369574 chr20-7059 12642 +

chr1 121482722 121485861 chr1-995 9070 +

chr20 55675229 55855175 chr20-6524 7592 +

chr3 64531319 64762040 chr3-4022 7213 +

chr20 49286444 49384563 chr20-4482 6165 +

library(ChIPpeakAnno)

bed "extdata", "MACS_output.bed", package="ChIPpeakAnno")

gr1 "BED", header=FALSE)

## one can also try import from rtracklayer

library(rtracklayer)

gr1.import import(bed, format="BED")

identical(start(gr1), start(gr1.import))

gr1[1:2]

gr1.import[1:2] #note the name slot is different from gr1

gff "extdata", "GFF_peaks.gff", package="ChIPpeakAnno")

gr2 "GFF", header=FALSE, skip=3)

ol

makeVennDiagram(ol)

##还可以用binOverFeature来根据特定的GRanges对象(通常是TSS)来画分布图

## Distribution of aggregated peak scores or peak numbers around transcript start sites.

data(TSS.human.GRCh37) ## 主要是借助于这个GRanges对象来做注释，也可以用getAnnotation来获取其它GRanges对象来做注释

## featureType ： TSS, miRNA, Exon, 5'UTR, 3'UTR, transcript or Exon plus UTR

peaks=MUT_rep1_peaks

macs.anno AnnotationData=TSS.human.GRCh37,

output="overlapping", maxgap=5000L)

## 得到的macs.anno对象就是已经注释好了的，每个peaks是否在基因上，或者距离基因多远，都是写的清清楚楚

if(require(TxDb.Hsapiens.UCSC.hg19.knownGene)){

aCR

precedence=c("Promoters", "immediateDownstream",

"fiveUTRs", "threeUTRs",

"Exons", "Introns"),

TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)

barplot(aCR$percentage)

}