PICRUSt2下载

PICRUSt2是一个用于预测微生物基因组功能的工具，它可以通过基因家族注释和基于参考基因组的比对来对16S rRNA或宏基因组数据进行功能预测。在R中使用PICRUSt2需要先安装相关的R包，然后按照以下步骤操作： 1. 安装依赖包和PICRUSt2 ```r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c("phyloseq", "dplyr", "stringr", "tidyr", "magrittr")) ``` ```r if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools") devtools::install_github("picrust/picrust2", build_vignettes = TRUE) ``` 2. 加载依赖包和PICRUSt2 ```r library(phyloseq) library(dplyr) library(stringr) library(tidyr) library(magrittr) library(picrust2) ``` 3. 读取OTU表和物种注释信息 ```r otu_file <- "otu_table.biom" otu_table <- import_biom(otu_file) tax_file <- "taxonomy.tsv" tax_table <- read.delim(tax_file, header = FALSE, stringsAsFactors = FALSE) colnames(tax_table) <- c("OTU_ID", "taxonomy") ``` 4. 处理物种注释信息 ```r tax_table %<>% separate(taxonomy, into = c("kingdom", "phylum", "class", "order", "family", "genus", "species"), sep = ";") tax_table$OTU_ID <- str_remove(tax_table$OTU_ID, ".*OTU_") ``` 5. 将OTU表和物种注释信息合并成phyloseq对象 ```r tax_table_phylo <- otu_table_to_phyloseq( otu_table = otu_table, tax_table = tax_table, refseq_col = "OTU_ID" ) ``` 6. 运行PICRUSt2 ```r picrust2_pipeline(tax_table_phylo, outdir = "picrust2_output") ``` 在运行过程中，PICRUSt2会下载和安装必要的参考数据库，并生成预测的功能注释文件。最后，可以使用以下命令读取并处理预测的功能注释文件： ```r ko_file <- "picrust2_output/metagenome_predictions/ko_metagenome_contributions.tsv.gz" ko_table <- read.delim(ko_file, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) colnames(ko_table)[1] <- "OTU_ID" ko_table$OTU_ID <- str_remove(ko_table$OTU_ID, ".*OTU_") ko_table %<>% pivot_longer(cols = -OTU_ID, names_to = "KO_ID", values_to = "abundance") %>% group_by(OTU_ID, KO_ID) %>% summarize(abundance = sum(abundance)) %>% ungroup() # 根据KO ID获取KEGG通路信息 ko_to_pathway_file <- "picrust2_output/pathways_out/predicted_metagenome_unstratified.tsv.gz" ko_to_pathway_table <- read.delim(ko_to_pathway_file, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE) colnames(ko_to_pathway_table) <- c("KO_ID", "pathway_ID", "abundance") ``` 处理完毕后，就可以使用R中的其它工具对预测的功能注释数据进行进一步分析和可视化。

<think>好的，用户之前已经询问过如何根据16S rRNA结果绘制KEGG代谢图，现在他们进一步提问，有了PICRUSt2和Tax4Fun2的结果后，具体如何在KEGG上绘制代谢图。我需要详细解释步骤，确保他们能理解如何从预测结果到可视化。 首先，用户可能已经完成了功能预测，现在需要知道如何处理这些结果，比如KO表的获取和整理。可能需要说明不同工具的输出格式，例如PICRUSt2生成的pred_metagenome_unstrat.tsv，而Tax4Fun2的KEGG输出可能需要转换。 接下来，用户需要将KO映射到KEGG通路。这里需要解释如何获取KEGG的映射关系，可能提到使用KEGG的API或者R包如KEGGREST，或者工具内置的功能。同时，要指导如何将KO丰度聚合到通路上，可能涉及求和或平均值，并注意数据库版本的问题。 然后是可视化工具的选择，用户可能需要不同的展示方式，比如整体代谢网络用iPath3，具体通路用Pathview，或者自定义图表用ggplot2。需要给出具体的例子，比如R代码使用Pathview生成通路图，或者Python生成iPath3的输入文件。 还需要注意注意事项，比如功能预测的局限性，数据库更新，以及统计方法的应用，如归一化和差异检验。用户可能需要这些信息来确保结果的准确性。 另外，用户可能对具体操作步骤不熟悉，需要分步说明，并给出资源链接供进一步学习。还要考虑用户可能遇到的实际问题，比如如何安装R包，如何处理文件路径等，但可能不需要太深入，除非用户后续提问。 最后，确保回答结构清晰，步骤明确，工具对比表格帮助用户选择合适的方法，并提供示例代码和注意事项，帮助用户顺利从预测结果过渡到可视化。</think>以下是基于PICRUSt2和Tax4Fun2的功能预测结果绘制KEGG代谢通路的完整流程： --- ### 一、数据预处理 #### 1. **功能预测结果提取** - **PICRUSt2输出文件** ```bash # 核心文件 pred_metagenome_unstrat.tsv # 未分层的KO丰度表 pathways_out/path_abun_unstrat.tsv # 预计算的通路丰度 ``` - **Tax4Fun2输出文件** ```R # R输出对象中包含： $KEGG_Abundance # KEGG通路丰度矩阵 $KO_Abundance # KEGG Ortholog丰度矩阵 ``` #### 2. **KO丰度表整理** - **格式要求** 矩阵需满足：**行=样本，列=KO ID**，例如： | Sample | K00001 | K00002 | ... | |--------|--------|--------|-----| | S1 | 15.2 | 8.7 | ... | --- ### 二、KEGG通路映射 #### 1. **获取KEGG元数据** - **手动下载映射表** 从KEGG官网下载层级关系： ```bash wget https://www.kegg.jp/kegg-bin/download_htext?htext=ko00001.keg -O kegg_hierarchy.txt ``` - **R语言自动化获取** ```R library(KEGGREST) # 获取KO到通路的映射 ko2pathway <- keggLink("pathway", "ko") # 保存为数据框 write.table(ko2pathway, "ko2pathway.tsv", sep="\t") ``` #### 2. **丰度聚合计算** - **计算公式** 对每个通路的KO丰度求和： $$ \text{Pathway}_i = \sum_{j \in \text{KO}_i} \text{KO}_j $$ - **Python实现示例** ```python import pandas as pd ko_table = pd.read_csv("ko_abundance.csv", index_col=0) pathway_map = pd.read_csv("ko2pathway.tsv", sep="\t", header=None) # 聚合到通路层级 pathway_table = ko_table.T.groupby(pathway_map[1]).sum().T ``` --- ### 三、可视化方法详解 #### 方法1：Pathview绘制通路拓扑图 **适用场景**：展示特定通路内各KO的丰度分布 ```R library(pathview) # 输入数据要求：命名向量（样本需单独绘制） ko_vector <- as.numeric(ko_table[1, ]) # 取第一个样本 names(ko_vector) <- colnames(ko_table) # 绘制糖酵解通路（ko00010） pv <- pathview( gene.data = ko_vector, pathway.id = "00010", species = "ko", # 使用KEGG通用通路 kegg.dir = "./kegg/", # 指定本地数据库路径 gene.idtype = "KEGG", plot.col.key = FALSE ) ``` #### 方法2：iPath3交互式可视化 **适用场景**：全局代谢网络探索 1. **准备输入文件** 生成包含KO丰度的TSV文件： ```text #KO Value K00001 15.2 K00002 8.7 ... ``` 2. **上传到iPath3** 访问 https://pathways.embl.de/ipath3 ，上传文件后调整参数： - **Color Mapping**：选择"Custom Data"映射丰度颜色 - **Layer Selection**：勾选"Enzymes"和"Metabolites" #### 方法3：ggplot2高级绘图 **示例：绘制通路丰度热图** ```R library(pheatmap) # 数据标准化（Z-score） pathway_z <- t(scale(t(pathway_table))) pheatmap(pathway_z, clustering_method = "ward.D2", color = colorRampPalette(c("navy", "white", "firebrick3"))(50), show_rownames = TRUE, main = "KEGG Pathway Abundance Heatmap") ``` --- ### 四、关键注意事项 1. **数据归一化** - 推荐使用**CSS标准化**（MetagenomeSeq包）消除测序深度差异 ```R library(metagenomeSeq) mrobj <- newMRexperiment(t(ko_table)) mrobj_norm <- cumNorm(mrobj, p=0.5) ko_norm <- MRcounts(mrobj_norm, norm=TRUE) ``` 2. **通路筛选策略** - **丰度过滤**：保留在至少30%样本中存在的通路 - **方差过滤**：保留top 20%高变异通路 - **生物学相关性**：聚焦碳代谢（ko01200）、氮代谢（ko00910）等核心通路 3. **差异分析整合** 结合LEfSe分析显著差异通路： ```bash lefse-format_input.py pathway_abundance.tsv lefse_input.in -c 1 -u 2 -o 1000000 lefse-run.py lefse_input.in lefse_results.res ``` --- ### 五、工具对比表 | 工具 | 输入数据 | 可视化类型 | 交互性 | 安装难度 | |------------|-------------------|------------------|--------|----------| | Pathview | KO丰度 + 通路ID | 静态通路拓扑图 | 无 | ★★☆☆☆ | | iPath3 | KO列表 + 丰度值 | 全局代谢网络 | 高 | 无需安装 | | Cytoscape | 通路-节点关系表 | 自定义网络图 | 中 | ★★★☆☆ | | ggplot2 | 整理后的数据框 | 任意统计图形 | 无 | ★★★★☆ | --- ### 六、常见问题解决 **Q1: 出现"pathway ID not found"错误** - 检查KEGG数据库版本是否匹配（PICRUSt2默认使用旧版） - 更新本地数据库： ```R download.kegg(pathway.id="ko", species="ko") ``` **Q2: 通路图中节点颜色不变化** - 确认输入的丰度值范围：Pathview默认仅对log2(fold-change)有明显颜色映射 - 调整颜色阈值： ```R pathview(..., limit=list(gene=c(-2,2))) ``` **Q3: iPath3显示"KO not mapped"** - 检查KO编号格式是否正确（需带"K"前缀，如K00001） - 确认使用的iPath3版本是否支持最新KEGG注释 如果需要处理具体数据案例，可提供部分KO丰度表，可针对性给出参数优化建议。