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

news2024/10/6 16:22:16

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_filtertaxonkit_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) -> tree

easy_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格式。
  • 输入:
    • 除了listtaxid-changelog之外,lineage, reformat, name2taxid, filterlca 均可从标准输入(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’,获取最新推送。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1829356.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

QuickLook最强大的C#开源文件预览神器

功能特点&#xff1a; 快速预览&#xff1a;用户只需选中文件并按下空格键&#xff0c;即可立即查看文件内容&#xff0c;无需打开特定应用程序或软件。 多格式支持&#xff1a;QuickLook支持预览几乎所有常见的文件类型&#xff0c;包括但不限于&#xff1a; 图像&#xff1…

【计算机网络】已解决:“‘ping‘ 不是内部或外部命令,也不是可运行的程序或批处理文件”报错

文章目录 一、问题分析背景二、可能出错的原因三、错误代码示例四、正确解决方法与示例五、注意事项 已解决“‘ping’ 不是内部或外部命令&#xff0c;也不是可运行的程序或批处理文件”报错 一、问题分析背景 在Windows操作系统中&#xff0c;ping 命令是一个常用的网络诊断…

STM32CubeMX配置-外部中断配置

一、简介 MCU为STM32G070&#xff0c;配置为上升沿触发外部中断&#xff0c;在上升沿外部中断回调函数中进行相关操作。 二、外部中断配置 查看规格书中管教描述&#xff0c;找到I/O对应的外部中断线&#xff0c;然后进行如下上升沿触发外部中断配置。 三、生成代码 调用上升沿…

【Linux】进程_6

文章目录 五、进程8. 进程地址空间 未完待续 五、进程 8. 进程地址空间 上图可能很多人都看过了&#xff0c;这里再来验证一下&#xff1a; 验证位置&#xff1a; 验证堆栈的生长方向&#xff1a; 在上面的空间布局图中&#xff0c;有一个疑问&#xff0c;画的空间是 内存…

【探索Linux】P.34(HTTPS协议)

阅读导航 引言一、HTTPS是什么1. 什么是"加密"2. 为什么要加密3. 常见的加密方式&#xff08;1&#xff09;对称加密&#xff08;2&#xff09;非对称加密 二、证书认证1. CA认证 三、HTTPS的加密底层原理✅非对称加密对称加密证书认证 温馨提示 引言 在上一篇文章中…

Ubuntu 22.04 解决 firefox 中文界面乱码

问题复现 在为Ubuntu 22.04 Server安装完整的GNOME 42.01桌面后&#xff0c;将桌面语言设置为中文时&#xff0c;打开Firefox可能会出现中文乱码的问题。经过网上调查发现&#xff0c;这个问题是由Snap软件包引起的。 解决方案 为了避免在Ubuntu 22.04中文模式下的乱码问题…

消息队列-RabbitMQ-消息确认机制

为了保证消息的不丢失&#xff0c;可靠抵达&#xff0c;可以使用事务消息&#xff0c;但是性能会下降250倍&#xff0c;为此引入确认机制。 1.ConfirmCallBack 服务器收到消息就回调 ● 被broker接收到只能表示Message已经到达服务器&#xff0c;并不能保证消息一定会投递到目…

【漏洞复现】红海云eHR PtFjk.mob 任意文件上传漏洞

免责声明&#xff1a; 本文内容旨在提供有关特定漏洞或安全漏洞的信息&#xff0c;以帮助用户更好地了解可能存在的风险。公布此类信息的目的在于促进网络安全意识和技术进步&#xff0c;并非出于任何恶意目的。阅读者应该明白&#xff0c;在利用本文提到的漏洞信息或进行相关测…

c++编程(18)——deque的模拟实现(2)容器篇

欢迎来到博主的专栏——c编程 博主ID&#xff1a;代码小豪 文章目录 deque的数据结构deque的构造默认构造填充构造 deque的其他操作deque的插入、删除push_back和push_frontpop_back和pop_frontclear、erase和insert操作 传送门 在上一篇中&#xff0c;我们已经实现了deque最核…

2024/6/16 英语每日一段

Nature has the means--to a degree--to limit the effects of climate change. Intact ecosystems such as forests, grasslands, oceans and peatlands are “carbon sinks”--natural storage systems that remove atmospheric carbon and other greenhouse gases--and are …

【Java03】Java中数组在内存中的机制

1. 内存中的数组 Java中的数组是一种引用类型&#xff0c;数组变量&#xff08;引用&#xff09;和数组元素在内存中是分开的。 Java中的数组变量其实就是指针。 如果想要访问数组元素&#xff0c;只能通过这个数组的引用变量&#xff08;指针&#xff09;来访问。 实际数组对…

华为OD机试 - 多段线数据压缩(Java 2024 D卷 100分)

华为OD机试 2024D卷题库疯狂收录中&#xff0c;刷题点这里 专栏导读 本专栏收录于《华为OD机试&#xff08;JAVA&#xff09;真题&#xff08;D卷C卷A卷B卷&#xff09;》。 刷的越多&#xff0c;抽中的概率越大&#xff0c;每一题都有详细的答题思路、详细的代码注释、样例测…

Hadoop+Spark大数据技术(微课版)总复习

图1 Hadoop开发环境 图2 HDFS 图3 MapReduce 图4 HBase 图5 Scala 图6 Spark 图7 Spark RDD 图8 &#xff08;不考&#xff09; 图9 Spark SQL 图10 Spark Streaming 图11 Spark GraphX 第一章 Hadoop大数据开发环境 hadoop是什么&#xff1f; &#xff08;判断题&#…

数字化转型中的数据资产运营:从数据资产的获取、存储、分析到应用的全流程管理策略

一、引言 随着信息技术的迅猛发展&#xff0c;数字化转型已成为企业提升竞争力、实现可持续发展的关键途径。数据资产作为数字化转型的核心要素&#xff0c;其运营与管理水平直接决定了企业能否在激烈的市场竞争中脱颖而出。本文将从数据资产的获取、存储、分析到应用的全流程…

【尚庭公寓SpringBoot + Vue 项目实战】看房预约管理(十三)

【尚庭公寓SpringBoot Vue 项目实战】看房预约管理&#xff08;十三&#xff09; 文章目录 【尚庭公寓SpringBoot Vue 项目实战】看房预约管理&#xff08;十三&#xff09;1、业务说明2、代码开发2.1、根据条件分页查询预约信息2.2、根据ID更新预约状态 1、业务说明 看房预约…

Sping源码(九)—— Bean的初始化(非懒加载)— Bean的创建方式(Supplier)

序言 目前介绍几种Spring创建对象的方式&#xff0c;其中包括 FactoryBean、Cglib动态代理、自定义BeanPostProcessor&#xff08;InstantiationAwareBeanPostProcessor&#xff09;。 这篇文章会继续扩展Spring中Bean的创建方式——Supplier。 Supplier 先看下我们的Suppli…

【动态规划】| 路径问题之最小路径和 力扣64

&#x1f397;️ 主页&#xff1a;小夜时雨 &#x1f397;️专栏&#xff1a;动态规划 &#x1f397;️如何活着&#xff0c;是我找寻的方向 目录 1. 题目解析2. 代码 1. 题目解析 题目链接: https://leetcode.cn/problems/minimum-path-sum/description/ 这道题目和之前一道…

大模型生成的常见Top-k、Top-p、Temperature参数

参考&#xff1a; https://zhuanlan.zhihu.com/p/669661536 topK&#xff0c;topP https://www.douyin.com/video/7380126984573127945 主要是softmax产生的词表每个词的概率分布后&#xff0c; topK&#xff0c;比如K3&#xff0c;表示采样概率最大的前3个&#xff0c;其他全…

# 梯影传媒T6投影仪刷机方法及一些刷机工具链接

梯影传媒T6投影仪刷机方法及一些刷机工具链接 文章目录 梯影传媒T6投影仪刷机方法及一些刷机工具链接1、安装驱动程序2、备份设备rom【boot、system】3、还原我要刷进设备的rom【system】4、打开开发者模式以便于安装apk5、root设备6、更多好链接&#xff1a; 梯影传媒T6使用的…

经验分享,如何去除文本中的空格

有时候我们需要去掉一窜文本中的空格&#xff0c;这里分享一个好用的免费网站&#xff0c;可实现在线去除 网址&#xff1a;http://www.txttool.com/t/?idMzM4 使用截图&#xff1a;