最近在看論文
論文中的部分代碼是公開的，代碼的鏈接是
https://github.com/CornilleAmandine/-apricot_evolutionary_history_2021
其中有一個(gè)畫弦圖的代碼
正好自己最近在學(xué)習(xí)circlize這個(gè)包，所以重復(fù)一下這個(gè)代碼
但是這個(gè)代碼只有一部分，數(shù)據(jù)也只公開了染色體長(zhǎng)度的部分，所以我們只能按照這個(gè)代碼畫出最外圈表示染色體的部分，也就是論文中Figure6 a 的最外圈
以下是代碼
首先是讀入染色體的長(zhǎng)度
ref<-read.table("circlize/Genome_len.chr",header = TRUE)
初始的一些參數(shù)設(shè)置
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))
畫表示染色體的矩形塊
這里我把顏色改動(dòng)了一下，我個(gè)人認(rèn)為這個(gè)原始論文中有點(diǎn)偏 屎黃 的配色不太好看
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)
如果想要實(shí)現(xiàn)內(nèi)圈的內(nèi)容 可以參考 小白魚的微信推文 https://mp.weixin.qq.com/s/KY9IZ91YYLNNXasJh2E2Ug
介紹的很詳細(xì)了
我按照這個(gè)推文模仿了基因密度，如何統(tǒng)計(jì)基因密度 可以參考推文
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()
最終出圖
這個(gè)示例數(shù)據(jù)還不能公開
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1、R語(yǔ)言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡(jiǎn)單小例子；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！
