NCycDB数据库使用（二）：丰度问题

最新推荐文章于 2025-05-28 23:38:48 发布

Ming314!

最新推荐文章于 2025-05-28 23:38:48 发布

阅读量1k

点赞数 18

分类专栏：生物信息学宏基因组/组学分析元素循环分析文章标签：数据库 conda python 经验分享

本文链接：https://blog.youkuaiyun.com/weixin_57129211/article/details/144067542

版权

生物信息学同时被 3 个专栏收录

7 篇文章

订阅专栏

宏基因组/组学分析

6 篇文章

订阅专栏

元素循环分析

2 篇文章

订阅专栏

NCycDB这个数据库我认为比别的数据库方便的就是可以直接得到丰度，那这个数据库进行丰度计算的原理是什么？

这个代码的运行原理是：

1、代码详解

这段代码是一个用于通过 BLAST、DIAMOND 或 USEARCH 对基因组数据进行比对，并根据比对结果计算基因丰度并进行随机抽样的 Perl 脚本。下面是对代码的详细分析：

1.1. 输入参数处理

通过 Getopt::Long 获取用户输入的命令行参数，支持以下参数：

-d ：指定序列文件所在目录。
-m ：指定使用的比对工具（diamond、usearch 或 blast）。
-f ：指定序列文件的类型（支持 .fastq、.fastq.gz、.fasta、.fasta.gz 等格式）。
-s ：指定序列类型（prot 表示蛋白质，nucl 表示核酸）。
-si ：指定样本信息文件，该文件包含每个样本的序列数量。
-rs ：指定随机抽样的大小（默认为最小的样本大小）。
-o ：指定输出文件路径，输出文件包含功能基因的丰度。

1.2. 比对工具选择

根据输入的比对工具（diamond、usearch 或 blast），代码分别进行处理：

DIAMOND：如果选择了 diamond，它会首先创建数据库（使用 diamond makedb），然后使用 blastx 或 blastp 对所有的序列文件进行比对。
USEARCH：如果选择了 usearch，它使用 usearch_global 进行序列比对，并输出到指定文件。
BLAST：如果选择了 blast，它使用 blastall 进行比对，并生成输出文件。

1.3. 基因标识映射 (`id2gene.map`)

使用一个映射文件（data/id2gene.map.2019Jul）将比对结果中的基因ID转换为基因名称。这些基因ID和基因名称通过映射表进行关联，并存储在 %id2gene 哈希表中。

1.4. 计算丰度

根据比对结果计算基因丰度：

对每个比对文件（如 diamond 输出文件），逐行读取比对结果。
比对结果的第一列是查询序列ID，第二列是目标数据库中的ID，第三列是比对的得分等信息。
通过映射表 %id2gene 将目标ID转换为基因名称，并记录每个样本的基因出现次数。

1.5. 样本信息处理

通过读取样本信息文件（sampleinfo），获取每个样本的序列数，并存储在 %size 哈希表中。然后根据最小样本数量或用户指定的随机抽样数量进行抽样。

如果指定了随机抽样大小（-rs），脚本会根据每个样本的大小（样本中序列的数量）进行随机抽样：

每个样本的序列会随机抽取 rs 个序列（默认情况下 rs 是最小样本的序列数）。
随机抽样后，丰度值会根据被抽中的序列数进行调整。

具体的抽样过程在 RandomSampling 函数中实现：

my @array = shuffle( 1 .. $size{$sample} );
@array = @array[ 0 .. $rs - 1 ];  # 保留抽取的部分
@array = grep { $_ <= $sum{$sample} } @array;  # 根据丰度进行筛选

1.6. 随机抽样 (RandomSampling)

在 RandomSampling 子程序中：

使用 List::Util 模块的 shuffle 函数对每个样本的序列进行随机抽样。
抽样的大小由用户指定，默认为样本中最小的序列数。
对每个基因，计算其在抽样后的丰度。

1.7. 输出功能基因丰度

最终，脚本会将每个基因的丰度信息输出到指定的文件（-o），格式如下：

每一行对应一个基因（基因名称），每列对应一个样本的丰度值。
如果某个基因在某个样本中没有出现，则该样本的丰度值为0。

代码的核心流程总结

选择比对工具：根据用户输入选择 diamond、usearch 或 blast。
比对过程：使用选择的工具对序列文件进行比对，并生成比对结果。
基因丰度计算：根据比对结果和基因映射表计算每个基因在不同样本中的丰度。
随机抽样：根据用户指定的随机抽样大小，进行基因丰度的抽样计算。
输出结果：将丰度数据输出到指定的文件。