How to plot chord diagram to show gene density with R language circlize package

created at 07-13-2021 views: 25

greate example on Cell

I was reading a paper recently

paper

Part of the code in the paper is public, and the link to the code is

https://github.com/CornilleAmandine/-apricot_evolutionary_history_2021

github project

There is a code for drawing a chord diagram

I happen to be learning the circlize package recently, so repeat this code

But this code is only part of the code, and the data only discloses the length of the chromosome, so we can only draw the outermost circle representing the part of the chromosome according to this code, which is the outermost circle of Figure 6 a in the paper.

chord diagram

code

The first is the length of the read chromosome

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

initial parameter settings

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

rectangular blocks representing chromosomes

Here I changed the color. I personally think that the color scheme of the original paper is a little bit too yellow

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)

rectangular blocks representing chromosomes

Add the outermost scale

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)

outermost scale

complete code

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

final result:

final result

created at:07-13-2021
edited at: 07-13-2021: