R语言Circos图可视化
准备数据
需要准备两个数据:一个是基因表达谱,另一个是基因的注释(可以为KO注释,也可以是别的什么注释)
基因表达谱
sample1 | sample2 | sample3 | ... | |
---|---|---|---|---|
gene1 | 1.0 | 2.0 | 2.0 | ... |
gene2 | 3.0 | 3.0 | 4.0 | ... |
gene3 | 5.0 | 5.0 | 5.0 | ... |
gene4 | 6.0 | 7.0 | 9.0 | ... |
... | ... | ... | ... | ... |
通路信息
gene | KO | pathway |
---|---|---|
gene1 | KO1 | pathway1 |
gene2 | KO2 | pathway1 |
gene3 | KO2 | pathway2 |
... | ... | ... |
模拟数据
library(tidyverse)
library(magrittr)
library(circlize)
#模拟数据
## Data1
fpkm <- rbind(cbind(matrix(rnorm(500*3, mean = 1), nr = 500), matrix(rnorm(500*3, mean = 2), nr = 500),matrix(rnorm(500*3, mean = 3), nr = 500)))
fpkm <- fpkm[sample(500, 500), ] # randomly permute rows
rownames(fpkm) <- paste0("gene", seq(500))
colnames(fpkm) <- c("A1", "A2", "A3", "B1", "B2", "B3", "C1", "C2", "C3")
fpkm %<>% as.data.frame() %>% mutate(gene = row.names(.))
# Data2
pathways <- rep(paste0("pathway", seq(6)), sample(12:20, size = 6)) %>% sample(70)
KOs <- rep(paste0("KO", seq(20)), sample(5:20, size = 20, replace = TRUE)) %>% sample(70)
KOannotation <- data.frame(KO=KOs, pathway=pathways)
KOannotation <- KOannotation[sample(70, 200, TRUE),]
KOannotation$gene <- sample(paste0("gene", seq(500)),200)# 假设你有富集到的想要可视化的通路
maps <- c("pathway1", "pathway2", "pathway3")
samples <- c("A1", "A2", "A3", "B1", "B2", "B3", "C1", "C2", "C3")
思考画图数据类型
首先简单介绍一下circos对象。正如普通的图有x轴和Y轴,你可以理解为一个circos图有很多个轴(具体数量由你数据确定)。那么每个轴上自然就有对应位置(如下图中的1、2、3、4)。
那么自然而然的很容易想象,如果要画link,需要一个类似这样的数据
from_axis | from_position | to_axis | to_position |
---|---|---|---|
A | 1 | B | 4 |
A | 2 | C | 5 |
A | 3 | D | 6 |
进一步的,如果要控制link的宽度,那应该需要规定每个link的起始和结束的位置
from_axis | from_position_start | from_position_end | to_axis | to_position_start | to_position_end |
---|---|---|---|---|---|
A | 0.5 | 1.5 | B | 3.5 | 4.5 |
A | 1.5 | 2.5 | C | 4.5 | 5.5 |
A | 2.5 | 3.5 | D | 5.5 | 6.5 |
考虑完link之后,考虑热图。首先确认是用circos.heatmap()
画的热图(如上图)。这个数据比较简单,因此考虑先画出热图,再考虑怎么画中间的和弦图。
热图
数据清洗
# 准备热图颜色
col_fun1 = colorRamp2(c(-2, 0, 2), c("#247ab5", "white", "#fda1a0"))
# 需要画图的基因
plot_gene <- KOannotation %>% filter(pathway %in% maps) %>% pull(gene) %>%{filter(fpkm, row.names(fpkm) %in% .)}
# 需要画图的KO
plot_KO <- plot_gene %>% left_join(KOannotation) %>%filter(pathway %in% maps) %>% # There are some unenriched mapgroup_by(KO) %>%summarise(across(samples,sum))
# 需要画图的Pathway
plot_map <- plot_gene %>% left_join(KOannotation) %>%filter(pathway %in% maps) %>% # There are some unenriched mapgroup_by(pathway) %>%summarise(across(samples,sum))plot_data1 <- bind_rows(plot_gene %>% rename(id=gene), plot_KO %>% rename(id=KO), plot_map %>% rename(id=pathway)) %>%`row.names<-`(.$id) %>%select(-id)
plot_data1 <- t(scale(t(plot_data1))) %>% as.data.frame()
# 在热图上划分为gene、KO、pathway
lev_split = row.names(plot_data1) %>% str_match("[a-zA-Z]+") %>% factor()
画图
circos.clear()
circos.par(gap.degree=10, track.height=0.1)
# 分多次画表达谱数据,更有层次
circos.heatmap(plot_data1[samples[1:3]], split = lev_split , col = col_fun1, rownames.side = "outside", cluster = TRUE)
circos.heatmap(plot_data1[samples[7:9]], split = lev_split , col = col_fun1)
circos.heatmap(plot_data1[samples[4:6]], split = lev_split , col = col_fun1)
根据lev_split
,这个热图划分为gene
,KO
,和pathway
三个轴。这里需要指出的是,每个gene、KO、pathway在每个轴上占的长度都是1,例如gene67在gene轴上的位置为0-1,pathway2在pathway轴上的位置为0-1.因此我们要画link的话,根据之前的讨论,如果有三个KO需要连接到pathway2,应该需要一个类似下面的数据:
from_axis | from_position_start | from_position_end | to_axis | to_position_start | to_position_end |
---|---|---|---|---|---|
KO | 0.5 | 1.5 | pathway | 0 | 0.333 |
KO | 1.5 | 2.5 | pathway | 0.333 | 0667 |
KO | 2.5 | 3.5 | pathway | 0.667 | 1 |
注意上表,可能存在多个KO连接到一个pathway的情况,因此我们要对start和end位置进行合理的拆分,以避免重叠。这样做的另一个好处是可以让link线有粗有细,看上去美观许多。所以
另一方面我们要注意,由于circos热图上的基因的是根据聚类结果排布的,所以和数据框里的数据顺序已经不一样了,因此我们首先需要获取画图后的每个gene、KO、pathway在其对应轴上的坐标,这时就需要用到circlize
的get.cell.meta.data
从图上获得相应信息
和弦图
数据清洗
根据上面的讨论结果,我们的数据清洗要实现两个目的:
获取热图生成后的gene、KO、pathway上每个轴的信息,可以整理成如下格式:
id sector position gene1 gene 23 KO1 KO 45 pathway1 pathway 56 表A
将一对多的gene-KO关系、KO-pathway关系进行计算,以获得相对的位置
from_axis from_position_start from_position_end to_axis to_position_start to_position_end gene 0.5 1.5 KO 5 5.33 gene 1.5 2.5 KO 5.33 5.66 gene 2.5 3.5 KO 5.66 6 表B
注意观察上表,由于3个基因连接到了相同的KO,所以我将他们分别连接到了KO不同的位置
进一步,想要获得这个表,可以把过程拆分为以下几步:
2.1生成一个连接对象表
from to gene1 KO2 gene2 KO1 KO1 pathway1 2.2 根据连接对象在表中出现的次数,再计算一个位移
from to from_start from_end to_start to_end gene1 KO2 0 1 0 0.333 gene2 KO1 0 1 0.333 0.667 KO1 pathway1 0 0.5 0.667 1 KO1 pathway2 0.5 1 0 0.0333 表C
2.3 结合表A和表C,计算得到表B
明确思路以后,下面是代码
# 1 获得gene、KO、pathway在每个轴上的位置
plot_data3 = data.frame()
for(lev in levels(lev_split)){a <- rownames(plot_data1)[lev_split==lev][get.cell.meta.data("row_order", sector.index = lev)]a <- seq(length(a)) %>% `names<-`(a) %>% enframe("id", "position")a$sector = levplot_data3 <- rbind(plot_data3, a)
}
plot_data3$position <- plot_data3$position - 1 # 因为每个元素的范围是0~1,所以统一减一方便后面相加
#2.1 获取点对点连接表
plot_data2 <- KOannotation %>%filter(gene %in% row.names(plot_data1), KO %in% row.names(plot_data1)) %>%select(gene, KO) %>%rename(from=gene, to=KO)
tmp_obj1 <- KOannotation %>%filter(KO %in% row.names(plot_data1), pathway %in% row.names(plot_data1)) %>%select(KO, pathway) %>%rename(from=KO, to=pathway)
plot_data2 %<>% bind_rows(tmp_obj1) %>% distinct(from, to)
#2.2 计算位移
tmp_obj1<-plot_data2 %>% group_by(from) %>% mutate(V1=1/n()) %>% group_by(from) %>% mutate(from_end = cumsum(V1), from_start = from_end-V1) %>%select(from, to, from_start, from_end)
tmp_obj2<-plot_data2 %>% group_by(to) %>% mutate(V1=1/n()) %>% group_by(to) %>% mutate(to_end = cumsum(V1), to_start = to_end-V1) %>%select(from, to, to_start, to_end)
plot_data2 %<>% left_join(tmp_obj1) %>% left_join(tmp_obj2)
# 2.3 合并上述两表
plot_data2 %<>% left_join(plot_data3, by=c('from'='id')) %>% rename('from_position'=position, 'from_sector'=sector) %>%left_join(plot_data3, by=c('to'='id')) %>%rename('to_position'=position, 'to_sector'=sector)
由于上述plot_data2
包含了全部点之间的连线关系,我们可能不需要这么多,所以可以挑选自己想要展示的数据。这一步可能需要反复画图几次确定需要展示的数据
#3. 进一步筛选想要展示的连线
# highlight the pathway I wanted
highlight_pathway <- c("pathway1", "pathway2")
highlight_KO <- c("KO5", "KO6", "KO20")
highlight_gene <- c("gene292", "gene256", "gene67", "gene146", "gene52", "gene391", "gene139", "gene327", "gene218", "gene142", "gene375", "gene194")plot_link <- plot_data2 %>% filter(to %in% c(highlight_pathway, highlight_KO)) %>% mutate(col="#dcdcdc80")
highlight_link <- plot_link %>% filter (from %in% c(highlight_gene, highlight_KO)) %>% mutate(col="#fb9b9a80")
plot_link <- plot_link %>% filter (!from %in% c(highlight_gene, highlight_KO))
plot_data2<-rbind(plot_link, highlight_link)
使用circos.link
循环作图
for ( idx in seq(nrow(plot_data2))){tmp_obj <- plot_data2[idx,]circos.link(tmp_obj[['from_sector']], c(tmp_obj[['from_position']] + tmp_obj[['from_start']], tmp_obj[['from_position']] + tmp_obj[['from_end']]), tmp_obj[['to_sector']], c(tmp_obj[['to_position']] + tmp_obj[['to_start']], tmp_obj[['to_position']] + tmp_obj[['to_end']]),col = tmp_obj[['col']],border = NA)
}
全部代码
library(tidyverse)
library(magrittr)
library(circlize)
#模拟数据
## Data1
fpkm <- rbind(cbind(matrix(rnorm(500*3, mean = 1), nr = 500), matrix(rnorm(500*3, mean = 2), nr = 500),matrix(rnorm(500*3, mean = 3), nr = 500)))
fpkm <- fpkm[sample(500, 500), ] # randomly permute rows
rownames(fpkm) <- paste0("gene", seq(500))
colnames(fpkm) <- c("A1", "A2", "A3", "B1", "B2", "B3", "C1", "C2", "C3")
fpkm %<>% as.data.frame() %>% mutate(gene = row.names(.))
# Data2
pathways <- rep(paste0("pathway", seq(6)), sample(12:20, size = 6)) %>% sample(70)
KOs <- rep(paste0("KO", seq(20)), sample(5:20, size = 20, replace = TRUE)) %>% sample(70)
KOannotation <- data.frame(KO=KOs, pathway=pathways)
KOannotation <- KOannotation[sample(70, 200, TRUE),]
KOannotation$gene <- sample(paste0("gene", seq(500)),200)# 假设你有富集到的想要可视化的通路
maps <- c("pathway1", "pathway2", "pathway3")
samples <- c("A1", "A2", "A3", "B1", "B2", "B3", "C1", "C2", "C3")
# 准备热图颜色
col_fun1 = colorRamp2(c(-2, 0, 2), c("#247ab5", "white", "#fda1a0"))
# 需要画图的基因
plot_gene <- KOannotation %>% filter(pathway %in% maps) %>% pull(gene) %>%{filter(fpkm, row.names(fpkm) %in% .)}
# 需要画图的KO
plot_KO <- plot_gene %>% left_join(KOannotation) %>%filter(pathway %in% maps) %>% # There are some unenriched mapgroup_by(KO) %>%summarise(across(samples,sum))
# 需要画图的Pathway
plot_map <- plot_gene %>% left_join(KOannotation) %>%filter(pathway %in% maps) %>% # There are some unenriched mapgroup_by(pathway) %>%summarise(across(samples,sum))plot_data1 <- bind_rows(plot_gene %>% rename(id=gene), plot_KO %>% rename(id=KO), plot_map %>% rename(id=pathway)) %>%`row.names<-`(.$id) %>%select(-id)
plot_data1 <- t(scale(t(plot_data1))) %>% as.data.frame()
# 在热图上划分为gene、KO、pathway
lev_split = row.names(plot_data1) %>% str_match("[a-zA-Z]+") %>% factor()
circos.clear()
circos.par(gap.degree=10, track.height=0.1)
# 分多次画表达谱数据,更有层次
circos.heatmap(plot_data1[samples[1:3]], split = lev_split , col = col_fun1, rownames.side = "outside", cluster = TRUE)
circos.heatmap(plot_data1[samples[7:9]], split = lev_split , col = col_fun1)
circos.heatmap(plot_data1[samples[4:6]], split = lev_split , col = col_fun1)
plot_data3 = data.frame()
for(lev in levels(lev_split)){a <- rownames(plot_data1)[lev_split==lev][get.cell.meta.data("row_order", sector.index = lev)]a <- seq(length(a)) %>% `names<-`(a) %>% enframe("id", "position")a$sector = levplot_data3 <- rbind(plot_data3, a)
}
plot_data3$position <- plot_data3$position - 1 # 因为每个元素的范围是0~1,所以统一减一方便后面相加
#2.1 获取点对点连接表
plot_data2 <- KOannotation %>%filter(gene %in% row.names(plot_data1), KO %in% row.names(plot_data1)) %>%select(gene, KO) %>%rename(from=gene, to=KO)
tmp_obj1 <- KOannotation %>%filter(KO %in% row.names(plot_data1), pathway %in% row.names(plot_data1)) %>%select(KO, pathway) %>%rename(from=KO, to=pathway)
plot_data2 %<>% bind_rows(tmp_obj1) %>% distinct(from, to)
#2.2 计算位移
tmp_obj1<-plot_data2 %>% group_by(from) %>% mutate(V1=1/n()) %>% group_by(from) %>% mutate(from_end = cumsum(V1), from_start = from_end-V1) %>%select(from, to, from_start, from_end)
tmp_obj2<-plot_data2 %>% group_by(to) %>% mutate(V1=1/n()) %>% group_by(to) %>% mutate(to_end = cumsum(V1), to_start = to_end-V1) %>%select(from, to, to_start, to_end)
plot_data2 %<>% left_join(tmp_obj1) %>% left_join(tmp_obj2)
# 2.3 合并上述两表
plot_data2 %<>% left_join(plot_data3, by=c('from'='id')) %>% rename('from_position'=position, 'from_sector'=sector) %>%left_join(plot_data3, by=c('to'='id')) %>%rename('to_position'=position, 'to_sector'=sector)
highlight_pathway <- c("pathway1", "pathway2")
highlight_KO <- c("KO5", "KO6", "KO20")
highlight_gene <- c("gene292", "gene256", "gene67", "gene146", "gene52", "gene391", "gene139", "gene327", "gene218", "gene142", "gene375", "gene194")plot_link <- plot_data2 %>% filter(to %in% c(highlight_pathway, highlight_KO)) %>% mutate(col="#dcdcdc80")
highlight_link <- plot_link %>% filter (from %in% c(highlight_gene, highlight_KO)) %>% mutate(col="#fb9b9a80")
plot_link <- plot_link %>% filter (!from %in% c(highlight_gene, highlight_KO))
plot_data2<-rbind(plot_link, highlight_link)
for ( idx in seq(nrow(plot_data2))){tmp_obj <- plot_data2[idx,]circos.link(tmp_obj[['from_sector']], c(tmp_obj[['from_position']] + tmp_obj[['from_start']], tmp_obj[['from_position']] + tmp_obj[['from_end']]), tmp_obj[['to_sector']], c(tmp_obj[['to_position']] + tmp_obj[['to_start']], tmp_obj[['to_position']] + tmp_obj[['to_end']]),col = tmp_obj[['col']],border = NA)
}
成品
注意点
如果画多层圈图,会根据行号对齐,但是不会根据行名对齐。所以所有组的表达谱数据必须一步准备。
热图和和弦图代码之间进行了大量数据清洗和计算,不要怕,没关系。因为只有拿到热图之后才能获取对应gene、KO、pathway的位置
每次画图之前习惯性的用
circos.clear()
清理一下缓存。
R语言Circos图可视化相关推荐
- R语言配对图可视化:配对图(pair plot)可视化(根据分类变量的值为散点图上的数据点添加颜色和形状、Add color and shape by variables)
R语言配对图可视化:配对图(pair plot)可视化(根据分类变量的值为散点图上的数据点添加颜色和形状.Add color and shape by variables) 目录
- R语言配对图可视化:pivot_longer函数将宽格式的数据重塑为长格式并进行数据全连接和左连接(left join)、配对图可视化(根据分类变量的值为散点图上的数据点添加颜色)
R语言配对图可视化:pivot_longer函数将宽格式的数据重塑为长格式并进行数据全连接和左连接(left join).配对图可视化(根据分类变量的值为散点图上的数据点添加颜色,Add color ...
- R语言plotly包可视化线图(line plot)、使用restyle参数自定义设置可视化结果中线条的颜色、使用按钮动态切换线条的颜色(change line color with button)
R语言plotly包可视化线图(line plot).使用restyle参数自定义设置可视化结果中线条的颜色.使用按钮动态切换线条的颜色(change line color with button i ...
- R语言使用ggplot2可视化凹凸图(bumps chart、凹凸图是一种特殊形式的线图,旨在探索随着时间的推移等级的变化)、并设置凹凸图的线条为曲线而不是直线(change into curves)
R语言使用ggplot2可视化凹凸图(bumps chart.凹凸图是一种特殊形式的线图,旨在探索随着时间的推移等级的变化).并设置凹凸图的线条为曲线而不是直线(change bumps chart ...
- R语言ggridges包可视化山脊图(Ridgeline Plots)并且在山脊图中添加均值竖线(Add Mean Line to RIdgeline Plot with ggridges in R)
R语言ggridges包可视化山脊图(Ridgeline Plots)并且在山脊图中添加均值竖线(Add Mean Line to RIdgeline Plot with ggridges in R) ...
- R语言sunburst图(sunburst plot)可视化实战:使用sunburstR包和ggplot2包进行可视化
R语言sunburst图(sunburst plot)可视化实战:使用sunburstR包和ggplot2包进行可视化 目录 R语言sunburst图
- R语言删除ggplot可视化图中的所有x轴轴标签实战:ggplot可视化默认包含所有x轴轴标签、删除ggplot可视化图中的所有x轴轴标签实战
R语言删除ggplot可视化图中的所有x轴轴标签实战:ggplot可视化默认包含所有x轴轴标签.删除ggplot可视化图中的所有x轴轴标签实战 目录
- R语言相关关系可视化函数梳理(附代码)
来源:R语言中文社区 作者:赵镇宁 本文约3177字,建议阅读6分钟. 本文为你介绍R语言相关关系可视化的函数进行了初步梳理,大家可根据个人需求及函数功能择优选择. 当考察多个变量间的相关关系时,通常 ...
- 送书|北大出版:R语言数据分析与可视化从入门到精通
生物信息学习的正确姿势 NGS系列文章包括NGS基础.高颜值在线绘图和分析.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流 ...
最新文章
- Python常用模块之configparser模块
- WebAssembly 系列(五)为什么 WebAssembly 更快? 1
- 截取图片生成头像插件
- tomcat(15)Digester库
- eclipse把tomcant用到一个项目里_聊一个镜头工艺里容易被忽略,但很重要的项目...
- html5显示字母的值,使用HTML5 Canvas API控制字体的显示与渲染的方法
- ELK6.0部署:Elasticsearch+Logstash+Kibana搭建分布式日志平台
- 洛谷1005 【NOIP2007】矩阵取数游戏
- 牛腩视频播放管理系统
- Crackme015
- 计算机网络专业学python_「非计算机专业」小白如何学好Python?
- 显示纯服务器_BBT三行代码搭建服务器,让Dynamo跳出IronPython的封锁
- c语言程序设计黄保和第二章,c语言程序设计答案(选择题+编程)黄保和、江戈版...
- 常见音频编码格式解析
- 李运华老师的一些经典见解收藏
- RFID中的天线技术-应用及设计现状
- cobar mysql 性能,Cobar + MySQL 技術驗證(li)
- 大数据平台架构有哪些
- 解决创建MAVEN工程速度慢的问题
- docker命令push,pull等设置代理
热门文章
- tensorflow conv2d()参数解析
- JS 复习(6)JavaScript对象
- Nexmo 短信平台接口 遇到的坑
- 升级安装win11 22H2(跳过TPM和CPU等检测)
- Oracle Forensics t00ls
- DDR2 sodimm + Flash + Triple-Speed Ethernet + IO in nios
- c++中数字与字符,字符与其ASCII转换
- SpringCloudSpringBoot集成Acivity6.0
- 2020 Python中文社区热门文章 Top 10
- Python27 No module named PIL解决方法