R语言ggplot2 | 三元图

Posted 2023-04-06 酷在前行

📋文章目录

• R语言绘图过程
• • R包及数据加载
  • 计算绘图数值
  • 数据可视化

   三元图，顾名思义，是一个等边三角形式的图像，它将本该是三维的x,y,z三轴转化为二维的三角形展示出来，三角形的三个角可以是一个或者一组样本，通过观察三角形中点的位置判断样本在三组间的分布状况。通常这类图用于展示组学数据（进行差异表达基因的分析），作为延伸，也可以用来分析微生物16S/ITS/18S扩增子宏基因组数据，并且可以结合火山图分析，探究微生物在三组样本间的富集状况。

R语言绘图过程

R包及数据加载

   本例中使用的R包为ggtern和limma（用于差异矩阵的构建与三元图绘图）

# load packages
packages = c("ggtern", "limma")
lapply(packages, function(x) if (!(x %in% installed.packages())) 
  install.packages(x)
)
lapply(packages, library, quietly = TRUE, character.only = TRUE)

例子数据情况如下：800行，105列。其中分为四组，包括ck、om、pp、others。行为otu id，列为sample id。

dori<- read.csv("volcano figure.csv", header = T, row.names = 1);head(dori)
group<- data.frame(site= colnames(dori),#构建测试分组
                   group=factor(c(rep('others',30), rep('ck',20),
                                  rep('om',29), rep('pp',26))))
data<- dori[, which(names(dori) %in% filter(group, group != 'others')[,1])]; head(data)#筛选用于分析绘图的数据
str(data); table(is.na(data))#检查数据格式是否正确，是否具有缺失值

核查数据，确保无误（必须要求为数值型向量）：

计算绘图数值

数据加载好后，我们需要计算

1. 每个otu在总表里的平均丰度，作为对应点（气泡）的大小。#size
2. 剔除others组后，分别计算ck、om、pp小组中各otu在三组总表中的丰度，作为对应点在三元图中的定位（x, y, z）。

sum<- data.frame(tsum= rowSums(dori),
                  dsum= rowSums(data))#计算总表otu和与ck,om,pp三表otu总和
#ck
cdata<- data[,which(colnames(data) %in%
                      filter(group, group== "ck")[,1])]
csum<- data.frame(site= rowSums(cdata) / sum$dsum, #位置参数
                  size= rowSums(cdata) / sum$tsum)#气泡大小
#om
odata<- data[,which(colnames(data) %in%
                      filter(group, group== "om")[,1])]
osum<- data.frame(site= rowSums(odata) / sum$dsum,
                  size= rowSums(odata) / sum$tsum)
#pp
pdata<- data[,which(colnames(data) %in%
                      filter(group, group== "pp")[,1])]
psum<- data.frame(site= rowSums(pdata) / sum$dsum,
                  size= rowSums(pdata) / sum$tsum)

合并气泡位置数据与大小数据：

site<- data.frame(otu= rownames(data), ck= csum$site, om= osum$site, pp= psum$site)
size<- data.frame(csum$size, osum$size, psum$size) %>%
  apply(1, max)
pdata<- data.frame(site, size); pdata[is.na(pdata)]<- 0; head(pdata)

构建差异分析矩阵并进行差异分析，寻找两组数据间的差异丰度微生物类群

design <- model.matrix(~0+factor(group$group[31:105]))
colnames(design) = levels(factor(group$group[31:105]))
rownames(design) = colnames(data);design#设计对比矩阵模型
#由于rep函数是有规律性的，实际构建design矩阵时要确保otu表里的样品名与group分组对应

contrast.matrix1<- makeContrasts(ck-pp, ck-om, levels = design)
contrast.matrix2<- makeContrasts(om-ck, om-pp, levels = design)
contrast.matrix3<- makeContrasts(pp-ck, pp-om, levels = design)#分别就三个组别进行差异分析

构建模型（ck）：

fit<- lmFit(data,design)#构建模型
fit1<- contrasts.fit(fit,contrast.matrix1)
fit2<- eBayes(fit1)
tempOutput = topTable(fit2,coef = 1,n = Inf)
nrDEG = na.omit(tempOutput)
nrDEG$Change<- ifelse(nrDEG$logFC>=2&nrDEG$P.Value<=0.05,"Up",
                           ifelse(nrDEG$logFC<2&nrDEG$P.Value<=0.05,"Down","Not Sig"))
ck<- nrDEG %>%
  filter(Change %in% c("Up", "Not Sig")) %>%
  select(Change)
ck<- mutate(ck, con= "3",
           otu= rownames(ck))
count(ck, ck$change)#使用富集出来的up与Not sig进行绘图，只给Up着色

对ck，om，pp逐个构建模型太过繁琐，可以使用循环简化这个过程：

clist<- list(contrast.matrix1, contrast.matrix2, contrast.matrix3)
rlist<- list()
for(i in 1:length(clist))
  nrDEG <- contrasts.fit(fit,clist[[i]]) %>%
    eBayes() %>% topTable(coef = 1,n = Inf) %>% na.omit()
  nrDEG$Change <- ifelse(nrDEG$logFC>=2&nrDEG$P.Value<=0.05,"Up",
                         ifelse(nrDEG$logFC<2&nrDEG$P.Value<=0.05,"Down",
                                "Not Sig"))
  res<- nrDEG %>%
    filter(Change %in% c("Up", "Not Sig")) %>%
    select(Change)
  res<- mutate(res, con= 4-i,
              otu= rownames(res))
  rlist[[i]]<- res

数据可视化

筛选出Up和Not sig分组进行可视化：

Up<- rbind(filter(rlist[[1]], Change == "Up"),
           filter(rlist[[2]], Change == "Up"),
           filter(rlist[[3]], Change == "Up")) %>%#筛选出Up组，准备后续着色
  group_by(otu) %>% filter(con == max(con)) #将重合的少数otu划归到ck中

No<- rbind(filter(rlist[[1]], Change == "Not Sig"),
           filter(rlist[[2]], Change == "Not Sig"),
           filter(rlist[[3]], Change == "Not Sig")) %>%
  group_by(otu) %>% filter(con == max(con))
Not<- data.frame(Change = "Not Sig",
                 con = "0",
                 otu = unique(No$otu)[-which(unique(No$otu) %in% Up$otu)])

phmdata<- rbind(Not, Up) %>%
  left_join(pdata, by= "otu")
colnames(phmdata)<- c("change", "con", "otu", "ck", "om", "pp", "size")

绘图：

ggtern(data=phmdata, aes(x=ck, y=om, z=pp, color= con))+
  geom_point(aes(size=size *3))+
  scale_color_manual(values = c("#B5B5B5", "#7E9D7E", "#811F1F", "#304D75" ))+
  theme(legend.position = "none")+
  labs(x= "ck(20)", y="om(29)", z="pp(26)")

三元图展示：

R语言ggplot2可视化：ggplot2可视化半小提琴图（Half Violin Plots）

R语言ggplot2可视化：ggplot2可视化半小提琴图（Half Violin Plots）

目录

R语言ggplot2可视化：ggplot2可视化半小提琴图（Half Violin Plots）