![]([https://files.mdnice.com/user/47194/8d13c60d-023e-44fb-8d96-26da2fcf53c3.gif](https://files.mdnice.com/user/47194/8d13c60d-023e-44fb-8d96-26da2fcf53c3.gif))
背景
展示感兴趣的基因附近的ChIP-seq、DNase/ATAC-seq或RNA-seq的bigwig信号图,不使用IGV或UCSC genome browser那种截图,画出同款矢量图,接着可以用AI等软件编辑。
![]([https://files.mdnice.com/user/47194/2604d117-7783-44a6-be59-a9bd4f8e5bf0.gif](https://files.mdnice.com/user/47194/2604d117-7783-44a6-be59-a9bd4f8e5bf0.gif))
使用场景
- 场景一:用DNase/ATAC-seq或H3K4me3、H3K4me1、H3K27ac等组蛋白修饰的ChIP-seq数据证实基因启动子、增强子的位置。
- 场景二:展示ChIP-seq数据,证实哪些转录因子调控我的基因。
- 场景三:展示RNA-seq和ChIP-seq的信号,证实转录因子结合对基因转录的影响。
输入数据
需要展示的区域loci.bed
bigwig文件TAL1.bw,POLR2A.bw
bigwig文件描述20230831.txt
开始画图
rm(list = ls()) #### 魔幻操作,一键清空~
getwd()
setwd('D:/rtmp/genomeview/')
library(data.table)
library(Gviz)
library(RColorBrewer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
grt<-GeneRegionTrack(TxDb.Hsapiens.UCSC.hg38.knownGene)
grt@dp@pars$transcriptAnnotation<-"symbol"
# 输入数据
bwInfo<-read.table("20230831.txt",header=F,row.names=1,as.is=T)
head(bwInfo)
V2
POLR2A POLR2A_FetalLiver_GSM970259
TAL1 TAL1_FetalLiver_GSM1427077
gloci<-read.table("loci.bed",header=F,as.is=T)
head(gloci)
V1 V2 V3 V4
1 chr1 37996770 38005505 -
#可调整的参数
genefold<-as.numeric("1.5")#放大、缩小展示的范围
# 展示的基因组范围
colnames(gloci)<-c("chr","start","end","strand")
chr<-gloci[rownames(gloci),]$chr
gloci$width<-with(gloci,end-start)
startpoint<-gloci[rownames(gloci),]$start-genefold*gloci[rownames(gloci),]$width
endpoint<-gloci[rownames(gloci),]$end+genefold*gloci[rownames(gloci),]$width
#下面将scale等track写入tracklist
tracklist<-list()
#写入chromosome
itrack <- IdeogramTrack(genome = "hg38", chromosome = chr,outline=T)
tracklist[["itrack"]]<-itrack
#写入比例尺
scalebar <- GenomeAxisTrack(scale=0.25,col="black",fontcolor="black",name="Scale",labelPos="above",showTitle=TRUE)
tracklist[["scalebar"]]<-scalebar
#写入基因组位置
axisTrack <- GenomeAxisTrack(labelPos="above",col="black",fontcolor="black",name=paste(chr,":",sep=""),exponent=0,showTitle=TRUE)
tracklist[["axisTrack"]]<-axisTrack
#写入bigwig
#配色
colpal<-rep(brewer.pal(12,"Paired"),20)
coldf<-data.frame(col=colpal[1:nrow(bwInfo)],row.names = rownames(bwInfo),stringsAsFactors = F)
for(index in rownames(bwInfo)){
bgFile<-file.path("D:/rtmp/genomeview/bw/",paste(index,".bw",sep=""))
tracklist[[index]]<-DataTrack(range = bgFile,genome="hg38",type="histogram",
name=chartr("_","\n",bwInfo[index,]),
col.histogram=coldf[index,])#每个track颜色不同才好看
}
#写入基因结构
tracklist[["grt"]]<-grt
#画图
plotTracks(tracklist, from = startpoint, to = endpoint,
chromosome=chr,background.panel = "white", background.title = "white",
col.title="black",col.axis="black",
rot.title=0,cex.title=0.9,margin=38,title.width=1.5,
collapseTranscripts = "longest")#同一个基因的多个transcript压缩成最长的一个
#输出pdf文件
pdf("20230831.pdf",height=length(tracklist)+1,width=12)
plotTracks(tracklist, from = startpoint, to = endpoint,
chromosome=chr,background.panel = "white", background.title = "white",
col.title="black",col.axis="black",
rot.title=0,cex.title=0.9,margin=38,title.width=1.5,
collapseTranscripts = "longest")#同一个基因的多个transcript压缩成最长的一个
dev.off()