1.Load Packages
library(GenomicRanges)
library(rtracklayer)
library(IRanges)
library(DiffBind)
library(tidyverse)
library(chipseq)
library(Gviz)
library(VennDiagram)
library(BSgenome.Mmusculus.UCSC.mm10)
library(biomaRt)
library(ggplot2)
setwd(“D:/CHIP/promoter”)
##########################################################################
Prepare Mm10 gene Annotation file
mart = useMart(biomart = “ensembl”,
dataset = “mmusculus_gene_ensembl”,
host=”uswest.ensembl.org”)
listAttributes(mart)[1:30,1]
chroms=c(1:19,’X’,’Y’)
Mm10Annot <- getBM(attributes = c(‘ensembl_gene_id’,’external_gene_name’,
‘chromosome_name’,’start_position’,
‘end_position’,’strand’),
filters=’chromosome_name’,
values=chroms,
mart=mart)
save(Mm10Annot,file=”Mm10Annot.RData”)
#####################################################################
Load Mm10Annot
load(“./Mm10Annot.RData”)
head(Mm10Annot)
Mm10Annot$TSS = ifelse(Mm10Annot$strand == “1”,Mm10Annot$start_position,Mm10Annot$end_position )
head(Mm10Annot)
promoter_regions = GRanges(seqnames = Rle( paste0(‘chr’, Mm10Annot$chromosome_name) ),
ranges = IRanges( start = Mm10Annot$TSS -2000,
end = Mm10Annot$TSS + 500),
strand = Rle( rep(“*”, nrow(Mm10Annot)) ),
gene = Mm10Annot$external_gene_name)
promoter_regions
save(Mm10Annot,file=”promoter_regions.RData”)
#####################################################################
Load promoter_regions
load(“./promoter_regions.RData”)
head(promoter_regions)
#####################################################################
Load bed files
flag_wt_rep1 <- import.bed(“./WT_FLAG.core1.bed”)
flag_wt_rep2 <- import.bed(“./WT_FLAG.core2.bed”)
#####################################################################
target <- import.bed(“./core.region.5k.bed”)
target
#####################################################################
load(“./siMm10.RData”)
siMm10
#####################################################################
center <- (start(target)+end(target))/2
name <- rep(paste0(‘S’, 1:length(target)))
chr <- seqnames(target)
value <- rep(1,length(target))
pos.center <- data.frame(name,center,value,chr)
pos.center[1:3,]
tiles = sapply( 1:nrow(pos.center), function(i)
if( pos.center$value[i] == “1” )
pos.center$center[i] + seq( -10000, 9900, length.out=200 )
else
pos.center$center[i] + seq( 9900, -10000, length.out=200 ) )
tiles = GRanges(tilename = paste(rep(pos.center$name, each=200), 1:200, sep=”_” ),
seqnames = Rle(rep(pos.center$chr,each=200)),
ranges = IRanges(start = as.vector(tiles),
width = 100),
strand = Rle(rep(“*”, length(as.vector(tiles)))),
seqinfo = siMm10)
tiles
#####################################################################
flag.WT.p = countOverlaps(tiles,flag_wt_rep1) +
countOverlaps(tiles,flag_wt_rep2)
flag.WT.p.matrix = matrix(flag.WT.p, nrow=nrow(pos.center),
ncol=200, byrow=TRUE )
x <- flag.WT.p.matrix
rowSums(x)
#####################################################################
library(tidyverse)
library(hrbrthemes)
library(viridis)
#####################################################################
name <- c(rep(“A”,4537),rep(“B”,4537))
value <- c(x,x)
data <- data.frame(name,value)
pbox <- ggplot(data,aes(x=name,y=value,fill=name))+
geom_boxplot()
pbox
#####################################################################
