画出igv款式的矢量图
# 背景
展示感兴趣的基因附近的ChIP-seq、DNase/ATAC-seq或RNA-seq的bigwig信号图，不使用IGV或UCSC genome browser那种截图，画出同款矢量图，接着可以用AI等软件编辑。
# 使用场景
- 场景一：用DNase/ATAC-seq或H3K4me3、H3K4me1、H3K27ac等组蛋白修饰的ChIP-seq数据证实基因启动子、增强子的位置。
- 场景二：展示ChIP-seq数据，证实哪些转录因子调控我的基因。
- 场景三：展示RNA-seq和ChIP-seq的信号，证实转录因子结合对基因转录的影响。
# 输入数据
需要展示的区域loci.bed
bigwig文件TAL1.bw，POLR2A.bw
bigwig文件描述20230831.txt
# 开始画图
```R
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)
```
```R
V2
POLR2A POLR2A_FetalLiver_GSM970259
TAL1 TAL1_FetalLiver_GSM1427077
```
```R
gloci<-read.table("loci.bed",header=F,as.is=T)
head(gloci)
```
```R
V1 V2 V3 V4
1 chr1 37996770 38005505-
```
```R
#可调整的参数
genefold<-as.numeric("1.5")#放大、缩小展示的范围
# 展示的基因组范围
colnames(gloci)<-c("chr","start","end","strand")
chr<-gloci$chr
gloci$width<-with(gloci,end-start)
startpoint<-gloci$start-genefold*gloci$width
endpoint<-gloci$end+genefold*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,row.names = rownames(bwInfo),stringsAsFactors = F)
for(index in rownames(bwInfo)){
bgFile<-file.path("D:/rtmp/genomeview/bw/",paste(index,".bw",sep=""))
tracklist[]<-DataTrack(range = bgFile,genome="hg38",type="histogram",
name=chartr("_","\n",bwInfo),
col.histogram=coldf)#每个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()
```
页:
[1]