R语言系统发育树专栏 | 复杂系统发育树教程 (二)

文摘   2024-12-17 18:35   云南  

一边学习,一边总结,一边分享!

由于微信改版,一直有同学反映。存在长时间接收不到公众号的推文。那么请跟随以下步骤,将小杜的生信筆記设置为星标,不错过每一条推文教程。

欢迎关注《小杜的生信笔记》!!

如何加入社群

小杜的生信笔记仅有微信社群

1. 微信群:付费社群。添加小杜好友,加友请知:加友须知!!,加入社群请查看小杜的生信笔记付费加友入群声明

入群声明

2. 小杜个人微信:若你有好的教程或想法,可添加小杜个人微信。值得注意的是,小杜个人微信并不支持免费咨询长时间咨询,但支持小问题2-3个免费咨询。

小杜微信:

知识星球:

本期教程图形

在后台回复关键词获得代码:20241217

关于的《R语言系统发育树专栏》

基于R语言绘制系统发育树,关于这块的内容自己一直想学习(PS:确实是学习),在自己做的这块中很少使用到类似的教程。但对于自己来说,不是一直不适用,只是使用的概率比较少而已。我想很多的同学也会存在与我自己一样的情况,因此,自己一边学习,一边记录,一边分享吧。对于使用ggtree、ggtreeExtra包,或是使用其他可视化软件绘制进化树的教程,在网上会有很多的教程(PS:只要你愿意搜索)。我们在自己学习过程中,也会引用各位博主大佬分享的教程脚本,大家一起学习。


Code

rm(list=ls())
#!加载R包
# BiocManager::install("ggtree")
# BiocManager::install("ggtreeExtra")
# BiocManager::install("treeio")
# install.packages("tidytree")
# install.packages("ggstar")
# install.packages("ggplot2")
# install.packages("ggnewscale")
######################1、加载包####################
library(ggtreeExtra)
library(ggtree)
library(treeio)
library(tidytree)
library(ggstar)
library(ggplot2)
library(ggnewscale)
  1. 加载数据
# 读取树形结构文件 (Newick格式)
tree <- read.tree("data/kegg.nwk") # 读取 Newick 格式的树形数据
# 读取各环的属性数据
dt1 <- read.csv("data/tippoint_attr.csv") # 读取提示点的属性数据
dt2 <- read.csv("data/firstring_attr.csv") # 读取第一环的属性数据
dt3 <- read.csv("data/secondring_attr.csv") # 读取第二环的属性数据
dt4 <- read.csv("data/barplot_attr.csv") # 读取条形图的属性数据
head(dt1)
> head(dt1)
ID Phyla
1 aae Aquificae
2 aap Spirochaetes/Proteobacteria
3 aas Bacteroidetes
4 aba Other
5 abc Spirochaetes/Proteobacteria
6 abi Euryarchaeota
> head(dt2)
ID ring Type1
1 pto 1 V/A-type ATPase
2 vmo 1 V/A-type ATPase
3 bvu 1 V/A-type ATPase
4 osp 1 V/A-type ATPase
5 hje 1 V/A-type ATPase
6 ttj 1 V/A-type ATPase
> head(dt3)
ID Abundance Type2
1 ecp 0.5000000 acyl-CoA synth
2 ica 0.3333333 acyl-CoA synth
3 bge 0.6153846 acyl-CoA synth
4 bbi 0.1666667 acyl-CoA synth
5 cab 0.4166667 acyl-CoA synth
6 ate 0.3333333 acyl-CoA synth
> head(dt4)
ID Length
1 pmm 0.3909040
2 syg 0.6531337
3 cyb 0.6526900
4 cya 0.6427066
5 syc 0.5728231
6 sye 0.5229063
  1. 数据处理
##'@数据处理
dt1 <- aggregate(.~ID, dt1, paste0, collapse="/")
##'@重新排列门类列的顺序
dt1$Phyla <- factor(dt1$Phyla, levels=c("Actinobacteria","Aquificae","Bacteroidetes",
"Chlamydiae","Chlorobi","Chloroflexi","Crenarchaeota",
"Cyanobacteria","Euryarchaeota","Firmicutes",
"Spirochaetes/Proteobacteria",
"Tenericutes","Thermi","Thermotogae","Other"))

#'@对Type2重新排序
dt3$Type2 <- factor(dt3$Type2, levels=c("FA synth init", "FA synth elong",
"acyl-CoA synth", "beta-Oxidation",
"Ketone biosynth"))
#'@提取支系层的节点标签

nodelab <- tree$node.label[nchar(tree$node.label)>0]
nodeids <- nodeid(tree, nodelab)

#'@支系标签的位置
textex <- c(1.0, 0.4, 0.2, 1.4, 1.4, 0.4, 1.4, 1.4, 0.4, 0.4,
0.8, 1, 0.6, 0.6, 0.4, 0.3, 0, 0.4, 0.1, 0.25,
0.2, 0.3, 0.8, 0.8, 0.8, 0.6, 2.4)

#‘@设置层级标签的属性
cladelabels <- mapply(function(x, y, z){geom_cladelabel(node=x, label=y, barsize=NA, extend=0.3,
offset.text=z, fontsize=1.3, angle="auto",
hjust=0.5, horizontal=FALSE, fontface="italic")},
nodeids, nodelab, textex, SIMPLIFY=FALSE)


##'@设置高亮图层的颜色
fills <- c("#808080", "#808080", "#808080", "#808080", "#808080",
"#191970", "#87CEFA", "#FFC125", "#B0171F", "#B0171F",
"#B0171F", "#B0171F", "#B0171F", "#B0171F", "#B0171F",
"#B0171F", "#B0171F", "#B0171F", "#B0171F", "#B0171F",
"#B0171F", "#B0171F", "#9ACD32", "#9ACD32", "#9ACD32",
"#006400", "#800000")

##'@设置节点的高亮
highlights <- mapply(function(x, y){geom_hilight(node=x, extendto=5.8, alpha=0.3,
fill=y, color=y, size=0.05)},
nodeids, fills, SIMPLIFY=FALSE)


##'@自定义颜色
colors <- c("#9ACD32", "#EE6A50", "#87CEFA", "#FFC125", "#D15FEE", "#8DEEEE", "#800000",
"#006400", "#800080", "#808080", "#B0171F", "#191970", "#7B68EE",
"#00CD00", "Black")

  1. 基础图形绘制
p1 <- ggtree(
tree,
layout="circular",
size=0.1
)

p1

ggsave("./output/fig1.jpg",width = 6, height = 6)
##'@设置高亮
p1 <- p1 +
highlights

p1
p2 <- p1 %<+% dt1 +
geom_tippoint(
mapping=aes(
fill=Phyla # 将数据列 "Phyla" 映射为点的填充颜色
),
shape = 21, # 点的形状,此处为带有边框的圆形
size = 1.2, # 点的大小
stroke = 0.05, # 点边框的宽度
position = "identity", # 直接绘制,无额外调整位置
show.legend = FALSE # 不在图例中显示这些点
)+
scale_fill_manual(values=colors) # 手动设置点填充颜色的映射,`colors` 是一个向量

p2

ggsave("./output/fig3.jpg",width = 6, height = 6)

添加额外的分支注释标签和新填充映射

p2 <- p2 +
cladelabels +
new_scale_fill()

p2

ggsave("./output/fig3-2.jpg", plot = p2, width = 6, height = 6)
p3 <- p2 + 
geom_fruit( # 在树旁添加图形(如热图)
data = dt2, # 关联的数据
geom = geom_tile, # 使用矩形块绘制热图
mapping = aes( # 映射数据到图形属性
y = ID, # 系统树节点 ID
x = ring, # 热图的列分组变量
fill = Type1 # 热图填充颜色的变量
),
offset = 0.01, # 热图与树的水平偏移距离
pwidth = 0.14 # 热图宽度占比
) +
scale_fill_manual( # 手动设置热图填充颜色
name = "ATP synthesis", # 图例标题
values = c("#339933", "#dfac03"), # 自定义颜色
guide = guide_legend( # 图例样式
keywidth = 0.35, # 图例键宽度
keyheight = 0.35, # 图例键高度
order = 1 # 图例显示顺序
)
) +
new_scale_fill() # 再次引入新的填充映射

p3

ggsave("./output/fig4.jpg", plot = p3, width = 6, height = 6)

添加热图

p4 <- p3 + 
geom_fruit(
data = dt3, # 第二个数据集,包含系统树叶节点的相关信息
geom = geom_tile, # 热图样式,使用矩形块表示
mapping = aes( # 数据映射
y = ID, # 系统树节点的 ID
alpha = Abundance, # 丰度值映射为透明度
x = Type2, # 分类变量,用于热图的列分组
fill = Type2 # 分类变量,用于填充颜色
),
offset = 0.001, # 热图与树的距离
pwidth = 0.18 # 热图宽度占比
) +
scale_fill_manual( # 手动设置分类变量的填充颜色
name = "Fatty Acid metabolism", # 图例标题
values = c("#b22222", "#005500", "#0000be", "#9f1f9f", "#793a07"), # 颜色
guide = guide_legend( # 自定义图例样式
keywidth = 0.35,
keyheight = 0.35,
order = 2
)
) +
scale_alpha_continuous( # 透明度映射
range = c(0, 0.4), # 透明度范围
guide = guide_legend( # 自定义透明度图例样式
keywidth = 0.35,
keyheight = 0.35,
order = 3
)
) +
new_scale_fill() # 新的填充映射,用于后续操作


p4

ggsave("./output/fig5.jpg", width = 7, height = 6)

添加柱状图

p5 <- p4 +
geom_fruit(
data = dt4, # 数据集,包含柱状图信息
geom = geom_bar, # 添加柱状图
mapping = aes(
y = ID, # 与树节点对应
x = Length, # 表示定量变量
fill = Phyla # 分类变量,用颜色表示
),
stat = "identity", # 数据按原值绘制
orientation = "y", # 水平方向柱状图
pwidth = 0.3, # 柱状图宽度
position = position_dodgex() # 横向分组排列
) +
scale_fill_manual(
values = colors, # 设置分类变量颜色
guide = guide_legend(
keywidth = 0.35, # 图例键宽度
keyheight = 0.35, # 图例键高度
order = 4 # 图例显示顺序
)
) +
geom_treescale(
fontsize = 1.2, # 比例尺文字大小
linesize = 0.3 # 比例尺线条粗细
) +
theme(
legend.position = c(0.93, 0.76), # 图例位置
legend.background = element_rect(fill = NA), # 图例背景透明
legend.title = element_text(size = 6), # 图例标题字体大小
legend.text = element_text(size = 4.5), # 图例内容字体大小
legend.spacing.y = unit(0.02, "cm") # 图例项之间的垂直间距
)

p5

ggsave("./output/fig6.jpg", plot = p5, width = 8, height = 6)

在后台回复关键词获得代码:20241217


若我们的教程对你有所帮助,请点赞+收藏+转发,大家的支持是我们更新的动力!!

往期部分文章

1. 最全WGCNA教程(替换数据即可出全部结果与图形)

推荐大家购买最新的教程,若是已经购买以前WGNCA教程的同学,可以在对应教程留言,即可获得最新的教程。(注:此教程也仅基于自己理解,不仅局限于此,难免有不恰当地方,请结合自己需求,进行改动。)


2. 精美图形绘制教程

3. 转录组分析教程

4. 转录组下游分析

小杜的生信筆記 ,主要发表或收录生物信息学教程,以及基于R分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!

小杜的生信筆記
小杜的生信筆記,生信小白,初来乍到请多指教。 主要学习分享,转录组数据分析,基于R语言数据分析和绘制图片等,以及相关文献的分享。
 最新文章