首页 > 组学教程 > 用ChIPseeker对ChIP-seq数据进行可视化,图表直观颜值高
2023
01-12

用ChIPseeker对ChIP-seq数据进行可视化,图表直观颜值高

ChIP-seq主要用来研究蛋白质和DNA的相互作用,ChIPseeker 可以用来对ChIP-seq数据进行注释与可视化,下面我们就来介绍一下如何用ChIPseeker对chip-seq数据进行可视化操作。

操作步骤

把所有sample_peaks文件放在工作路径下,格式为

#安装程序

#source("http://bioconductor.org/biocLite.R") #biocLite("ChIPseeker") #biocLite("TxDb.Athaliana.BioMart.plantsmart28")

#开始分析

library("ChIPseeker") library(ggplot2) library("TxDb.Athaliana.BioMart.plantsmart28") txdb <- TxDb.Athaliana.BioMart.plantsmart28 get_files_data<-function () { dir <- getwd() files <- list.files(dir) protein <- gsub(pattern = "GSM\\d+_(\\w+_\\w+)_.*", replacement = "\\1", files) protein <- sub("_Chip.+", "", protein) res <- paste(dir, files, sep = "/") res <- as.list(res) names(res) <- protein return(res)} files<-get_files_data()

#读入启动子信息

promoter <- getPromoters(TxDb=txdb,upstream=3000, downstream=3000)

#画覆盖图

covplot(files[[4]])

#画两个样品的覆盖图

peak1=GenomicRanges::GRangesList(sample1=readPeakFile(files[[1]]),sample2=readPeakFile(files[[2]])) covplot(peak1, weightCol="V5") + facet_grid(chr ~ .id)

#放大一个区域

covplot(peak1, weightCol="V5", chrs=c("1", "2"), xlim=c(1, 100000)) + facet_grid(chr ~ .id)

#画其中一个热图

peak <- readPeakFile(files[[4]]) tagMatrix <- getTagMatrix(peak, windows=promoter) tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

#画结合区域热图

peakHeatmap(files,weightCol="V5",TxDb=txdb,upstream=3000,downstream=3000,color=rainbow(length(files)))

如图所示,前三个样品结合位点不在转录起始区域,后两个在。说明前三个样品不是调控转录的。

#结合强度图

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),xlab="Genomic Region (5’->3’)",ylab = "Read Count Frequency")

可以看到转录起始区域结合强度最高

#直接以文件做图,结果与上图一致

plotAvgProf2(files[[4]], TxDb=txdb,upstream=3000, downstream=3000,xlab="Genomic Region (5’->3’)",ylab = "Read Count Frequency")

#加上bootstrap值

#所有数据一起做图

tagMatrixList <- lapply(files,getTagMatrix,windows=promoter) plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

#分面,同时计算置信区间

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000),conf=0.95,resample=500, facet="row")

#画饼状图

peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000),TxDb=txdb, annoDb="org.At.tair.db") plotAnnoPie(peakAnno)

#画柱状图

plotAnnoBar(peakAnno)

#画维恩饼图

vennpie(peakAnno)

#画 upset图

upsetplot(peakAnno)

#结合起来

upsetplot(peakAnno, vennpie=TRUE)

#离最近基因的距离

plotDistToTSS(peakAnno,title="Distribution of transcription factor-binding loci\nrelative to TSS")

#同时画多个柱状图

peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb,tssRegion=c(-3000, 3000), verbose=FALSE) plotAnnoBar(peakAnnoList)

#同时画多个距离图

plotDistToTSS(peakAnnoList)

#画维恩图

genes <- lapply(peakAnnoList, function(i)as.data.frame(i)$geneId) vennplot(genes[2:4], by="Vennerable")

#输出注释信息

peakAnnotable<-as.data.frame(peakAnno) write.table(peakAnnotable,file="peakanno.result",sep="\t",quote=F)



最后编辑:
作者:萌小白
一个热爱网络的青年!

发布评论

表情