这里是佳奥，让我们继续作图的学习！

第八步，peaks-overlap。

1 韦恩图（作者处理的.bed，参考dm3）

rm(list=ls())
require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
(bedFiles=list.files(pattern = '*bed'))

library(ChIPpeakAnno)
fs=lapply(bedFiles,function(x){
  peak <- readPeakFile( x )  
  keepChr= !grepl('Het',seqlevels(peak)) 
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  peak
})
library(stringr)
ol <- findOverlapsOfPeaks(fs[[1]],fs[[2]],fs[[3]] )
png('factors_overlapVenn.png')
makeVennDiagram(ol,
                NameOfPeaks=str_split(bedFiles,'_',simplify = T)[,1],
                TxDb=txdb)
dev.off()

3_factors_overlapVenn.png

有时候会遇见我们实际使用的参考基因组与作者使用的参考基因组不一样的情况。

本例子使用的是dm6参考基因组，而作者使用的是dm3，这时候就需要基因组版本转换。

2 版本转换、重画韦恩图

回到.bam文件：

macs2 callpeak - t antiH3K27ry-1.bam -C antiH3K27ry-1-Input.bam -f BAM -g dm -n H3K27 -q 0.05

参考人的方法：必应搜索 dm6 to dm3 chromosome position

dm3转dm6：需要上传.bed文件

https://genome.ucsc.edu/cgi-bin/hgLiftOver

第九步，liftover。

dm3转其他基因组的网站：下载dm3ToDm6.over.chain.gz，在R中会用到。

https://hgdownload-test.gi.ucsc.edu/goldenPath/dm3/liftOver/

# https://genome.ucsc.edu/cgi-bin/hgLiftOver
library(rtracklayer)
path='dm3ToDm6.over.chain'
ch = import.chain(path)
ch
 
require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
library(ChIPpeakAnno)
peak_dm3 <- readPeakFile('Pho_WT.narrowPeak.bed')  

seqlevelsStyle(peak_dm3) = "UCSC"  # necessary
peak_dm6 = liftOver(peak_dm3, ch)
class(peak_dm6)
peak_dm3
peak_dm6

重新画韦恩图：

rm(list=ls())
require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
(bedFiles=list.files(pattern = '*dm6'))

library(ChIPpeakAnno)
fs=lapply(bedFiles,function(x){
  peak <- readPeakFile( x )  
  keepChr= !grepl('Het',seqlevels(peak)) 
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  peak
})
library(stringr)
ol <- findOverlapsOfPeaks(fs[[1]],fs[[2]],fs[[3]] )
png('new3_factors_overlapVenn.png')
makeVennDiagram(ol,
                NameOfPeaks=str_split(bedFiles,'_',simplify = T)[,1],
                TxDb=txdb)
dev.off()

no.2_new3_factors_overlapVenn.png

有了转化dm6以后的.bed文件，我们开始画最后的韦恩图。

第十步，veen-double。

##都在的peaks
rm(list=ls())
require(ChIPseeker) 
require(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
require(clusterProfiler) 
(bedFiles=list.files(pattern = '*dm6'))
bedFiles=bedFiles[-2]

library(ChIPpeakAnno)
fs=lapply(bedFiles,function(x){
  peak <- readPeakFile( x )  
  keepChr= !grepl('Het',seqlevels(peak)) 
  seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
  peak
})
library(stringr)
str_split(bedFiles,'_',simplify = T)[,1]
H3K27 <- readPeakFile('H3K27_peaks.narrowPeak')  
keepChr= seqlevels(H3K27) %in% c('2L','2R','3L','3R','4','X','Y','M')
seqlevels(H3K27, pruning.mode="coarse") <- seqlevels(H3K27)[keepChr]
seqlevels(H3K27, pruning.mode="coarse") <- paste0('chr',seqlevels(H3K27))
seqlevels(H3K27)
tmp1=findOverlapsOfPeaks(fs[[1]], H3K27)
tmp2=findOverlapsOfPeaks(fs[[2]], H3K27)
tmp3=findOverlapsOfPeaks(fs[[3]], H3K27)

ol <- findOverlapsOfPeaks(tmp1$peaklist$`fs..1..///H3K27`,
                          tmp2$peaklist$`fs..2..///H3K27`,
                          tmp3$peaklist$`fs..3..///H3K27`)
png('3_factors_overlapVenn_in_domain.png')
makeVennDiagram(ol,
                NameOfPeaks=str_split(bedFiles,'_',simplify = T)[,1],
                TxDb=txdb)
dev.off()

no.3_3_factors_overlapVenn_in_domain.png

还是与文章数目略有差别，比例相似。

##都不在的peaks
lapply(1:7,function(i){
  gr=ol$peaklist[[i]]
  dat=as.data.frame(gr)[,1:3]
  dat[,1]=gsub('chr','',dat[,1])
  write.table(dat,sep = '\t',
              col.names = F,file = paste0('G',i,'_in.bed')
               ,quote = F,row.names = F)
})

ol <- findOverlapsOfPeaks(tmp1$peaklist$fs..1..,
                          tmp2$peaklist$fs..2..,
                          tmp3$peaklist$fs..3..)
png('3_factors_overlapVenn_not_domain.png')
makeVennDiagram(ol,
                NameOfPeaks=str_split(bedFiles,'_',simplify = T)[,1],
                TxDb=txdb)
dev.off()

no.4_3_factors_overlapVenn_not_domain.png

这样的功能bedtools也可以做。通过三个.bed文件做bed overlap。

具体教程：

https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

需要文件：具体格式还需要sort，最后得到的就是上面的两张图，一张在一张不在。

$ ls
cg.dm6.3field.bed H3K27_peaks.narrowPeak  pho.dm6.3field.bed  spps.dm6.3field.bed

下一步我们对.bed文件进行可视化。

我们下一篇再见！