/ 小明的数据分析笔记本_跟着NatureCommunications学画图:R语言circlize包画弦图展示基因密度

小明的数据分析笔记本_跟着NatureCommunications学画图:R语言circlize包画弦图展示基因密度

由Punicagranatum发表于2021-07-11 21:42:04

最近在看论文

image.png

论文中的部分代码是公开的,代码的链接是

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

image.png

其中有一个画弦图的代码

正好自己最近在学习circlize这个包,所以重复一下这个代码

但是这个代码只有一部分,数据也只公开了染色体长度的部分,所以我们只能按照这个代码画出最外圈表示染色体的部分,也就是论文中Figure6 a 的最外圈

image.png

以下是代码

首先是读入染色体的长度
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)
image.png
添加最外圈的刻度
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)

image.png

如果想要实现内圈的内容 可以参考 小白鱼的微信推文 

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=" ",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=" ") %>% 
  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=" ")
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()

最终出图

image.png

这个示例数据还不能公开

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享:1、R语言和python做数据分析和数据可视化的简单小例子;2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记;3、生物信息学入门学习资料及自己的学习笔记!


作者:Punicagranatum

公众号:小明的数据分析笔记本

原文链接:http://mp.weixin.qq.com/s?__biz=MzI3NzQ3MTcxMg==&mid=2247489530&idx=1&sn=e073e3df5a79460f29a99770537e491f&chksm=eb649f75dc13166340bfe72b834c8b14ac8f8507a4a4d125109271e9bf283ea5b9e83f569496&scene=0&xtrack=1#rd

发布时间:2021-07-11 21:42:04

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注