跟着Nature Communications学画图：R语言circlize包画弦图展示基因密度

用户7010445

发布于 2021-07-30 10:42:34

1.4K0

发布于 2021-07-30 10:42:34

以下是代码

首先是读入染色体的长度

ref<-read.table("circlize/Genome_len.chr",header = TRUE)

初始的一些参数设置

library(circlize)
circos.clear()
col_text <- "grey20"
circos.par("track.height"=0.8,gap.degree=5,start.degree =86,clock.wise = T,
           cell.padding=c(0,0,0,0))
circos.initialize(factors=ref$Genome,          
                  xlim=matrix(c(rep(0,8),ref$Length),ncol=2))

画表示染色体的矩形块

这里我把颜色改动了一下，我个人认为这个原始论文中有点偏屎黄的配色不太好看

circos.track(ylim=c(0,1),panel.fun=function(x,y) {
  Genome=CELL_META$sector.index
  xlim=CELL_META$xlim
  ylim=CELL_META$ylim
  circos.text(mean(xlim),mean(ylim),Genome,cex=0.5,col=col_text,
              facing="bending.inside",niceFacing=TRUE)
},bg.col="#00ADFF",bg.border=F,track.height=0.06)

添加最外圈的刻度

brk <- c(0,0.5,1,1.5,2,2.5,3,3.5,4,4.5)*10^7
circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {
  circos.axis(h="top",major.at=brk,labels=round(brk/10^7,1),labels.cex=0.4,
              col=col_text,labels.col=col_text,lwd=0.7,labels.facing="clockwise")
},bg.border=F)

如果想要实现内圈的内容可以参考

https://mp.weixin.qq.com/s/KY9IZ91YYLNNXasJh2E2Ug

介绍的很详细了

我按照这个推文模仿了基因密度，如何统计基因密度可以参考推文

使用Tbtools根据gtf文件统计基因密度

代码

library(ComplexHeatmap)

library(circlize)
col_text <- "grey40"
lncRNA_density<-read.csv("fruit_ripening/data/gene_density/lncRNA_gene_density.tsv",
                         sep="\t",header = F) %>% 
  arrange(V1,V2)
head(lncRNA_density)
summary(lncRNA_density$V4)

mRNA_density<-read.csv("fruit_ripening/data/gene_density/mRNA_gene_density.tsv",
                       header=F,sep="\t") %>% 
  arrange(V1,V2)

head(mRNA_density)
summary(mRNA_density$V4)

color_assign <- colorRamp2(breaks = c(1, 10, 21), 
                           col = c('#00ADFF', 'orange', 'green2'))

chr<-read.csv("fruit_ripening/data/gene_density/chr_len.txt",
              header=F,sep="\t")
chr
circos.par("track.height"=0.8,gap.degree=5,cell.padding=c(0,0,0,0))
circos.initialize(factors=chr$V1,
                  xlim=matrix(c(rep(0,8),chr$V2),ncol=2))

circos.track(ylim=c(0,1),panel.fun=function(x,y) {
  chr=CELL_META$sector.index
  xlim=CELL_META$xlim
  ylim=CELL_META$ylim
  circos.text(mean(xlim),mean(ylim),chr,cex=0.5,col=col_text,
              facing="bending.inside",niceFacing=TRUE)
},bg.col="grey90",bg.border=F,track.height=0.06)
brk <- c(0,10,20,30,40,50,55)*1000000
brk_label<-paste0(c(0,10,20,30,40,50,55),"M")
circos.track(track.index = get.current.track.index(), 
             panel.fun = function(x, y) {
               circos.axis(h="top",
                           major.at=brk,
                           labels=brk_label,
                           labels.cex=0.4,
                           col=col_text,
                           labels.col=col_text,
                           lwd=0.7,
                           labels.facing="clockwise")
             },
             bg.border=F)

circos.genomicTrackPlotRegion(
  lncRNA_density, track.height = 0.12, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
  } )


circos.genomicTrackPlotRegion(
  mRNA_density, track.height = 0.12, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
  } )

gene_legend <- Legend(
  at = c(1, 10, 21), 
  labels = c(1,10,21),
  labels_gp = gpar(fontsize = 8),
  col_fun = color_assign,
  title = 'gene density', 
  title_gp = gpar(fontsize = 9), 
  grid_height = unit(0.4, 'cm'), 
  grid_width = unit(0.4, 'cm'), 
  type = 'points', pch = NA, 
  background = c('#00ADFF', 'orange', 'green2'))

pushViewport(viewport(x = 0.5, y = 0.5))
grid.draw(gene_legend)

upViewport()

circos.clear()

最终出图