画出igv款式的矢量图

生信喵实验柴

发布于 2023-09-04 08:06:10

3950

发布于 2023-09-04 08:06:10

文章被收录于专栏：生信喵实验柴

背景

展示感兴趣的基因附近的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

开始画图

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()