Commit fee2ae63 authored by San Leon Granado's avatar San Leon Granado
Browse files

Upload New File

parent dc05039a
---
title: "R Notebook"
output: html_notebook
---
# Create plots for every peak in Chipseek
```{r echo=FALSE}
library(trackViewer)
library(ggplot2)
library(rtracklayer)
library(Gviz)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
txdb=TxDb.Mmusculus.UCSC.mm10.knownGene
```
# Load peaks with 1Kb to peak sides
```{r}
chipPeaks=import("Zbtb24mESCtrack2KbAroundPeak.bed")
```
```{r loadTracks}
Zbtb24mESCtrack=import("Zbtb24mESCtrack.bed")
Zbtb24MotifPositions=import("Zbtb24MotifPositions.bed")
CpG=import("mm10.CpG.onlyExtendedRegions.bed")
GpGislands=import("mm10.CpGisland.extended.bed")
CpGi=AnnotationTrack(GpGislands,genome="mm10",
name="CpG Islands",
fill="black",
background.title="transparent",
col.frame="transparent",
col.axis="black",
col="black",
col.title="black")
options(ucscChromosomeNames=FALSE)
methAverage=read.table("methAverage.csv",sep="\t",header=FALSE)
Gene=GeneRegionTrack(txdb,genome="mm10",
name="Gene",
background.title="transparent",
fill="black",collapseTranscripts = "longest",
col.frame="transparent", col.axis="black",transcriptAnnotation = "symbol",
col="black", col.title="black",rot.title=0)
symbols <- unlist(mapIds(org.Mm.eg.db, gene(Gene), "SYMBOL", "ENTREZID", multiVals = "first"))
symbol(Gene) <- symbols[gene(Gene)]
# methylAverage=AnnotationTrack(start=to-200,width=200,chromosome=mychr,
# id=as.character(round(methAverage$V4[index],digits=2)),
# featureAnnotation="id",
# genome="mm10",fontcolor.feature="black",fill="white",
# background.title="transparent",
# col.frame="transparent")
```
```{r}
sizeAxis=0.8
axis=TRUE
#par(mfrow=c(4,4),mar=c(3,1,1,1))
#options(Gviz.scheme = "default")
#pdf("tmp/A4.allpeaks.pdf",paper="a4r",family="sans")
for (index in seq(1,40)){
print(index)
mychr= as.character(seqnames(chipPeaks[index])@values)
myRegion=ranges(chipPeaks[index])
gr=chipPeaks[index]
from=myRegion@start
to=myRegion@start+myRegion@width
axisTrack <- GenomeAxisTrack()
myMethylationPositive=import("WTvsZbtb24.UCSC.onlyuExtendedPeaks.positives.bedGraph",which=gr ,format="bedGraph")
myMethylationNegative=import("WTvsZbtb24.UCSC.onlyuExtendedPeaks.negatives.bedGraph",which=gr ,format="bedGraph")
myChipSeq=import("Zb24_mES_bwa_mm10q30_treat_pileup.bdg.bw",which=gr)
myH3K4me3=import("H3K4me3.serum.bw",which=gr)
myATAC=import("WT.ATACseq.mESC.serum.bw",which=gr)
H3K4me3TrackHist <- DataTrack(myH3K4me3, ylim=c(0,150),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="hist",
fill.mountain=c("black","black"),
col.mountain="black",
fill.histogram="black", col.histogram="black",
chromosome=mychr,
name="H3K4me3",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
H3K4me3TrackMountain <- DataTrack(myH3K4me3, ylim=c(0,150), cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="mountain",degree=3,
fill.mountain=c("black","black"),
col.mountain="black",
chromosome=mychr,
name="H3K4me3",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
myH3K27me3=import("H3K27me3.serum.bw",which=gr)
H3K27me3TrackHist <- DataTrack(myH3K27me3, ylim=c(0,10),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="hist",
fill.mountain=c("black","black"),
fill.histogram="black", col.histogram="black",
col.mountain="black",
chromosome=mychr,
name="H3K27me3",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
H3K27me3TrackMountain <- DataTrack(myH3K27me3, ylim=c(0,10),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="mountain",degree=3,
fill.mountain=c("black","black"),degree=3,
col.mountain="black",
chromosome=mychr,
name="H3K27me3",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
myPOL2=import("POL2.BRUCE.ENCODE.bigwig",which=gr)
POL2TrackHist <- DataTrack(myPOL2, ylim=c(0,10),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="hist",
fill.mountain=c("black","black"),
col.mountain="black",
fill.histogram="black", col.histogram="black",
chromosome=mychr,
name="POL2",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
POL2TrackMountain <- DataTrack(myPOL2, ylim=c(0,10),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="mountain",degree=3,
fill.mountain=c("black","black"),degree=3,
col.mountain="black",
chromosome=mychr,
name="POL2",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
ATACTrackHist <- DataTrack(myATAC, ylim=c(0,80),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",degree=3,
type="mountain",
fill.mountain=c("black","black"), col.mountain="black",
chromosome=mychr,
name="ATAC-seq",
background.title="transparent", col.axis="black",rot.title=0,window=-1,
col="black", col.title="black")
chipPeaksTrack=AnnotationTrack(Zbtb24mESCtrack,showAxis=axis,
genome="mm10",
chromosome=mychr,
name="Zbtb24 peak",
fill="black",
background.title="transparent",
col.frame="transparent",
col.axis="black",
col="black", rot.title=0,
col.title="black")
CpGTrack=AnnotationTrack(CpG,
genome="mm10",
chromosome=mychr,
name="CpG",
fill="black",
background.title="transparent",
col.frame="transparent",
col.axis="black",
col="black", rot.title=0,
col.title="black",stacking="dense")
chipPeaksMotif=AnnotationTrack(Zbtb24MotifPositions,
genome="mm10",
chromosome=mychr,
name="Zbtb24 motif",
fill="black",
background.title="transparent",
col.frame="transparent",
col.axis="black",
col="black",
col.title="black",rot.title=0
)
chipseq <- DataTrack(myChipSeq,ylim=c(0,80),cex.axis=sizeAxis,showAxis=axis,
genome="mm10",
type="mountain",
fill.mountain=c("black","black"),
col.mountain="black",
chromosome=mychr,
name="ChIPseq",rot.title=0,
background.title="transparent", col.axis="black",
col="black", col.title="black")
dt <- DataTrack(myMethylationPositive,ylim=c(-100,100),cex.axis=sizeAxis,showAxis=axis,
chromosome=mychr,
genome="mm10",
type=c("hist"),
name="Zbtb24",
fill.histogram="red", col.histogram="red",
background.title="transparent",
col.frame="transparent", col.axis="black",
col="black", col.title="black")
displayPars(dt)=list(groups=factor("Hypermethylated",levels=c("Hypermethylated","Hypomethylated")),legend=FALSE,fill="transparent",col="red")
dt2 <- DataTrack(myMethylationNegative,ylim=c(-100,100),cex.axis=sizeAxis,showAxis=axis,
chromosome=mychr,
genome="mm10",
type=c("hist"),
name="Zbtb24",
fill.histogram="blue", col.histogram="blue",
background.title="transparent",
col.frame="transparent", col.axis="black",
col="black", col.title="black")
ot=OverlayTrack(trackList = list(dt,dt2),
background.title="transparent",
col.frame="transparent", col.axis="black",
col.title="black")
cdca7positives=DataTrack("WTvsCdca7.UCSC.onlyuExtendedPeaks.positives.bedGraph",ylim=c(-100,100),cex.axis=sizeAxis,showAxis=axis,
chromosome=mychr,
genome="mm10",
type=c("hist"),
name="Cdca7",
fill.histogram="orange", col.histogram="orange",
background.title="transparent",
col.frame="transparent", col.axis="black",
col="orange", col.title="black")
displayPars(cdca7positives)=list(groups=factor("Hypermethylated",levels=c("Hypermethylated","Hypomethylated")),legend=FALSE,fill="transparent",col="orange")
cdca7negatives=DataTrack("WTvsCdca7.UCSC.onlyuExtendedPeaks.negatives.bedGraph",ylim=c(-100,100),showAxis=axis,
chromosome=mychr,
genome="mm10",
type=c("hist"),
name="Cdca7",
fill.histogram="darkgreen", col.histogram="darkgreen",
background.title="transparent",
col.frame="transparent", col.axis="black",
col="black", col.title="black")
cdca7=OverlayTrack(trackList = list(cdca7positives,cdca7negatives),
background.title="transparent",
col.frame="transparent", col.axis="black",
col.title="black")
#pdf(paste("peak",index,"pdf",sep="."))
errorH3K4=tryCatch({
pdf("tmp.pdf")
a=plotTracks(H3K4me3TrackMountain, chromosome=mychr,from=from, to=to,fontcolor="black")
dev.off()
},error=function(e){e})
errorH3K27=tryCatch({
pdf("tmp.pdf")
b=plotTracks(H3K27me3TrackMountain, chromosome=mychr,from=from, to=to,fontcolor="black")
dev.off()},error=function(e){e})
errorPOL2=tryCatch({
pdf("tmp.pdf")
c=plotTracks(POL2TrackMountain, chromosome=mychr,from=from, to=to,fontcolor="black")
dev.off()},error=function(e){e})
if(any(class(errorH3K4)=="error")){
H3K4me3=H3K4me3TrackHist
}else{
H3K4me3=H3K4me3TrackMountain
}
if(any(class(errorH3K27)=="error")){
H3K27me3=H3K27me3TrackHist
} else{
H3K27me3=H3K27me3TrackMountain
}
if(any(class(errorPOL2)=="error")){
POL2=POL2TrackHist
} else{
POL2=POL2TrackMountain
}
pdf(paste0("tmpFinal/peak",index,".pdf"),width=8,height=8,family="sans")
if(index!=34){
plotTracks(list(axisTrack,CpGTrack,ot,cdca7,chipseq,chipPeaksTrack,chipPeaksMotif,CpGi,Gene,H3K4me3,H3K27me3), chromosome=mychr,from=from, to=to,fontcolor="black",sizes=c(0.5,0.5,2,2,1,0.5,0.5,0.5,0.5,1,1),main=paste0(mychr,":",from,"-",to),scale=1000)
}else{
plotTracks(list(axisTrack,CpGTrack,ot,cdca7,chipseq,chipPeaksTrack,chipPeaksMotif,CpGi,CpGi,H3K4me3,H3K27me3), chromosome=mychr,from=from, to=to,fontcolor="black",sizes=c(0.5,0.5,2,2,1,0.5,0.5,0.5,0.5,1,1),main=paste0(mychr,":",from,"-",to),scale=1000)
}
dev.off()
}
#dev.off()
```
```{r trackViewer}
library(trackViewer)
H3K4me3 <- importScore("H3K4me3.serum.bw",ranges=gr,format="BigWig")
H3K27me3 <- importScore("H3K27me3.serum.bw",ranges=gr,format="BigWig")
POL2 <-importScore("POL2.BRUCE.ENCODE.bigwig",ranges=gr,format="BigWig")
trs <- geneModelFromTxdb(TxDb.Mmusculus.UCSC.mm10.knownGene,
org.Mm.eg.db,
gr=gr)
H3K4me3$dat <- coverageGR(H3K4me3$dat)
H3K27me3$dat <- coverageGR(H3K27me3$dat)
POL2$dat <- coverageGR(POL2$dat)
viewerStyle <- trackViewerStyle()
setTrackViewerStyleParam(viewerStyle, "margin", c(.1, .05, .02, .02))
vp <- viewTracks(trackList(trs,H3K4me3, H3K27me3, POL2),
gr=gr, viewerStyle=viewerStyle,
autoOptimizeStyle=TRUE,smooth=TRUE)
```
```{r correlation}
library(ggpubr)
library(pheatmap)
negatives=read.table("correlationMethylation/CorrelationCdca7_Zbtb24.negatives.bed",sep="\t")
negatives=negatives[negatives$V4< -5 | negatives$V8< -5, ]
positives=read.table("correlationMethylation/CorrelationCdca7_Zbtb24.positives.bed",sep="\t")
positives=positives[positives$V4>5 | positives$V8>5, ]
pheatmap(log2(abs(negatives[,c(4,8)])),cluster_rows = FALSE,cluster_cols = FALSE)
negatives$Zbtb24=log2(abs(negatives$V4))
negatives$Cdca7=log2(abs(negatives$V7))
ggscatter(negatives,x="Zbtb24",y="Cdca7",add="reg.line",add.params = list(color="blue",fill="lightgray",conf.int=TRUE))+
stat_cor(method="pearson",label.x=3,label.y=40)
positives$Zbtb24=log2(abs(positives$V4))
positives$Cdca7=log2(abs(positives$V7))
ggscatter(positives,x="Zbtb24",y="Cdca7",add="reg.line",add.params = list(color="blue",fill="lightgray",conf.int=TRUE))+
stat_cor(method="pearson",label.x=3,label.y=30)
```
# Figure with some examples
```{r}
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment