用 Tax4Fun2 对 16S 微生物组数据进行功能预测

生信菜鸟团

发布于 2020-06-24 16:26:24

4.7K0

发布于 2020-06-24 16:26:24

16S rRNA 基因的高通量测序已被广泛用于研究各种海洋，地表和宿主相关环境中微生物群落的组成和结构。但许多生物学问题更需要我们研究其功能变化，而不仅仅是微生物分类组成。近年来，有不少研究团队开发了几款预测工具，例如 PICRUSt，Tax4Fun，Piphillin，Faprotax 和 paprica。尽管这些工具并不能替代宏基因组测序，但它们仍在一定程度上为我们提供了独特见解。

这些工具的预测能力往往取决于公共数据库中基因组的功能信息。但这些公开的基因组并不一定代表所研究的生态系统中存在的总功能多样性。鉴于公开基因组的数量迅速增加，尤其是一些通过宏基因组分箱得到的数据，未发表的数据，亦或是某些特定环境的基因组信息，所以纳入这些数据可以大大提高功能预测的准确性。

基于这些想法，Tax4Fun 迎来了新的版本升级 Tax4Fun2。Tax4Fun2 是一个快速且易用的 R 包，默认的参考数据集包含 275 个古细菌和 12,002 个细菌基因组。Tax4Fun2 最大的一个升级之处在于可以合并我们的自定义数据集，以增强功能预测的鲁棒性和特异性。

开发团队还将 Tax4Fun2 与 Tax4Fun、PICRUSt 进行比较。结果表明在所有数据集中，Tax4Fun2 的预测性能均优于 PICRUSt 和 Tax4Fun。

Tax4Fun2 预测的功能谱与从宏基因组数据分析得到的功能谱高度相关。尽管 Kelp 相关群落的预测特征与功能特征显著相关，但中位 Spearman 相关系数仅为 0.72，这表明缺少合适的参考数据可能限制了 Tax4Fun2 的性能。

为了解决这个问题，开发团队使用了来自 90 个海带宏基因组的 68 个 MAGs ，构建特定于海带的宏基因组数据集。结果表明这大大提高了预测的准确性（中位 Spearman 相关系数为 0.86），也减少了预测中未使用的序列比例。

此外，添加了自定义的海带数据集后，甚至还可以成功预测那些用参考数据集预测失败的样本。这也证实了添加特定环境参考数据库的优势之处，这种方法也使 Tax4Fun2 更加灵活，并使之在一众类似软件中脱颖而出。

微生物生态学中另一个主要的研究问题是，微生物群落是否包含功能冗余的成分，及其在何种程度上包含功能冗余，从而使之在面对多变的环境时保持生态系统的相对稳定。在 Tax4Fun2 中，针对单个功能通路引入了功能冗余指数（FRI）的概念。我们可以基于包含特定功能的物种比例及其相互之间的系统发育关系来计算 FRI。高 FRI 表示该功能几乎在所有成员中普遍存在，而低 FRI 则表示该功能仅存在于一些密切相关的物种中或仅在一个物种中被检测到，FRI 为 0 表示该功能不存在。Tax4Fun2 会计算相对 FRI（rFRI）和绝对 FRI（aFRI），前者会用环境中的平均系统发生距离进行归一化，后者则用 Tax4Fun2 参考数据中所有原核生物的平均系统发生距离进行归一化。简而言之，rFRI 可用于比较某一项研究中的样本，而 aFRI 允许比较不同生态系统中的功能冗余指数。

为了检测 FRI 的准确性，开发团队模拟了 1,000 个环境，每个环境都包含 100 个原核基因组。接着从每个模拟环境中提取了 16S rRNA 基因序列，将它们以 97％的相似度聚类并计算了 FRI 值。随后根据模拟环境的实际基因组信息将这些值与 FRI 值进行了比较。结果表明， Tax4Fun2 对微生物群落中功能冗余的进行了较好的估计（Spearman 等级相关性> 90％）。此外，用海水样本进行进一步验证。其中六个样本是为赤潮期间采集，而另外三个为正常海水样本，作为参照数据。在参照样本中，将近 7,000 个功能显示出更高的功能冗余，而只有 1,468 个功能在赤潮样本中有更高的冗余度。

这表明在赤潮期间功能冗余发生了很大变化。浮游植物水华通常以基质控制的演替为特征，即在水华期间和之后的不同阶段，不同的细菌进化枝在浮游植物群落中占主导地位。因此，在特定阶段参与某些底物更新的环境成员占主导地位，因此它们的基因组和相关功能将更加冗余，而对于其他所有环境成员，情况则相反。

综上所述，Tax4Fun2 不但具有较高的预测能力，还可通过结合用户自定义数据集来进一步提高其准确性。除此之外，Tax4Fun2 还可以计算特定功能的冗余度，这对于预测在多变环境中定功能丢失的可能性非常重要。

Tax4Fun2 安装及使用

GitHub 地址：https://github.com/bwemheu/Tax4Fun2

Step1. 安装 Tax4Fun2 和下载参考数据库

安装 Tax4Fun2

下载：

wget https://github.com/bwemheu/Tax4Fun2/releases/download/1.1.5/Tax4Fun2_1.1.5.tar.gz

用 install.packages() 安装：

install.packages(pkgs = "Tax4Fun2_1.1.5.tar.gz", repos = NULL, source = TRUE)

加载 Tax4Fun2：

library(Tax4Fun2)

下载参考数据库

在 Tax4Fun2 中我们可用自带的 buildReferenceData() 命令来安装配置参考数据集：

buildReferenceData(path_to_working_directory = ".", use_force = FALSE, install_suggested_packages = TRUE)

参数

•path_to_working_directory 表示数据库安装位置，默认为当前目录；•use_force 表示是否覆盖，默认为 FALSE；•install_suggested_packages 表示是否安装依赖包 ape 和 seqinr，默认为 TURE。

在 Windows 下可能出现数据库已下载但无法解压的情况，我们也可以手动下载然后解压，数据库下载地址：https://owncloud.gwdg.de/index.php/s/pii5Hoitr5SgDqK/download

安装依赖包

此命令将下载 blast（v2.9.0），并将文件夹放在 Tax4Fun2 参考数据库文件夹中。此外，它将测试 R 包 ape 和 seqinr 是否已正确安装。

buildDependencies(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", install_suggested_packages = TRUE)

也可以直接手动下载 blast，然后将 blast_bin 目录放到 Tax4Fun2_ReferenceData_v2 文件夹下，blast 下载地址：ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

下载测试数据

getExampleData(path_to_working_directory = ".")

Step 2. 生成自定义参考数据集

1. 提取 SSU 序列（16S rRNA 和 18S rRNA）

# Option A) 从单个基因组中提取 SSU 序列
extractSSU(genome_file = "OneProkaryoticGenome.fasta", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2")

# Option B) 从多个基因组中提取 SSU 序列
extractSSU(genome_folder = "MoreProkaryoticGenomes", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2")

参数

•genome_file 或 genome_folder：指定单个基因组文件或包含多个基因组的文件夹。必须将多个基因组放在一个文件夹中（一个基因组一个文件），且所有文件为相同的文件扩展名；•file_extension = "fasta"：指定 fasta格式的基因组文件；•path_to_reference_data = "Tax4Fun2_ReferenceData_v2"：指定参考数据集路径。

2. 为原核基因组进行功能注释

# Option A) 对单个基因组进行功能注释
assignFunction(genome_file = "OneProkaryoticGenome.fasta", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", num_of_threads = 1, fast = TRUE)

# Option B) 对多个基因组进行功能注释
assignFunction(genome_folder = "MoreProkaryoticGenomes/", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", num_of_threads = 1, fast = TRUE)

参数

•num_of_threads = 1：设置 diamond 运行的线程数；•fast = TRUE：用默认参数运行 diamond；若设置为 FALSE 则表示运行 diamond --sensitive，可提高灵敏度，但速度要慢得多。

Mac 可能需要自定义 diamond 路径：

# Option A) 对单个基因组进行功能注释
assignFunction(genome_file = "OneProkaryoticGenome.fasta", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", num_of_threads = 1, fast = T, path_to_diamond_binary_mac = "diamond")

# Option B) 对多个基因组进行功能注释
assignFunction(genome_folder = "MoreProkaryoticGenomes/", file_extension = "fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", num_of_threads = 1, fast = T, path_to_diamond_binary_mac = "diamond")

3. 生成自定义参考数据集

提取 SSU 序列并注释功能后，每个基因组都至少包括两个文件：SSU序列文件和功能谱文件。下一步就是用 generateUserData()命令生成自定义参考数据集：

# 1) 根据单个基因组生成参考数据集
generateUserData(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_user_data = ".", name_of_user_data = "User_Ref0", SSU_file_extension = "_16SrRNA.ffn", KEGG_file_extension = "_funPro.txt")

# 2) 根据多个基因组生成参考数据集
generateUserData(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_user_data = "MoreProkaryoticGenomes", name_of_user_data = "User_Ref1", SSU_file_extension = "_16SrRNA.ffn", KEGG_file_extension = "_funPro.txt")

# 3) 根据多个基因组生成参考数据集并用 uclust 去冗余
generateUserDataByClustering(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_user_data = "MoreProkaryoticGenomes", name_of_user_data = "User_Ref2", SSU_file_extension = "_16SrRNA.ffn", KEGG_file_extension = "_funPro.txt", use_force = T)

参数

•path_to_refernce_data：指定参考数据路径；•path_to_user_data：指定自定义数据集路径；•name_of_user_data：自定义数据集名称。

Step 3. 进行功能预测

这里需要把你的 OTU 表格式化成和示例数据中一样的格式。

PlanA 只用默认的参考数据集进行功能预测

# 1. 先进行 blast
runRefBlast(path_to_otus = "KELP_otus.fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR", database_mode = "Ref99NR", use_force = T, num_threads = 6)

# 2. 再进行功能预测，normalize_pathways = FALSE
makeFunctionalPrediction(path_to_otu_table = "KELP_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR", database_mode = "Ref99NR", normalize_by_copy_number = TRUE, min_identity_to_reference = 0.97, normalize_pathways = FALSE)

# 或 normalize_pathways = TRUE

makeFunctionalPrediction(path_to_otu_table = "KELP_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR", database_mode = "Ref99NR", normalize_by_copy_number = TRUE, min_identity_to_reference = 0.97, normalize_pathways = TRUE)

参数

•path_to_otus：指定 OTU 的 fasta 文件；•path_to_otu_table：指定 OTU 矩阵路径；•path_to_refernce_data：指定参考数据库路径；•path_to_temp_folder：指定临时目录，结果将会存放在该目录下；•database_mode = "Ref99NR"：选择参考数据库，可选用 ' Ref99NR 或 Ref100NR，数值表示用 uclust 聚类的阈值（分别为 99% 和 100%）；•num_threads：设置 diamond 运行的线程数；•use_force：表示是否覆盖；•normalize_by_copy_number：是否进行归一化；•normalize_pathways = FALSE ：将每个 KO 的相对丰度关联到它所属的每个通路。将其设置为 true，相对丰度将平均分配到所属的每个通路。

PlanB 结合默认数据库和自定义数据库（未聚类）进行功能预测

# 1. 生成自定义数据集，命名为 KELP1
generateUserData(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_user_data = "KELP_UserData", name_of_user_data = "KELP1", SSU_file_extension = ".ffn", KEGG_file_extension = ".txt")

# 2. 用 include_user_data = TRUE，path_to_user_data，name_of_user_data 三个参数添加自定义数据集进行 blast
runRefBlast(path_to_otus = "KELP_otus.fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR_withUser1", database_mode = "Ref99NR", use_force = T, num_threads = 6, include_user_data = T, path_to_user_data = "KELP_UserData", name_of_user_data = "KELP1")

# 3. 添加自定义数据集进行功能预测
makeFunctionalPrediction(path_to_otu_table = "KELP_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR_withUser1", database_mode = "Ref99NR", normalize_by_copy_number = T, min_identity_to_reference = 0.97, normalize_pathways = F, include_user_data = T, path_to_user_data = "KELP_UserData", name_of_user_data = "KELP1")

参数

•include_user_data = T：包括自定义数据集；•name_of_user_data：自定义数据集的名称；•path_to_user_data ：自定义数据集的路径。

PlanC 结合默认数据库和自定义数据库（聚类）进行功能预测

# 1. 生成自定义参考数据集并用 uclust 去冗余，命名为 KELP2
generateUserDataByClustering(path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_user_data = "KELP_UserData", name_of_user_data = "KELP2", SSU_file_extension = ".ffn", KEGG_file_extension = ".txt", similarity_threshold = 0.99)

# 2. 用 include_user_data = TRUE，path_to_user_data，name_of_user_data 三个参数添加自定义数据集进行 blast
runRefBlast(path_to_otus = "KELP_otus.fasta", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR_withUser2", database_mode = "Ref99NR", use_force = T, num_threads = 6, include_user_data = T, path_to_user_data = "KELP_UserData", name_of_user_data = "KELP2")

# 3. 添加自定义数据集进行功能预测
makeFunctionalPrediction(path_to_otu_table = "KELP_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Kelp_Ref99NR_withUser2", database_mode = "Ref99NR", normalize_by_copy_number = T, min_identity_to_reference = 0.97, normalize_pathways = F, include_user_data = T, path_to_user_data = "KELP_UserData", name_of_user_data = "KELP2")

Step 4. 计算（多）功能冗余指数 FRI

计算 KEGG 功能通路的系统发育分布（高 FRI -> 高冗余，低 FRI -> 低冗余度）：

# 1. 进行 blast
runRefBlast(path_to_otus = "Water_otus.fna", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Water_Ref99NR", database_mode = "Ref99NR", use_force = T, num_threads = 6)

# 2. 计算 FRI
calculateFunctionalRedundancy(path_to_otu_table = "Water_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Water_Ref99NR", database_mode = "Ref99NR", min_identity_to_reference = 0.97)

# New in the latest pre-release (v1.1.6): prevalence_cutoff (see comment on pre-release)
calculateFunctionalRedundancy(path_to_otu_table = "Water_otu_table.txt", path_to_reference_data = "Tax4Fun2_ReferenceData_v2", path_to_temp_folder = "Water_Ref99NR", database_mode = "Ref99NR", min_identity_to_reference = 0.97, prevalence_cutoff = 1.0)

Have Fun！