R调用Taxonkit展示系统发育信息

2024-06-17 00:52
文章标签 展示 调用 信息 系统发育 taxonkit

本文主要是介绍R调用Taxonkit展示系统发育信息，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

Introduction

TaxonKit是一个用于处理生物分类学数据的命令行工具。
它的主要功能是处理NCBI的生物分类学数据，包括对分类单元（如物种、属、科等）的查找、分类单元的上下位关系查询、分类单元名称的标准化等。

为了方便R社区用户（自己）使用和流程整合，我把Taxonkit工具整合进了R包pctax，也开发了一些配套的系统发育分析和可视化方法。

R调用Taxonkit

准备工作

  1. 安装pctax
    pctax稳定版本可在CRAN上获得:
install.packages("pctax")

或者你可以通过以下方式从GitHub安装pctax的开发版本:

# install.packages("devtools")
devtools::install_github("Asa12138/pctax")
  1. 安装taxonkit：
library(pctax)
pctax::install_taxonkit(make_sure = TRUE)#成功后taxonkit会安装在下面这个目录👇
tools::R_user_dir("pctax")
  1. 下载NCBI Taxonomy数据文件:
pctax::download_taxonkit_dataset(make_sure = TRUE)#成功后Taxonomy数据文件会在下面这个目录👇
file.path(Sys.getenv("HOME"), ".taxonkit")

该函数会下载官网最新版本的Taxonomy数据库，如果需要制定版本的数据库，可以自己在官网下载：https://ftp.ncbi.nih.gov/pub/taxonomy/，然后指定位置：

pctax::download_taxonkit_dataset(make_sure = TRUE,taxdump_tar_gz = "~/Downloads/taxdump.tar.gz")

使用

# 下列命令不报错说明可以正常使用
check_taxonkit(print = FALSE)

主要功能与taxonkit一致：

函数功能
taxonkit_list列出指定TaxId下所有子单元的的TaxID
taxonkit_lineage根据TaxID获取完整谱系（lineage）
taxonkit_reformat将完整谱系转化为“界门纲目科属种株”的自定义格式
taxonkit_name2taxid将分类单元名称转化为TaxID
taxonkit_filter按分类学水平范围过滤TaxIDs
taxonkit_lca计算最低公共祖先(LCA)

并且help(taxonkit_*)可查看详细使用说明。

# 列出[genus] Homo下的所有子单元
taxonkit_list(ids = c(9605), indent = "-", show_name = TRUE, show_rank = TRUE)

##  [1] "9605 [genus] Homo"                                    
##  [2] "-9606 [species] Homo sapiens"                         
##  [3] "--63221 [subspecies] Homo sapiens neanderthalensis"   
##  [4] "--741158 [subspecies] Homo sapiens subsp. 'Denisova'" 
##  [5] "-1425170 [species] Homo heidelbergensis"              
##  [6] "-2665952 [no rank] environmental samples"             
##  [7] "--2665953 [species] Homo sapiens environmental sample"
##  [8] "-2813598 [no rank] unclassified Homo"                 
##  [9] "--2813599 [species] Homo sp."                         
## [10] ""

taxonkit_lineage, taxonkit_reformat, taxonkit_name2taxid, taxonkit_filter 与 taxonkit_lca 默认从文件中读取数据，也可通过指定text = TRUE从字符串输入读取输入数据：

# 查询9606和63221的完整谱系
taxonkit_lineage("9606\n63221", show_name = TRUE, show_rank = TRUE, text = TRUE)%>%pcutils::strsplit2(split = "\t",colnames = c("taxid","lineage","name","level"))

##   taxid
## 1  9606
## 2 63221
##                                                                                                                                                                                                                                                                                                                                                                                          lineage
## 1                               cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens
## 2 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens;Homo sapiens neanderthalensis
##                            name      level
## 1                  Homo sapiens    species
## 2 Homo sapiens neanderthalensis subspecies

从文件中读取数据：

names <- system.file("extdata/name.txt", package = "pctax")
taxonkit_name2taxid(names, name_field = 1, sci_name = FALSE, show_rank = FALSE)%>%pcutils::strsplit2(split = "\t",colnames = c("name","taxid"))

##                                              name   taxid
## 1                                    Homo sapiens    9606
## 2            Akkermansia muciniphila ATCC BAA-835  349741
## 3                         Akkermansia muciniphila  239935
## 4                 Mouse Intracisternal A-particle   11932
## 5                                        Wei Shen        
## 6 uncultured murine large bowel bacterium BAC 54B  314101
## 7                       Croceibacter phage P2559Y 1327037

系统发育树

如果是做16S测序的话，在分析过程中就会得到一个带距离的系统发育树。宏基因组分析如果组装MAG后用GTDB-Tk比对数据库后也可以获得有距离的系统发育树。

但有时候我们想要从物种名或taxid获取整齐的谱系信息，用来一个构建系统发育树（层级树，没有真实的距离，只展示包含关系）。这是一个常见的需求，很多文章都会画一个这样的树图来展示自己的数据。

可以实现这个需求的工具有一些：

  • iPhylo：https://iphylo.net/，免费，快速，支持NCBI taxonomy和一些化学物质分类树，赞
  • R包taxtree，很慢
  • PhyloT：https://phylot.biobyte.de/，收费

当然可以使用pctax包快速完成，对于分析流程都在R里做的人来说非常方便:

names <- system.file("extdata/name.txt", package = "pctax")%>%readLines()# 首先通过`name_or_id2df`获取整齐的系统发育分类：
tax_df=name_or_id2df(names,mode = "name")# 去除部分NA，原因可能是学名不标准，或者在新数据库里删除了，因为taxonomy数据库是不断变化的
tax_df=na.omit(tax_df)#用`df2tree`将分类层级表转化为树对象
tax_tree=pctax::df2tree(tax_df[,3:9])# tax_tree是phylo对象，可以用ape包直接简单绘图
ape:::plot.phylo(tax_tree)

可视化

pctax还提供了一些系统发育信息展示方法：

  1. 系统发育树
data(otutab, package = "pcutils")
#otutab是丰度数据，taxonomy是分类层级表（可通过name_or_id2df获得）
ann_tree(taxonomy, otutab) -> treeeasy_tree(tree, add_abundance = TRUE) -> p
p

添加主要Phylum的strip：

easy_tree(tree, add_abundance = TRUE,add_tiplab = FALSE) -> p
some_tax <- table(taxonomy$Phylum) %>%sort(decreasing = TRUE) %>%head(5) %>%names()
add_strip(p, some_tax)

当然，更多系统发育树的绘制可以参考我之前写的R绘制优美的进化树（基础）和R绘制优美的进化树（进阶），或者使用iPhylo网站来交互式绘图：iPhylo 生成并绘制优美的分类树

  1. 桑基图：
sangji_plot(tree)

3.旭日图

sunburst(tree)

TaxonKit 使用

TaxonKit是采用Go语言编写的命令行工具， 提供Linux, Windows, macOS操作系统不同架构（x86-64/arm64）的静态编译的可执行二进制文件。
发布的压缩包不足3Mb，除了Github托管外，还提供国内镜像供下载，同时还支持conda和homebrew安装。

用户只需要下载、解压，开箱即用，无需配置，仅需下载解压NCBI Taxonomy数据文件解压到指定目录即可。

  • 源代码 https://github.com/shenwei356/taxonkit ，
  • 文档 http://bioinf.shenwei.me/taxonkit （介绍、使用说明、例子、教程）

选择系统对应的版本下载最新版 https://github.com/shenwei356/taxonkit/releases ，解压后添加环境变量即可使用。或可选conda安装

conda install taxonkit -c bioconda -y

表格数据处理，推荐使用 csvtk 更高效：

conda install csvtk -c bioconda -y

测试数据下载可直接 https://github.com/shenwei356/taxonkit 下载项目压缩包，或使用git clone下载项目文件夹，其中的example为测试数据

git clone https://github.com/shenwei356/taxonkit

TaxonKit为命令行工具，采用子命令的方式来执行不同功能， 大多数子命令支持标准输入/输出，便于使用命令行管道进行流水作业， 轻松整合进分析流程中。

  • 输出：
    • 所有命令输出中包含输入数据内容，在此基础上增加列。
    • 所有命令默认输出到标准输出（stdout），可通过重定向（>）写入文件。
    • 或通过全局参数-o或--out-file指定输出文件，且可自动识别输出文件后缀（.gz）输出gzip格式。
  • 输入：
    • 除了list与taxid-changelog之外，lineage, reformat, name2taxid, filter 与 lca 均可从标准输入（stdin）读取输入数据，也可通过位置参数（positional arguments）输入，即命令后面不带 任何flag的参数，如 taxonkit lineage taxids.txt
    • 输入格式为单列，或者制表符分隔的格式，输入数据所在列用-i或--taxid-field指定。

TaxonKit直接解析NCBI Taxonomy数据文件（2秒左右），配置更容易，也便于更新数据，占用内存在500Mb-1.5G左右。 数据下载：

# 有时下载失败，可多试几次；或尝试浏览器下载此链接
wget -c https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz 
tar -zxvf taxdump.tar.gz# 解压文件存于家目录中.taxonkit/，程序默认数据库默认目录
mkdir -p $HOME/.taxonkit
cp names.dmp nodes.dmp delnodes.dmp merged.dmp $HOME/.taxonkit

Taxonkit的作者大大贴心地提供了中文文档：https://bioinf.shenwei.me/taxonkit/chinese/，非常详细，大家可以参考使用。

关注公众号，获取最新推送

关注公众号 ‘bio llbug’,获取最新推送。

这篇关于R调用Taxonkit展示系统发育信息的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/1068031。

相关文章

Java调用DeepSeek API的最佳实践及详细代码示例

Java调用DeepSeek API的最佳实践及详细代码示例

《Java调用DeepSeekAPI的最佳实践及详细代码示例》：本文主要介绍如何使用Java调用DeepSeekAPI,包括获取API密钥、添加HTTP客户端依赖、创建HTTP请求、处理响应、... 目录1. 获取API密钥2. 添加HTTP客户端依赖3. 创建HTTP请求4. 处理响应5. 错误处理6.
阅读更多...
Python如何实现PDF隐私信息检测

Python如何实现PDF隐私信息检测

《Python如何实现PDF隐私信息检测》随着越来越多的个人信息以电子形式存储和传输,确保这些信息的安全至关重要,本文将介绍如何使用Python检测PDF文件中的隐私信息,需要的可以参考下... 目录项目背景技术栈代码解析功能说明运行结php果在当今，数据隐私保护变得尤为重要。随着越来越多的个人信息以电子形
阅读更多...
Deepseek R1模型本地化部署+API接口调用详细教程(释放AI生产力)

Deepseek R1模型本地化部署+API接口调用详细教程(释放AI生产力)

《DeepseekR1模型本地化部署+API接口调用详细教程(释放AI生产力)》本文介绍了本地部署DeepSeekR1模型和通过API调用将其集成到VSCode中的过程,作者详细步骤展示了如何下载和... 目录前言一、deepseek R1模型与chatGPT o1系列模型对比二、本地部署步骤1.安装oll
阅读更多...
一分钟带你上手Python调用DeepSeek的API

一分钟带你上手Python调用DeepSeek的API

《一分钟带你上手Python调用DeepSeek的API》最近DeepSeek非常火,作为一枚对前言技术非常关注的程序员来说,自然都想对接DeepSeek的API来体验一把,下面小编就来为大家介绍一下... 目录前言免费体验API-Key申请首次调用API基本概念最小单元推理模型智能体自定义界面总结前言最
阅读更多...
JAVA调用Deepseek的api完成基本对话简单代码示例

JAVA调用Deepseek的api完成基本对话简单代码示例

《JAVA调用Deepseek的api完成基本对话简单代码示例》：本文主要介绍JAVA调用Deepseek的api完成基本对话的相关资料,文中详细讲解了如何获取DeepSeekAPI密钥、添加H... 获取API密钥首先，从DeepSeek平台获取API密钥，用于身份验证。添加HTTP客户端依赖使用Jav
阅读更多...
redis防止短信恶意调用的实现

redis防止短信恶意调用的实现

《redis防止短信恶意调用的实现》本文主要介绍了在场景登录或注册接口中使用短信验证码时遇到的恶意调用问题,并通过使用Redis分布式锁来解决,具有一定的参考价值,感兴趣的可以了解一下... 目录1.场景2.排查3.解决方案3.1 Redis锁实现3.2 方法调用1.场景登录或注册接口中，使用短信验证码场
阅读更多...
使用C/C++调用libcurl调试消息的方式

使用C/C++调用libcurl调试消息的方式

《使用C/C++调用libcurl调试消息的方式》在使用C/C++调用libcurl进行HTTP请求时,有时我们需要查看请求的/应答消息的内容（包括请求头和请求体）以方便调试,libcurl提供了多种... 目录1. libcurl 调试工具简介2. 输出请求消息使用 CURLOPT_VERBOSE使用 C
阅读更多...
C#实现系统信息监控与获取功能

C#实现系统信息监控与获取功能

《C#实现系统信息监控与获取功能》在C#开发的众多应用场景中,获取系统信息以及监控用户操作有着广泛的用途,比如在系统性能优化工具中,需要实时读取CPU、GPU资源信息,本文将详细介绍如何使用C#来实现... 目录前言一、C# 监控键盘1. 原理与实现思路2. 代码实现二、读取 CPU、GPU 资源信息1.
阅读更多...
Python调用另一个py文件并传递参数常见的方法及其应用场景

Python调用另一个py文件并传递参数常见的方法及其应用场景

《Python调用另一个py文件并传递参数常见的方法及其应用场景》：本文主要介绍在Python中调用另一个py文件并传递参数的几种常见方法,包括使用import语句、exec函数、subproce... 目录前言1. 使用import语句1.1 基本用法1.2 导入特定函数1.3 处理文件路径2. 使用ex
阅读更多...
在C#中获取端口号与系统信息的高效实践

在C#中获取端口号与系统信息的高效实践

《在C#中获取端口号与系统信息的高效实践》在现代软件开发中,尤其是系统管理、运维、监控和性能优化等场景中,了解计算机硬件和网络的状态至关重要,C#作为一种广泛应用的编程语言,提供了丰富的API来帮助开... 目录引言1. 获取端口号信息1.1 获取活动的 TCP 和 UDP 连接说明：应用场景：2. 获取硬
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2