搜索
查看: 4581|回复: 4

[software] 使用ChippeakAnno的通用下游分析流程

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-5-26 10:32:34 | 显示全部楼层 |阅读模式
[Python] 纯文本查看 复制代码
#1.	单个重组transcription factor的实验流程(整个流程中的数据为GRanges)
#提取重复的peaks
#准备注释信息
#查看绑定位点分布
#注释peaks
#GO,KEGG富集分析
#提取peaks附近的序列
#得到一个peaks的summary
#用双向启动子寻找peaks
#通过DNA interaction data 得到可能的增强子


source("https://bioconductor.org/biocLite.R")
biocLite("ChIPpeakAnno")
library(ChIPpeakAnno)
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")#提取gff文件的全目录
gr2<-toGRanges(gff,format="GFF",header=FALSE,skip=3)#将gff文件转换为GRanges
gr2$score<-as.numeric(gr2$score)#将score变为数值型
ol<-findOverlapsOfPeaks(gr1,gr2)#找到重合的peak,显示的项目有:
#$venn_cnt
#$peaklist
#$peaklist$gr2
#$peaklist$gr1
#$peaklist$`gr1///gr2`
#$uniquePeaks
#$mergedPeaks
#$peaksInMergedPeaks
#$overlappingPeaks
#$overlappingPeaks$`gr1///gr2`
#$all.peaks
#$all.peaks$gr1
#$all.peaks$gr2
ol<-addMetadata(ol,colNames = "score",FUN=mean)#将score加进去
ol$peaklist[["gr1///gr2"]][1:2]#选前两行的peak
makeVennDiagram(ol)#绘制veen图

#1.2 准备注释文件
biocLite("EnsDb.Hsapiens.v75")#注释文件h19
library(EnsDb.Hsapiens.v75)#更新了所有的包才可以
#建立annotation文件
annoData<-toGRanges(EnsDb.Hsapiens.v75,feature="gene")#基因的注释文件,将注释的形式设置为基因
annoData[1:2]

#1.3可视化绑定站点
overlaps<-ol$peaklist[["gr1///gr2"]]#提取overlap peak
binOverFeature(overlaps,annotationData = annoData,#计算来自位点的总体的峰值
               radius = 5000,nbins =20,FUN=length,errFun=0,#位点数,bin数
               ylab="count",
               main="distrbution of aggregated peak numbers around TSS")
source("https://bioconductor.org/biocLite.R")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
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)#返回的是比例
barplot(aCR$percentage,las=3)#绘制比例图

#1.4 注释peak
overlaps.anno<-annotatePeakInBatch(overlaps,
                                   AnnotationData = annoData,
                                   output = "nearestBiDirectionalPromoters",
                                   bindingRegion = c(-2000,500))#注释peak
library(org.Hs.eg.db)
overlaps.anno<-addGeneIDs(overlaps.anno,
                          "org.Hs.eg.db",
                          IDs2Add="entrez_id")#加入entrez_id,别的也可以,比如symbol
head(overlaps.anno)
#饼图显示features中peaks的分布
pie1(table(overlaps.anno$insideFeature))

#1.5GO和pathway富集
#这里的ensembl_gene_id默认的
over<-getEnrichedGO(overlaps.anno,orgAnn = "org.Hs.eg.db",
                    maxP=.05,minGOterm = 10,
                    multiAdjMethod = "BH",condense = TRUE)
head(over[["bp"]][,-c(3,10)])
source("https://bioconductor.org/biocLite.R")
biocLite("reactome.db")
library(reactome.db)
path<-getEnrichedPATH(overlaps.anno,"org.Hs.eg.db","reactome.db",maxP=.05)#path富集分析
head(path)

#1.6 得到peak周围的序列
library(BSgenome.Hsapiens.UCSC.hg19)
seq <- getAllPeakSequence(overlaps, upstream=20, downstream=20, genome=Hsapiens)
write2FASTA(seq, "test.fa")#写入fasta格式文件

#1.7输出峰值的summary
freqs<-oligoFrequency(Hsapiens$chr1,MarkovOrder=3)#得到寡核苷酸的频率
os<-oligoSummary(seq,oligoLength = 6,MarkovOrder = 3,
                 quickMotif = FALSE,freqs=freqs)#统计所有的长度为6的寡核苷酸,Markov chain顺序为3
#绘制结果
zscore<-sort(os$zscore)#根据zscore进行排序
h<-hist(zscore,breaks = 100,xlim=c(-50,50),main="histogram of zscore")#画图
text(zscore[length(zscore)],max(h$counts)/10,labels=names(zscore[length(zscore)]),adj=1)#图中显示的是最后一个寡核苷酸的名字
###尝试模拟数据
seq.sim.motif<-list(c("t","g","c","a","t","g"),
                    c("g","c","a","t","g","c"))#仿真数据
set.seed(1)#只是作为一个标记
seq.sim<-sapply(sample(c(2,1,0),1000,replace=TRUE,prob = c(0.07,0.1,0.83)),#生成1000条序列
                function(x){
                  s<-sample(c("a","c","g","t"),
                            sample(100:1000,1),replace = TRUE)#设置单个序列长度
          if(x>0){
            si<-sample.int(length(s),1)
            if(si>length(s)-6) si<-length(s)-6
              s[si:(si+5)]<-seq.sim.motif[[x]]
          }
                  paste(s,collapse="")
                })#模拟数据1000条序列
os<-oligoSummary(seq.sim,oligoLength = 6,MarkovOrder = 3,
                 quickMotif = TRUE)#其他同上
zscore<-sort(os$zscore,decreasing = TRUE)
h<-hist(zscore,breaks=100,main="Histogram of zscore")
text(zscore[1:2],rep(5,2),
     labels=names(zscore[1:2]),adj=0,srt=90)
##生成motifs
source("https://bioconductor.org/biocLite.R")
biocLite("motifStack")
biocLite("seqLogo")
## generate the motifs
library(motifStack)
pfms <- mapply(function(.ele, id)
  new("pfm", mat=.ele, name=paste("SAMPLE motif", id)), 
  os$motifs, 1:length(os$motifs))
motifStack(pfms[[1]])
plot(pfms[[1]])######画一个DNA序列的logo

#1.8用双向启动子寻找peak
#双向启动子是两个相邻基因的TSS之间的DNA区域,他们在相反的方向被转录并且通常共享的启动子区域共同调节。
bdp<-peaksNearBDP(overlaps,annoData,maxgap=5000)#找到双向启动i子,并且给出接近双向启动子的peaks的比例
c(bdp$percentPeaksWithBDP,
  bdp$n.peaks,
  bdp$n.peaksWithBDP)#比例,peaks总数,临近的peaks


#1.9根据DNA的相互作用数据找到可能的增强子
#找到潜在的与增强子区域结合的peaks
DNA5C <- system.file("extdata", 
                     "wgEncodeUmassDekker5CGm12878PkV2.bed.gz",
                     package="ChIPpeakAnno")#得到实例bed文件
DNAinteractiveData<-toGRanges(gzfile(DNA5C))#格式转换
findEnhancers(overlaps,annoData,DNAinteractiveData )#找到增强子(通过染色体构象捕获技术)




[AppleScript] 纯文本查看 复制代码
#2.比较多个转录因子(TFS)的绑定谱的流程
#给定来自来自不同的TFS的两个或多个peaks列表,我们可能感兴趣在发现DNA绑定表达谱是否与TFs相关
#如果相关,common binding模式。下面的例子测试三个TFS的绑定表达谱的相关性,并且找到common binding模式
#2.1数据输入
library(ChIPpeakAnno)
path<-system.file("extdata",package="ChIPpeakAnno")#找到路径
files<-dir(path,"broadPeak")
data<-sapply(file.path(path,files),toGRanges,format="broadPeak")#named list
names(data)<-gsub(".broadPeak","",files)#gsub是替换所有符合条件的,替换为空格

#2.2查看peaks集有无显著重复
#2.2.1超几何测试
#当我们基于几何分布测试两组数据之间的关联时,需要所有的可能结合位点的数量。
#函数makeVennDiagram中的参数totalTest表示总共有多少潜在的peaks用于几何测试。该数值应该大于peaks中的最大数值
#设置越小,测试越严格。计算p-value取决于totalTest的数值。下面的例子假设存在3%的编码区加启动子区。因为样本只是染色体2的
#一部分,我们估计所有的结合位点时基因组中可能结合的区域的1/24
ol<-findOverlapsOfPeaks(data,connectedPeaks = "keepAll")#找到重复的peaks,结果是列表
averagePeakWidth<-mean(width(unlist(GRangesList(ol$peaklist))))#求peaks的平均宽度
tot<-ceiling(3.3e+9*.03/averagePeakWidth/24)#大于ceiling后面的那个数的整数/平均宽度/24
makeVennDiagram(ol,totalTest = tot,connectedPeaks = "keepAll")#画出两个或多个peak set的venn图,并计算p-value查看是否显著

#2.2.2置换测试
#上面的几何测试需要输入一个对于给定的TFs的估计的所有的潜在的绑定位点数。为了规避这个需要,我们提供一个置换测试peakPermTest.
#在测试之前,用户需要生成随机的peaks基于给定的特征类型的peaks分布。保证绑定位置与特征相关,例如TSS和GENEend,随机peaks的宽度follow
#输入peaks的分布
#利用preparePpool可以建立一个peak pool代表所有的潜在的绑定位点(与绑定相关概率对于随机peak样本)
#下面的例子是建立一个人类基因组的peak pool 利用转录因子绑定位点簇(V3),hot spot是一个很容易被TFs限制的基因组区域,建议去除
data("HOT.spots")
data("wgEncodeTfbsV3")
hotGR<-reduce(unlist(HOT.spots))#
removeol<-function(.ele){
  ol<-findOverlaps(.ele,hotGR)#找到重复,是个列表
  if(length(ol)>0).ele<-.ele[-unique(queryHits(ol))]
  .ele
}#移除hotspot防止过估计两个输入Peak list
TAF<-removeol(data[["TAF"]])
TEAD4<-removeol(data[["Tead4"]])
YY1<-removeol(data[["YY1"]])
#SUBSET the pool为了节省时间
set.seed(1)
wgEncodeTfbsV3.subset<-wgEncodeTfbsV3[sample.int(length(wgEncodeTfbsV3),2000)]
pool<-new("permPool",grs=GRangesList(wgEncodeTfbsV3.subset),N=length(YY1))
pt1<-peakPermTest(YY1,TEAD4,pool=pool,seed=1,force.parallel=FALSE)#求出这两个peaks list置换检验的结果
plot(pt1)

#2.3 绑定模式的对比和可视化
features<-ol$peaklist[[length(ol$peaklist)]]#$`TAF///Tead4///YY1`
feature.recentered<-reCenterPeaks(features,width=4000)#基于peak中心重新建立一个peak list
library(rtracklayer)
files<-dir(path,"bigWig")#给文件的后缀都变为.bigWig
#rtracklayer包不接受bigWig文件在windows,
if(.Platform$OS.type!="windows"){
  cvglists<-sapply(file.path(path,file),import,
                   format="BigWig",
                   which=feature.recentered,
                   as="RleList")
}else {
  load(file.path(path,"cvglist.rds"))
}
names(cvglists)<-gsub(".bigWig","",files)#去除后缀
feature.center<-reCenterPeaks(features,width=1)
sig<-featureAlignedSignal(cvglists,feature.center,
                          upstream = 2000,downstream=2000)#在区间内的信号提取出来
#由于bw文件只是原始数据的子集,所以不是每个Peak都有信号
keep<-rowSums(sig[[2]])>2#Tead4  各行相加,大于2
sig<-sapply(sig,function(.ele).ele[keep,],simplify = FALSE)
feature.center<-feature.center[keep]
heatmap<-featureAlignedHeatmap(sig,feature.center,upstream = 2000,downstream=2000,upper.extreme=c(3,.5,4))#绘制区间内的热图

sig.rowsums<-sapply(sig,rowSums,na.rm=TRUE)#对上面的数据进行聚类
d<-dist(sig.rowsums)
hc<-hclust(d)
feature.center$order<-hc$order#加一列order

heatmap <- featureAlignedHeatmap(sig, feature.center, 
                                 upstream=2000, downstream=2000,
                                 upper.extreme=c(3,.5,4),
                                 sortBy="order")#聚类绘图

featureAlignedDistribution(sig, feature.center, 
                           upstream=2000, downstream=2000,
                           type="l")#绘制分布图


回复

使用道具 举报

365

主题

512

帖子

1713

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1713
发表于 2017-5-26 12:01:39 | 显示全部楼层
可以试试看写成Rmarkdown格式,我看群主都是这样写的,条理非常清晰适宜阅读
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
 楼主| 发表于 2017-5-26 23:06:15 | 显示全部楼层
ydchen 发表于 2017-5-26 12:01
可以试试看写成Rmarkdown格式,我看群主都是这样写的,条理非常清晰适宜阅读 ...

好的好的
回复 支持 反对

使用道具 举报

0

主题

6

帖子

273

积分

中级会员

Rank: 3Rank: 3

积分
273
发表于 2017-7-26 22:03:29 | 显示全部楼层
我跑到这儿总是出现错误,楼主帮我!

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

2

帖子

96

积分

注册会员

Rank: 2

积分
96
发表于 2017-8-17 19:18:47 | 显示全部楼层
你好,ChippeakAnno怎样才能得到每个peak对应的基因组元件信息呢?这一步只能显示比例
aCR<-assignChromosomeRegion(overlaps, nucleotideLevel=FALSE, #多个feature总结peaks
                            precedence=c("Promoters", "immediateDownstream",
                                         "fiveUTRs", "threeUTRs",
                                         "Exons", "Introns"), #设定优先级
                            TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)#返回的是比例
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-19 08:50 , Processed in 0.042233 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.