跟着Nature Communications学画图：R语言ci
2021-07-11 本文已影响0人
小明的数据分析笔记本
最近在看论文
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
介绍的很详细了
我按照这个推文模仿了基因密度，如何统计基因密度 可以参考推文
https://mp.weixin.qq.com/s/KerMMCTjWzso4yKq_zzNew
代码
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()
最终出图
image.png
这个示例数据还不能公开
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
