R语言实现CHIP-seq数据分析
ChIP-Seq是将ChIP(Chromatin Immuno precipitation)与二代测序技术相结合的技术，高效地在全基因组范围内检测与组蛋白、转录因子等互作的DNA区域。ChIP也称为结合位点分析法，是研究体内蛋白质与DNA相互作用的有力工具，通常用于修饰组蛋白、转录因子、辅因子以及其他染色质蛋白在染色质上的定位及丰度研究。
我们今天为大家介绍一个R语言中可以实现峰注释以及下游功能分析的包ChIPpeakAnno。这个包的安装目前有两种方式：
一种是利用新出的一个R包进行安装，前提是你的R语言版本>3.5.0：
install.packages("BiocManager")
BiocManager::install("ChIPpeakAnno",version = "3.8")
另一种就是不受版本限制，可以直接bioconductor安装：
source("https://bioconductor.org/biocLite.R")
biocLite("ChIPpeakAnno")
同时我们需要将相关的R包安装：
biocLite("EnsDb.Hsapiens.v75")#注释文件h19
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")#注释文件h19信息
biocLite("org.Hs.eg.db")#注释文件h19
biocLite("reactome.db")#通路分析
biocLite("motifStack")#motif注释
biocLite("seqLogo")#序列Logo
对包的载入部分我们就不赘述了，接下来直接看相关的功能。
下面是数据的导入，主要是需要将原始文件的格式转化为Ranges格式：
bed <- system.file("extdata","MACS_output.bed", package="ChIPpeakAnno")#读取文件所的所在路径
gr1 <- toGRanges(bed,format="BED", header=FALSE) #将bed文件转换为 Ranges格式
gff<-system.file("extdata","GFF_peaks.gff",package="ChIPpeakAnno")
gr2<-toGRanges(gff,format="GFF",header=FALSE,skip=3)#将gff文件转换为GRanges
利用函数进行峰的重合：
其中主要的参数ignore.strand来设置是否忽视来自哪条链；maxgap主要指的在函数findOverlaps中的最大匹配单元。
ol<-findOverlapsOfPeaks(gr1,gr2)#找到重合的peaks
此处也可以绘制相应的韦恩图，当然绘制韦恩图方式很多，我们在此只展现形式，美观嘛，暂时先不考虑，如果感兴趣可以深入研究下：
makeVennDiagram(ol)#绘制veen图
接下来我们导入基因注释文件：
annoData<-toGRanges(EnsDb.Hsapiens.v75,feature="gene")#基因的注释文件
获取数据中重叠的peaks：
overlaps<-ol$peaklist[["gr1///gr2"]]#提取overlap peaks
获取其注释信息：
overlaps.anno<-annotatePeakInBatch(overlaps,
AnnotationData = annoData,
output ="nearestBiDirectionalPromoters",
bindingRegion = c(-2000,500))#注释peak
对注释的信息进行entrez_id注释：
library(org.Hs.eg.db)
overlaps.anno<-addGeneIDs(overlaps.anno,
"org.Hs.eg.db",
IDs2Add="entrez_id")#加入entrez_id，其他也可以
注释信息的可视化：
binOverFeature(overlaps,annotationData =annoData,#计算来自位点的总体的峰值
radius = 5000,nbins=20,FUN=length,errFun=0,#位点数,bin数
ylab="count",
main="distrbution ofaggregated peak numbers around TSS")
pie1(table(overlaps.anno$insideFeature)) #饼图显示peaks的分布
接下来获取peaks在染色体上各区域所占的比例：
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
aCR<-assignChromosomeRegion(overlaps,nucleotideLevel=FALSE, #多个feature总结peaks
precedence=c("Promoters", "immediateDownstream",
"fiveUTRs", "threeUTRs",
"Exons", "Introns"), #设定优先级
TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)#返回的是比例
同样可以利用我们前面的rectome.db进行一个GO/pathway的富集分析，我们在此不再细讲，安装好那个包后大家可以直接利用下面的代码对我们注释好的信息进行富集分析：
path<-getEnrichedPATH(overlaps.anno,"org.Hs.eg.db","reactome.db",maxP=.05)#path富集分析