首页 > 其他分享 >R：三元图

R：三元图

时间：2023-06-12 10:58:38浏览次数：35  
标签：enrich OTUs value 三元 data otu size

setwd("D:\\Desktop")
library(tidyverse)
# 数据处理函数
data_clean <- function(otu, design, type=c("relative", "absolute"), threshold=0.001, times=100){

# 函数测试数据
# library(amplicon)
# otu=otutab
# metadata$SampleID=rownames(metadata)
# design=metadata[,c("SampleID","Group")]
# type="absolute"
# threshold=0.0005
# times=100

# 绝对丰度转相对丰度1
if (type == "absolute"){
otu_relative <- apply(otu, 2, function(x){x/sum(x)})
}else {otu_relative <- otu}

# 至少有一个样本大于阈值即保留
idx <- rowSums(otu_relative > threshold) >= 1
otu_threshold <- as.data.frame(otu_relative[idx, ])
otu_threshold$OTUs <- row.names(otu_threshold)

#转换宽表格为长表格
otu_longer <- pivot_longer(data=otu_threshold,
cols=-OTUs,
names_to="SampleID",
values_to="value")

# 按"SampleID"对应添加元数据中的分组Group
merge_data <- merge(otu_longer, design, by ="SampleID")
# 去除样本列
# otu <- subset(merge_data, select=-SampleID)
# 元数据不只有样本列，直接筛选OTUs、Group和value更稳健
otu <- subset(merge_data, select=c("Group","OTUs","value"))
# 按OTUs和Group求均值
otu_mean <- otu %>% group_by(OTUs, Group) %>%
summarise(value=mean(value))
# 转换回宽表格
otu_tern <- otu_mean %>%
group_by(Group, OTUs) %>%
mutate(index=row_number()) %>%
pivot_wider(names_from=Group,values_from=value) %>%
select(-index)
# 此处用group_by可以直接合并均值，不用长、宽转换两次

# 调整点大小均值，可缩放 size of points
otu_tern$size <- (apply(otu_tern[2:4], 1, mean))*times
return(otu_tern)
}

# 读取输入文件
otutab <- read.delim("otutab.txt", header=T, row.names=1)
design <- read.delim("metadata.txt", header=T, row.names=NULL)
# 只提取元数据中的样本名和分组列，要求名称为SampleID和Group
design = design[,c("SampleID","Group")]

# 计算三元图输入数据：各组相对丰度均值
otu_tern <- data_clean(otutab, design, type="absolute", threshold=0.001, times=100)
head(otu_tern,n=3)

library(ggtern)
p <- ggtern(data=otu_tern, aes(x=KO, y=OE, z=WT)) +
geom_point(aes(size=size), alpha=0.8, show.legend=T) +
scale_size(range=c(0, 6)) + geom_mask() +
guides(colour="none") + theme_bw() +
theme(axis.text=element_blank(), axis.ticks=element_blank())
p
ggsave(paste0("p1.png"), p, width=89, height=59, units="mm")
ggsave(paste0("p1.pdf"), p, width=89, height=59, units="mm")


###高丰度特征着色
top_OTUs <- function(data, rank=10){
# 按丰度size降序排列
data_order <- data[order(data$size, decreasing=T), ]
# if (missing(rank))
# rank=10
# 提取前N行
top <- data_order[1:rank, ]
# 其余部分
otu_el <- data_order[-(1:rank), ]
# 返回top和其它的列表
return(list(top, otu_el))
}
plot_data <- top_OTUs(otu_tern, rank=10)

# 配色方案
platte <- c('#a6cee3', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99',
'#e31a1c', '#fdbf6f', '#ff7f00', '#cab2d6', '#6a3d9a')

p <- ggtern(data=otu_tern,
aes(x=KO, y=OE, z=WT)) +
geom_mask() +
geom_point(data=plot_data[[2]], aes(size=size), color="grey",
alpha=0.8, show.legend=F) +
geom_point(data=plot_data[[1]], aes(size=size, color=OTUs),
show.legend=T) +
scale_colour_manual(values=platte) +
scale_size(range=c(0, 6)) +
# legend
guides(size="none") +
theme_bw() +
theme(axis.text=element_blank(),
axis.ticks=element_blank())
p
ggsave(paste0("p2.png"), p, width=89, height=59, units="mm")
ggsave(paste0("p2.pdf"), p, width=89, height=59, units="mm")


###富集显著性分析
library(edgeR)
# 计算3组比较各组特有的特征
enrich_data <- function(otu, design, p.value=0.05, adjust.method="fdr"){

# 函数测试数据
# library(amplicon)
# otu=otutab
# metadata$SampleID=rownames(metadata)
# design=metadata[,c("SampleID","Group")]
# p.value=0.05
# adjust.method="fdr"

dge_list <- DGEList(counts=otu, group=design$Group)
# Remove the lower abundance/(cpm, rpkm)
keep <- rowSums(dge_list$counts) >= 0
dge_keep <- dge_list[keep, ,keep.lib.sizes=F]
# scale the raw library sizes dgelist
dge <- calcNormFactors(dge_keep)
# fit the GLM
design.mat <- model.matrix(~ 0 + dge$samples$group)
d2 <- estimateGLMCommonDisp(dge, design.mat)
d2 <- estimateGLMTagwiseDisp(d2, design.mat)
fit <- glmFit(d2, design.mat)
#######
# if (missing(adjust.method))
# adjust.method="fdr"
# if (missing(p.value))
# p.value=0.05
group_index <- as.character(design$Group[!duplicated(design$Group)])
# enrich groups
lrt_1_2 <- glmLRT(fit, contrast=c(1, -1, 0))
lrt_1_3 <- glmLRT(fit, contrast=c(1, 0, -1))

de_1_2 <- decideTestsDGE(lrt_1_2, adjust.method=adjust.method,
p.value=p.value)
de_1_3 <- decideTestsDGE(lrt_1_3, adjust.method=adjust.method,
p.value=p.value)

rich_1 <- rownames(otu)[de_1_2 == 1 & de_1_3 == 1]
enrich_1 <- data.frame(OTUs=rich_1,
enrich=rep(group_index[1], length(rich_1)))
###############################
lrt_2_3 <- glmLRT(fit, contrast=c(0, 1, -1))
lrt_2_1 <- glmLRT(fit, contrast=c(-1, 1, 0))

de_2_3 <- decideTestsDGE(lrt_2_3, adjust.method=adjust.method,
p.value=p.value)
de_2_1 <- decideTestsDGE(lrt_2_1, adjust.method=adjust.method,
p.value=p.value)

rich_2 <- rownames(otu)[de_2_3 == 1 & de_2_1 == 1]
enrich_2 <- data.frame(OTUs=rich_2,
enrich=rep(group_index[2], length(rich_2)))
###################
lrt_3_1 <- glmLRT(fit, contrast=c(-1, 0, 1))
lrt_3_2 <- glmLRT(fit, contrast=c(0, -1, 1))

de_3_1 <- decideTestsDGE(lrt_3_1, adjust.method=adjust.method,
p.value=p.value)
de_3_2 <- decideTestsDGE(lrt_3_2, adjust.method=adjust.method,
p.value=p.value)

rich_3 <- rownames(otu)[de_3_1 == 1 & de_3_2 == 1]
enrich_3 <- data.frame(OTUs=rich_3,
enrich=rep(group_index[3], length(rich_3)))
enrich_index <- rbind(enrich_1, enrich_2, enrich_3)
return(enrich_index)
}
enrich_index <- enrich_data(otutab, design, p.value=0.05)
plot_data <- merge(otu_tern, enrich_index, by="OTUs", all.x=T)
p <- ggtern(data=plot_data,
aes(x=KO, y=OE, z=WT)) +
geom_mask() + # 可将超出边界的点正常显示出来
geom_point(aes(size=size, color=enrich),alpha=0.8) +
guides(size="none") +theme_bw() +
theme(axis.text=element_blank(),
axis.ticks=element_blank())
p
ggsave(paste0("p3.png"), p, width=89, height=59, units="mm")
ggsave(paste0("p3.pdf"), p, width=89, height=59, units="mm")

#######################################################################################################################################

 

这段代码是一个R脚本，用于进行数据处理和绘制三元图。解释：

    1. setwd("D:\\Desktop")：设置当前工作目录为桌面（Windows系统）。
    2. library(tidyverse)：加载tidyverse包，包含了一组用于数据处理和可视化的R包。
    3. data_clean <- function(otu, design, type=c("relative", "absolute"), threshold=0.001, times=100)：定义了一个名为data_clean的函数，用于数据处理。
    4. 函数data_clean中的代码是数据处理的具体实现，包括相对丰度转换、阈值筛选、长宽表格转换等操作。
    5. otutab <- read.delim("otutab.txt", header=T, row.names=1)：从名为"otutab.txt"的文件中读取OTU表格数据，设置第一行为列名。
    6. design <- read.delim("metadata.txt", header=T, row.names=NULL)：从名为"metadata.txt"的文件中读取元数据，设置第一行为列名。
    7. design = design[,c("SampleID","Group")]：保留元数据中的"SampleID"和"Group"两列。
    8. otu_tern <- data_clean(otutab, design, type="absolute", threshold=0.001, times=100)：调用data_clean函数，对输入数据进行处理，生成三元图的输入数据。
    9. head(otu_tern,n=3)：显示otu_tern数据的前3行。
    10. library(ggtern)：加载ggtern包，用于绘制三元图。
    11. p <- ggtern(data=otu_tern, aes(x=KO, y=OE, z=WT)) + ...：创建一个三元图对象p，设置数据来源和图形属性，包括点的大小、颜色等。
    12. ggsave(paste0("p1.png"), p, width=89, height=59, units="mm")：将图形p保存为名为"p1.png"的PNG图像文件。
    13. ggsave(paste0("p1.pdf"), p, width=89, height=59, units="mm")：将图形p保存为名为"p1.pdf"的PDF文件。
    14. top_OTUs <- function(data, rank=10)：定义了一个名为top_OTUs的函数，用于提取高丰度特征。
    15. 函数top_OTUs中的代码用于对数据进行排序和提取高丰度特征。
    16. plot_data <- top_OTUs(otu_tern, rank=10)：调用top_OTUs函数，提取otu_tern数据中的前10个高丰度特征。
    17. 定义了一个颜色方案platte，包含多个颜色值。
    18. p <- ggtern(data=otu_tern, aes(x=KO, y=OE, z=WT)) + ...：重新创建一个三元图对象p，设置数据来源和图形属性，包括点的大小、颜色等。
    19. ggsave(paste0("p2.png"), p, width=89, height=59, units="mm")：将图形p保存为名为"p2.png"的PNG图像文件。
    20. ggsave(paste0("p2.pdf"), p, width=89, height=59, units="mm")：将图形p保存为名为"p2.pdf"的PDF文件。
    21. enrich_data <- function(otu, design, p.value=0.05, adjust.method="fdr")：定义了一个名为enrich_data的函数，用于富集显著性分析。
    22. 函数enrich_data中的代码是进行富集显著性分析的具体实现，使用edgeR包进行差异分析。
    23. enrich_index <- enrich_data(otutab, design, p.value=0.05)：调用enrich_data函数，对otutab数据进行富集显著性分析，生成富集特征的索引。
    24. plot_data <- merge(otu_tern, enrich_index, by="OTUs", all.x=T)：将otu_tern数据和enrich_index数据按OTUs列进行合并，生成用于绘图的数据。
    25. p <- ggtern(data=plot_data, aes(x=KO, y=OE, z=WT)) + ...：重新创建一个三元图对象p，设置数据来源和图形属性，包括点的大小、颜色等。
    26. ggsave(paste0("p3.png"), p, width=89, height=59, units="mm")：将图形p保存为名为"p3.png"的PNG图像文件。
    27. ggsave(paste0("p3.pdf"), p, width=89, height=59, units="mm")：将图形p保存为名为"p3.pdf"的PDF文件。

    28. 这段代码主要用于进行数据处理、绘制三元图和进行富集显著性分析。下面是对代码的功能和生成的图片的解释：

      1. 数据处理：

        • 从"otutab.txt"文件读取OTU表格数据和"metadata.txt"文件读取元数据。
        • 进行数据清洗、转换和筛选，生成用于三元图绘制和富集分析的数据。

      2. 生成的图片：

        • p1.png和p1.pdf：展示了原始数据的三元图，其中每个点代表一个样本，KO、OE、WT分别表示不同条件下的特征值。
        • p2.png和p2.pdf：展示了筛选后的高丰度特征的三元图，用不同颜色表示不同的特征。
        • p3.png和p3.pdf：在p2的基础上，通过富集显著性分析标注了富集的特征，用不同颜色表示富集程度。

      综上所述，该代码利用输入的OTU表格数据和元数据进行数据处理，生成了原始数据三元图、高丰度特征三元图以及标注了富集特征的三元图。这些图像可以帮助研究人员进行数据可视化和富集分析，以便更好地理解数据并发现其中的模式和关联。

    29. 在上述代码中，KO、OE和WT分别代表不同的条件或组别。这些条件或组别可以根据实际研究设计来确定，具体含义可能取决于数据的特定上下文。通常，它们可以表示以下内容：

      • KO：代表基因敲除（Knockout）组，通常是指在实验中进行基因敲除或功能抑制的样本组。
      • OE：代表基因过表达（Overexpression）组，通常是指在实验中进行基因过度表达的样本组。
      • WT：代表野生型（Wild-Type）组，通常是指未进行任何基因操作或干预的对照组，代表自然状态或正常表达水平。

      这些条件的选择和命名可以根据具体实验的目的和研究对象进行调整。在三元图中，每个点的位置表示不同条件下的特征值，并且通过颜色或大小来表示特征的丰度或其他相关信息。

标签：enrich,OTUs,value,三元,data,otu,size
From： https://www.cnblogs.com/wzbzk/p/17474321.html

相关文章

  • 三元图
    三元图（TernaryPlot）广泛用于三个分组数据比较、筛选，通过三元图可以直观展示数据在三个分组的分布情况，高效率地筛选离群元素，同时配合方差分析等统计检验方法可以找到不同分组中显著富集的元素。三元图由三个坐标轴组成等边三角形，轴上数值代表对应分组占比数值，三个顶点标注的信息代......
  • 二分法 三元表达式 生成式 匿名函数 内置函数
    目录二分法三元表达式生成式列表生成式字典生成式集合生成式元组生成式(生成器）匿名函数内置函数二分法二分法思路1.二分法的使用前提条件：列表中得数字必须要有序2.将对象整除2分成两部分3.将目标数值与分割后的对象做比较来确定目标数值在哪一部分4.继续重复这两个步骤直至......
  • 算法之二分法、三元表达式、列表生成式、字典生成式(了解)、匿名函数、常见的内置函数
    算法之二分法二分概念二分算法，又称折半查找，即在一个单调有序的集合中查找一个解。每次分为左右两部分，判断解在哪个部分中并调整上下界，直到找到目标元素，每次二分后都将舍弃一半的查找空间。定义and实现：算法就是解决问题的高效办法常见的算法：二分法算法还可以锻炼我们的......
  • day16 Python下的三元运算符
    Python下的三元运算符【一】引言三元表达式（三目运算符）能够简洁我们的代码三元表达式其实是将if...else...判断语句的简化表达，代替很多ifelse和if-else一样，只有一个表达式会被执行。因此，三元表达式中的if和else可以包含大量的计算，但只有True的分支会被执行在Java、C......
  • HDU3939(毕达哥拉斯三元组的解)
    对于方程：，满足条件：x，y，z两两互素的正整数解为： ，其中m>n>0，gcd(m,n)=1，m，n一奇一偶。典型题目：POJ1305，HDU3939 对于POJ1305很简单，下面重点来解析HDU3939题。题目：SticksandRightTriangle#include<iostream>#include<string.h>#include<algorithm>#include<stdio.h>#include&l......
  • 推导式，三元运算符，匿名函数（lambda）
    set集合只有key的字典{}set中的元素必须是可hash的，也是不可变的元素是无序，不重复的set()转为集合，可以用来去重增.add()添加.update()迭代更新删.remove()删，返回值是none.clear()清空，空集合是set()改先删再加查for循环集合的交......
  • （七) 关系运算符，逻辑运算符，赋值运算符，三元运算符
    目录关系运算符一览逻辑运算符一览逻辑运算规则逻辑与和短路与区别逻辑或和短路或的区别关系运算符一览关系运算符的结果都是boolean类型,要么为true,要么为false关系运算符组成的表达式称为关系表达式,如a>=b逻辑运算符一览用于连接多个条件(多个关系表达式),最终......
  • Day 28 28.1 JS进阶之三元运算符
    JS工具之三元运算符【1】格式三元运算符:条件表达式?语句1:语句2;leta=10;letb=20;letd=a>b?a:bconsole.log(d);注释：条件运算符在执行时，首先对条件表达式进行求值，如果该值为true，则执行语句1，并返回执行结果如果该值为false，则......
  • 2017ACM/ICPC广西邀请赛-重现赛（感谢广西大学）C - Counting Stars 暴力三元环
    LittleAisanastronomylover,andhehasfoundthattheskywassobeautiful!Soheiscountingstarsnow!Therearenstarsinthesky,andlittleAhasconnectedthembymnon-directionaledges.Itisguranteedthatnoedgesconnectonestarwithi......
  • go模拟三元表达式最简单的方式
    众所周知，Go语言本身并没有提供内置的三元表达式，但是我们可以使用一些技巧来模拟实现。下面是使用最短的代码实现Go的三元表达式：packagemainimport"fmt"funcmain(){x:=10y:=20max:=map[bool]int{true:x,false:y}[x>y]fmt.Println(max)}......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99