chip画图

library(ChIPseeker)
library(ggplot2)
library(dplyr)
library(org.Mm.eg.db)

require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene#TxDb文件

#f <- getSampleFiles()[[4]]#.bed文件,这里用包里的测试数据
DM1 <- "C:\\Users\\admin\\Desktop\\sbed\\sum\\WT_DM1_rmdup_summits.bed"

DM2 <- "C:\\Users\\admin\\Desktop\\sbed\\sum\\WT_DM2_rmdup_summits.bed"

GM1 <- "C:\\Users\\admin\\Desktop\\sbed\\sum\\WT_GM1_rmdup_summits.bed"

GM2 <- "C:\\Users\\admin\\Desktop\\sbed\\sum\\WT_GM2_rmdup_summits.bed"

files <-list(DM1=DM1,DM2=DM2,GM1=GM1,GM2=GM2)

library(devtools)
#install_github("js229/Vennerable",force = TRUE)

#require(Vennerable, character.only = TRUE)
#library(devtools)
#install_github("js229/Vennerable",force = TRUE)

#绘制第1个样本peaks在染色体上的分布
covplot(files[[1]])

#分别绘制不同样本上的peaks分布图

require(GenomicFeatures)

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

peak2<-GenomicRanges::GRangesList(GM1=readPeakFile(files[[3]]),
                                 GM2=readPeakFile(files[[4]]))
covplot(peak2, weightCol="V5") + facet_grid(chr ~ .id)

#关注部分染色体
covplot(peak1, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)

covplot(peak2, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)

## GRangesList is also supported and can be used to compare coverage of multiple bed files.##GRangesList可用来比较多个bed文件。
peakDM=GenomicRanges::GenomicRangesList(DM1=readPeakFile(files[[1]]),DM2=readPeakFile(files[[2]]))
p=covplot(peakDM)
print(p)

#Peak可视化
peak <-readPeakFile(files[[1]])

##getPromoters函数准备启动子窗口
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
##getTagMatrix函数把peaks比对到启动子窗口上
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim =c(-3000, 3000), color="red") #一个图

peakHeatmap(files, TxDb=txdb, upstream=5000, downstream=5000, color=rainbow(length(files)))#正确的 有四个图

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency") #结合强度图

#用peakheatmap函数查看多个样本TSS上下游3kb区域内peak的分布情况：
#tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)

#peakHeatmap(tagMatrixList, TxDb=txdb, upstream=3000, downstream=3000, color=rainbow(length(files)))

#peakHeatmap(files, TxDb=txdb, upstream=3000, downstream=3000, color=rainbow(length(files)))

##lapply函数进行批量处理
tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
##conf定义置信区间，facet决定从上到下还是从左往右
plotAvgProf(tagMatrixList, xlim =c(-3000, 3000), conf=0.95, resample=500, facet="row")

#注释主体

x <- annotatePeak(files[[1]], tssRegion=c(-1000, 1000), TxDb=txdb)#x为输出对象

as.GRanges(x) %>% head(5)
y <- annotatePeak(files[[1]], tssRegion=c(-1000, 1000),addFlankGeneInfo = TRUE,flankDistance = 5000,annoDb='org.Mm.eg.db',TxDb=txdb)#x为输出对象

as.GRanges(y) %>% head(5)
##write.table(y,"C:\\Users\\admin\\Desktop\\DM1.csv",row.names=FALSE,col.names=TRUE,sep=",")

peakAnno <- annotatePeak(files[[1]], tssRegion =c(-3000, 3000),TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)

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

####附加功能
'''
library(clusterProfiler)
##另外，之前的功能注释，只是注释一个蛋白ChIP峰，其实也可以多个ChIP实验进行比较：
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)

names(genes) = sub("_", "\n", names(genes))
compKEGG <- compareCluster(geneCluster   = genes,    
                           data=?
                           fun           = "enrichKEGG",
                           pvalueCutoff  = 0.05,
                           pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")
'''
###

##以上就是ChIPseeker的全部可视化功能了，下面画个不同ChIP实验的峰的重叠图（韦恩图）就结束了。

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)

vennplot(genes,by="gplots")
 vennplot(genes[1:4],by="Vennerable") #注意这里没有"Vennerable"包/函数。我没安装上。所以画不出来，画出来后应该是彩色的。

 #ID转换 获取基因名字
 
 s.SYMBOL<- bitr(z$geneID, fromType="ENTREZID", toType= "SYMBOL",OrgDb="org.Mm.eg.db")
 
 head(s.SYMBOL)
 
 DM1.SYMBOL<- bitr(genes$ DM1, fromType="ENTREZID", toType= "SYMBOL",OrgDb="org.Mm.eg.db")
 head(DM1.SYMBOL)
 
 DM2.SYMBOL<- bitr(genes$ DM2, fromType="ENTREZID", toType= "SYMBOL",OrgDb="org.Mm.eg.db")
 head(DM2.SYMBOL)
 
 GM1.SYMBOL<- bitr(genes$ GM1, fromType="ENTREZID", toType= "SYMBOL",OrgDb="org.Mm.eg.db")
 head(GM1.SYMBOL)
 
 GM2.SYMBOL<- bitr(genes$ GM2, fromType="ENTREZID", toType= "SYMBOL",OrgDb="org.Mm.eg.db")
 head(GM2.SYMBOL)
 

#####go富集合和kegg富集
 library(DOSE)
 #GO富集分析
 
 ego_cc <- enrichGO(gene = DM2.SYMBOL[,1], #使用entrezID作为输入
                    
                    OrgDb=org.Mm.eg.db,
                    
                    ont = "CC",
                    
                    pAdjustMethod = "BH",
                    
                    minGSSize = 1,
                    
                    pvalueCutoff = 0.05,
                    
                    qvalueCutoff = 0.05,
                    
                    readable = TRUE
                    
 )
 
 barplot(ego_cc, showCategory=15, title="GO_Enrichment") #条状图，按p从小到大排的
 
 dotplot(ego_BP,title="EnrichmentGO_CC_dot") #点图，按富集的数从大到小的
 
 
 #write.table(as.data.frame(ego_cc@result),file=".......test_CC.txt",sep="\t")
 
 #KEGG富集分析
 
 kk <- enrichKEGG(gene = geneEntrezID[,2],
                  
                  organism ="mouse",
                  
                  pvalueCutoff = 0.05,
                  
                  qvalueCutoff = 0.01,
                  
                  minGSSize = 1,
                  
                  use_internal_data =FALSE
                  
 )
 
# write.table(as.data.frame(kk@result), file="test_kk.txt",sep="\t")
 dotplot(kk, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")
 
 #下面是另外一个版本的,这个命令是人的数据可以运行
 ##另外，之前的功能注释，只是注释一个蛋白ChIP峰，其实也可以多个ChIP实验进行比较：
 genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
 
 names(genes) = sub("_", "\n", names(genes))
 compKEGG <- compareCluster(geneCluster   = genes,
                            fun           = "enrichKEGG",organism="mmu",
                            pvalueCutoff  = 0.05,
                            pAdjustMethod = "BH")###这个命令是人的数据可以运行##小鼠的要加上mmu才行，默认是人的
 dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")
 
 
 compGO <- compareCluster(geneCluster   = genes,
                            fun           = "enrichGO",OrgDb='org.Mm.eg.db',
                            pvalueCutoff  = 0.05,
                            pAdjustMethod = "BH")###这个命令是人的数据可以运行##小鼠的要加上mmu才行，默认是人的
 dotplot(compGO, showCategory = 15, title = "GO Enrichment Comparison")
 
 
 
 
 
 
 
 
 
 
 
 
 

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

plotDistToTSS(peakAnnoList)

genes <- lapply(peakAnnoList, function(i) 
  as.data.frame(i)$geneId)
vennplot(genes[1:4])

#######part2

##由于某些注释重叠，因此用户可能希望查看具有重叠的完整注释，这可以通过vennpie功能部分解决。
vennpie(peakAnno)
##另外用upsetR可以查看完整注释和重叠
upsetplot(peakAnno)
##两张图可以结合在一起
 upsetplot(peakAnno, vennpie=TRUE) #此图老报错说我Rstudio画板太小。。。。

## TF结合位点相对于TSS的可视化分布
## 前面的annotatePeak已经计算了从峰（结合位点）到最近基因的TSS的距离。接下来可以用plotDistToTSS函数可视化其分布。
#plotDistToTSS(peakAnno,
#              title="Distribution of transcription factor-binding loci\nrelative to TSS")
## 从下图可以看出20%多的peaks（结合位点）离邻近基因的TSS距离小于1kb，而50%的结合位点离TSS的距离大于100kb。

##说实话这几个分布图，我很难理解其意义，如何从图中找出此图的目的，及作为实验研究的切入点，我感到困惑，还需多看文献。

##一旦获得带注释的最近基因，我们就可以通过整合生物本体论提供的生物学知识，进行功能富集分析，以鉴定这些基因中的主要生物学主题。例如例如，GO将基因注释为生物过程、分子功能和细胞成分。KEGG将基因注释为途径，DO将基因与人类疾病联系起来，而Reactome则将基因注释为途径和反应。下面的例子是以Reactome进行演示，其余GO、KEGG、DO亦可。
 
 #BiocManager::install('ReactomePA')
 
 

#library(ReactomePA)####这个是人类的才有用
##注意，enrichPathway是ReactomePA包里的一个函数，因此富集的是Reactome，而非KEGG。
#pathway1 <- enrichPathway(as.data.frame(peakAnno)$geneId)
#head(pathway1, 2)
#dotplot(pathway1)

###用测序片段进行功能富集！！此处的原理类似seq2pathway这个包，详细情况可以去Y叔的公众号看一下
gene <- seq2gene(peak, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=txdb)
pathway2 <- enrichPathway(gene)
head(pathway2, 2)

###前面的“plotAvgProf”和“tagHeatmap”函数可以接受tagMatrix文件并可视化ChIP实验。而“plotAvgProf2”和“peakHeatmap”可以接受bed文件，并作同样的工作。

# promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
# tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
##
## to speed up the compilation of this vigenette, we load a precaculated tagMatrixList

#data("tagMatrixList")
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

#读入具体文件
DM1bed <- readPeakFile(DM1)
DM1V5<-DM1bed@elementMetadata@listData$V5

DM2bed <- readPeakFile(DM2)
DM2V5<-DM2bed@elementMetadata@listData$V5

GM1bed <- readPeakFile(GM1)
GM1V5<-GM1bed@elementMetadata@listData$V5

GM2bed <- readPeakFile(GM2)
GM2V5<-GM2bed@elementMetadata@listData$V5

pcagmdm <- list(DM1V5=DM1V5,DM2V5=DM2V5, GM1V5=GM1V5,GM2V5=GM2V5)

#package ggbiplot
#install.packages("devtools", repos="http://mirror.lzu.edu.cn/CRAN/") #指定镜像
#library(devtools)
#install_github("vqv/ggbiplot",force=T)
library(ggbiplot)

#用prcomp函数进行PCA分析，需要对数据进行scale即均一化。
pcagmdm.pca <- prcomp(pcagmdm, scale. = TRUE)
summary(pacgmdm.pca)