library("karyoploteR")
library("IRanges")
library("GenomicFeatures")
library(Biobase)
library(dbplyr)
library(AnnotationDbi)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(BiocFileCache)
library(org.Hs.eg.db)
#------------------------------
pdf('1.pdf',width = 11,height = 8.5)
pic<-function(regi){
#文件有两列，第一列是enhancer位置，第二列是enhancer扩展后的要展示的染色体区域
r<-regi[1]#enhancer就在添加region的时候用到了
reg<-regi[2]
enh.region <- toGRanges(reg)
#----------------------------------------
histone.marks <- c(H3K4me3="H3K4me3.bigWig",
H3K4me1="H3K4me1.bigWig",
H3K27ac="H3K27ac.bigWig"
)
DNA.binding <- c(CTCF="CTCF.bigWig",
DNase="DNase.bigWig",
POL2="POLR2A.bigWig",
YY1="YY1.bigWig")
pp <- getDefaultPlotParams(plot.type=1)
pp$leftmargin <- 0.15
pp$topmargin <- 15
pp$bottommargin <- 15
pp$ideogramheight <- 5
pp$data1inmargin <- 10
pp$data1outmargin <- 0
kp <- plotKaryotype(zoom = enh.region, cex=0.7, plot.params = pp)
genes.data <- makeGenesDataFromTxDb(TxDb.Hsapiens.UCSC.hg38.knownGene,
karyoplot=kp,
plot.transcripts = TRUE,
plot.transcripts.structure = TRUE)
genes.data <- addGeneNames(genes.data)
genes.data <- mergeTranscripts(genes.data)
kpAddBaseNumbers(kp, tick.dist = 10000, minor.tick.dist = 2000,
add.units = TRUE, cex=0.7, tick.len = 3)
kpAddMainTitle(kp, "Epigenetic Regulation in HepG2", cex=1)
kpPlotGenes(kp, data=genes.data, r0=0, r1=0.14, gene.name.cex = 0.5)
#添加区域
kpPlotRegions(kp, data=c(r), col="#C1FFC1", r0=0.24, r1=1)
kpPlotRegions(kp, data=c(r), col="#FFC125", r0=0.22, r1=0.24)
#添加文本
kpText(kp,labels = "Enhancer chr22:30212552-30213117", chr="chr22", x=30212552,y=0.2,r0=0.20, r1=0.22,cex=0.6)
#Histone marks
total.tracks <- length(histone.marks)+length(DNA.binding)
out.at <- autotrack(1:length(histone.marks), total.tracks, margin = 0.3, r0=0.25)
kpAddLabels(kp, labels = "Histone marks", r0 = out.at$r0, r1=out.at$r1, cex=1,
srt=90, pos=1, label.margin = 0.14)
for(i in seq_len(length(histone.marks))) {
bigwig.file <- paste0(histone.marks[i])
at <- autotrack(i, length(histone.marks), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
r0=at$r0, r1=at$r1, col = "#3A5FCD")
computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax,
r0=at$r0, r1=at$r1, cex=0.7)
kpAddLabels(kp, labels = names(histone.marks)[i], r0=at$r0, r1=at$r1,
cex=0.7, label.margin = 0.035)
}
#DNA binding proteins
out.at <- autotrack((length(histone.marks)+1):total.tracks, total.tracks, margin = 0.3, r0=0.2)
kpAddLabels(kp, labels = "DNA-binding proteins", r0 = out.at$r0, r1=out.at$r1,
cex=1, srt=90, pos=1, label.margin = 0.14)
for(i in seq_len(length(DNA.binding))) {
bigwig.file <- paste0(DNA.binding[i])
at <- autotrack(i, length(DNA.binding), r0=out.at$r0, r1=out.at$r1, margin = 0.1)
kp <- kpPlotBigWig(kp, data=bigwig.file, ymax="visible.region",
r0=at$r0, r1=at$r1, col = "#FF4040")
computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
kpAxis(kp, ymin=0, ymax=computed.ymax, tick.pos = computed.ymax,
r0=at$r0, r1=at$r1, cex=0.7)
kpAddLabels(kp, labels = names(DNA.binding)[i], r0=at$r0, r1=at$r1,
cex=0.7, label.margin = 0.035)
}
}
#enhancer file
a<-read.csv('0.txt',header=FALSE,sep = '\t')
b<-apply(a,1,pic)
dev.off()
效果图：
0001.jpg
感觉就一个不好理解的参数：r0和r1
r0：是起始的高度
r1：是结束的高度
r1-r0才是我们想要展示内容的高度（想象这个画布本身就有刻度线）
传入的增强子文件格式：前面是enhancer，后面是enhancer上下5kb
enhancer.png
