這里是佳奧，讓我們繼續(xù)作圖的學(xué)習(xí)！

第八步，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

有時(shí)候會(huì)遇見(jiàn)我們實(shí)際使用的參考基因組與作者使用的參考基因組不一樣的情況。

本例子使用的是dm6參考基因組，而作者使用的是dm3，這時(shí)候就需要基因組版本轉(zhuǎn)換。

2 版本轉(zhuǎn)換、重畫(huà)韋恩圖

回到.bam文件：

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

參考人的方法：必應(yīng)搜索 dm6 to dm3 chromosome position

dm3轉(zhuǎn)dm6：需要上傳.bed文件

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

第九步，liftover。

dm3轉(zhuǎn)其他基因組的網(wǎng)站：下載dm3ToDm6.over.chain.gz，在R中會(huì)用到。

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

重新畫(huà)韋恩圖：

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

有了轉(zhuǎn)化dm6以后的.bed文件，我們開(kāi)始畫(huà)最后的韋恩圖。

第十步，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

還是與文章數(shù)目略有差別，比例相似。

##都不在的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也可以做。通過(guò)三個(gè).bed文件做bed overlap。

具體教程：

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

需要文件：具體格式還需要sort，最后得到的就是上面的兩張圖，一張?jiān)谝粡埐辉凇?/p>

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

下一步我們對(duì).bed文件進(jìn)行可視化。

我們下一篇再見(jiàn)！