运行你代码出错了> # 硅藻群落-环境因子关联分析（优化版）
> # 修改日期：2025年10月20日
> # 主要优化：错误处理、日志记录、代码健壮性
>
> # ============================== 1. 加载必需包 ==============================
> library(ggplot2)
> library(readxl)
> library(dplyr)
> library(vegan)
> library(ggrepel)
> library(plotrix)
> library(ape)
> library(psych)
> library(corrplot)
> library(openxlsx)
> library(phyloseq)
> library(ggcorrplot)
> theme_set(theme_bw()) # 设置默认主题
>
> # 安装缺失包的函数（按需启用）
> # if (!require("pacman")) install.packages("pacman")
> # pacman::p_load(ggplot2, readxl, dplyr, vegan, ggrepel, plotrix, ape, psych, corrplot, openxlsx, phyloseq, ggcorrplot)
>
> # ============================== 2. 读取硅藻数据 ==============================
> # 工作目录
> setwd("D:/p1")
>
> # 创建输出目录（递归创建，避免报错）
> output_dir <- "result/alpha_diatom/env_corr"
> if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
>
> # ① 读取硅藻抽平OTU表
> otu <- read.delim(
+ file = "result/raw/otutab_diatom_rare.txt",
+ header = TRUE,
+ sep = "\t",
+ row.names = 1,
+ check.names = FALSE
+ )
>
> # ② 读取硅藻分类表
> tax <- read.delim(
+ file = "result/raw/diatom_otus_tax.txt",
+ header = FALSE,
+ sep = "\t",
+ col.names = c("OTUID", "Taxonomy"),
+ check.names = FALSE
+ )
>
> # ③ 读取元数据
> factors <- read.delim(
+ file = "result/metadata.txt",
+ header = TRUE,
+ sep = "\t",
+ row.names = "ID",
+ check.names = FALSE
+ )
>
> # ============================== 3. 分类表处理（优化版） ==============================
> # 使用tax_table构造函数直接处理，避免手动矩阵操作
> # 步骤：拆分Taxonomy列，生成分类矩阵
> tax_split <- strsplit(tax$Taxonomy, ";")
> # 确定最大层级数（通常为7，但以防不一致）
> max_levels <- max(sapply(tax_split, length))
> tax_levels <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")[1:max_levels]
>
> # 填充分类矩阵
> tax_mat <- matrix(NA, nrow = nrow(tax), ncol = max_levels, dimnames = list(tax$OTUID, tax_levels))
> for (i in 1:nrow(tax)) {
+ current_tax <- tax_split[[i]]
+ if (length(current_tax) > 0) {
+ # 去除前缀（如"p:"）
+ cleaned_tax <- gsub("^[a-z]:", "", current_tax)
+ tax_mat[i, 1:length(cleaned_tax)] <- cleaned_tax
+ }
+ }
>
> # 添加OTUID列和原始Taxonomy列，构建数据框
> tax_df <- data.frame(
+ OTUID = tax$OTUID,
+ Taxonomy = tax$Taxonomy,
+ tax_mat,
+ stringsAsFactors = FALSE
+ )
> rownames(tax_df) <- tax_df$OTUID
>
> # 生成NMstr列（非NA分类层级的合并）
> tax_df$NMstr <- apply(tax_mat, 1, function(x) paste(na.omit(x), collapse = "|"))
>
> # ============================== 4. 构建phyloseq对象（强化匹配） ==============================
> # 4.1 转换OTU表
> OTU <- otu_table(otu, taxa_are_rows = TRUE)
>
> # 4.2 转换分类表（只保留分类层级列和NMstr）
> tax_columns <- c(tax_levels, "NMstr")
> tax_sub <- tax_df[, tax_columns, drop = FALSE]
错误于`[.data.frame`(tax_df, , tax_columns, drop = FALSE):
选择了未定义的列
最新发布
这篇博客介绍了如何在R语言中利用tidyr包的drop_na函数来过滤数据框中包含空值（NA）的行。示例代码展示了针对特定列col_a和col_b，以及所有列的过滤方法，帮助读者掌握数据预处理的关键步骤。
扫一扫
1582
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
