658 宏基因组功能冗余分析¶

# 基于GCN网络计算功能冗余
import numpy as np  # 数值计算
import pandas as pd  # 数据处理
from scipy.spatial.distance import jaccard  # Jaccard距离

# 读取基因存在/缺失矩阵（行=物种，列=基因/KO）
gene_matrix = pd.read_csv("species_gene_matrix.csv", index_col=0)  # 物种-基因矩阵
abundance = pd.read_csv("species_abundance.csv", index_col=0)  # 物种丰度表

def calculate_genome_fr(gene_matrix, abundance_vector):
    """计算基于基因组的功能冗余指数"""
    n_species = gene_matrix.shape[0]  # 物种数量
    n_genes = gene_matrix.shape[1]  # 基因数量

    # 计算每个基因被多少物种共享（加权）
    gene_redundancy = []  # 存储每个基因的冗余度
    for gene in gene_matrix.columns:  # 遍历每个基因
        carriers = gene_matrix[gene] > 0  # 携带该基因的物种
        n_carriers = carriers.sum()  # 携带物种数
        if n_carriers > 0:
            # 加权冗余：考虑携带物种的丰度
            weighted_r = abundance_vector[carriers].sum()  # 携带物种的总丰度
            gene_redundancy.append(weighted_r)
        else:
            gene_redundancy.append(0)

    # 总FR = 所有基因冗余度的均值
    fr_index = np.mean(gene_redundancy)  # 群落整体功能冗余指数
    return fr_index

# 计算每个样本的FR
for sample in abundance.columns:  # 遍历每个样本
    fr = calculate_genome_fr(gene_matrix, abundance[sample])  # 计算FR
    print(f"样本 {sample}: FR = {fr:.4f}")  # 输出结果

# 基于相对熵的功能冗余量化
from scipy.stats import entropy  # 信息熵函数
from scipy.special import rel_entr  # 相对熵（KL散度）

def single_trait_fr(trait_profile, species_abundance):
    """
    单性状水平的功能冗余量化
    trait_profile: 每个物种对某性状的贡献向量
    species_abundance: 物种相对丰度向量
    """
    # 归一化性状贡献
    trait_contrib = trait_profile * species_abundance  # 加权贡献
    trait_contrib = trait_contrib / trait_contrib.sum()  # 归一化为概率分布

    # 均匀分布作为最大冗余参考
    n_contributing = (trait_profile > 0).sum()  # 贡献物种数
    uniform = np.ones(n_contributing) / n_contributing  # 均匀分布

    # 功能冗余 = 1 - 归一化KL散度
    # KL散度越小 = 贡献越均匀 = 冗余越高
    contributing_idx = trait_profile > 0  # 有贡献的物种
    kl_div = np.sum(rel_entr(trait_contrib[contributing_idx], uniform))  # KL散度
    max_kl = np.log(n_contributing)  # 最大可能KL散度
    fr = 1 - (kl_div / max_kl) if max_kl > 0 else 0  # 归一化FR指数

    return fr  # 0=无冗余，1=完全冗余

# 示例
species_abund = np.array([0.3, 0.2, 0.15, 0.15, 0.1, 0.1])  # 6个物种的丰度
butyrate_trait = np.array([1, 1, 0, 1, 0, 1])  # 哪些物种能产丁酸（1=能，0=不能）
fr_butyrate = single_trait_fr(butyrate_trait, species_abund)
print(f"丁酸产生功能冗余度: {fr_butyrate:.3f}")  # 输出FR值

# 群落内FR和群落间FR
def within_community_fr(gene_matrix, abundance):
    """群落内功能冗余：同一群落中不同物种共享同一功能"""
    n_genes = gene_matrix.shape[1]  # 基因总数
    fr_scores = []  # 存储每个基因的FR

    for gene in gene_matrix.columns:
        carriers = gene_matrix[gene] > 0  # 携带该基因的物种
        n_carriers = carriers.sum()  # 携带物种数
        if n_carriers > 1:  # 至少2个物种共享才算冗余
            fr_scores.append(n_carriers - 1)  # 冗余度 = 携带数-1
        else:
            fr_scores.append(0)  # 无冗余

    return np.mean(fr_scores)  # 平均群落内FR

def between_community_fr(gene_matrices, labels):
    """群落间功能冗余：不同群落中不同物种执行相同功能"""
    # gene_matrices: 多个群落的基因矩阵列表
    all_genes = set()
    for gm in gene_matrices:
        all_genes.update(gm.columns)

    shared_functions = 0  # 共享功能计数
    total_functions = len(all_genes)  # 总功能数

    for gene in all_genes:
        communities_with_gene = 0  # 拥有该功能的群落数
        for gm in gene_matrices:
            if gene in gm.columns and (gm[gene] > 0).any():
                communities_with_gene += 1
        if communities_with_gene > 1:
            shared_functions += 1  # 跨群落共享

    return shared_functions / total_functions if total_functions > 0 else 0

# 使用COBRApy进行代谢冗余分析
# pip install cobra
import cobra  # 基于约束的代谢建模

def metabolic_fr_analysis(model_files, target_metabolite):
    """
    基于GSMM分析代谢功能冗余
    model_files: 各物种的代谢模型文件列表
    target_metabolite: 目标代谢物ID
    """
    producers = []  # 能产生目标代谢物的物种

    for model_file in model_files:
        model = cobra.io.read_sbml_model(model_file)  # 读取代谢模型

        # 检查是否能产生目标代谢物
        for rxn in model.reactions:  # 遍历所有反应
            metabolite_ids = [m.id for m in rxn.metabolites]  # 反应涉及的代谢物
            if target_metabolite in metabolite_ids:
                # FBA检验该反应是否可行
                with model:
                    model.objective = rxn  # 设置优化目标
                    sol = model.optimize()  # 求解
                    if sol.status == 'optimal' and sol.objective_value > 1e-6:
                        producers.append(model_file)  # 记录能产生该代谢物的物种
                        break

    fr = len(producers) / len(model_files)  # FR = 能产生的物种比例
    print(f"目标代谢物 {target_metabolite}:")
    print(f"  能产生的物种数: {len(producers)}/{len(model_files)}")
    print(f"  功能冗余度: {fr:.2f}")
    return producers, fr

# R语言：FR与群落稳定性关联分析
library(vegan)    # 群落生态学分析
library(ggplot2)  # 绑图

# 读取数据
otu <- read.csv("otu_table.csv", row.names = 1)  # OTU丰度表
func <- read.csv("function_table.csv", row.names = 1)  # 功能丰度表

# 计算物种多样性
alpha_div <- diversity(t(otu), index = "shannon")  # Shannon多样性

# 计算功能冗余（简化版：物种数/功能数）
n_species <- apply(otu > 0, 2, sum)  # 每个样本的物种数
n_functions <- apply(func > 0, 2, sum)  # 每个样本的功能数
fr_index <- n_species / n_functions  # 简化FR指数

# 计算群落稳定性（时间序列变异系数的倒数）
# stability = 1 / CV(abundance)
stability <- 1 / apply(otu, 1, function(x) sd(x)/mean(x))

# 可视化FR与稳定性的关系
df <- data.frame(FR = fr_index, Stability = stability, Diversity = alpha_div)
ggplot(df, aes(x = FR, y = Stability, color = Diversity)) +
  geom_point(size = 3) +  # 散点
  geom_smooth(method = "lm") +  # 线性拟合
  labs(x = "功能冗余指数", y = "群落稳定性", color = "Shannon多样性") +
  theme_minimal()

# 功能冗余分析工作流
1. 物种注释 → Kraken2/MetaPhlAn
2. 功能注释 → HUMAnN3/eggNOG-mapper
3. 构建物种-功能矩阵
4. 计算群落内FR（GCN/信息熵方法）
5. 比较组间FR差异（Wilcoxon/PERMANOVA）
6. FR与多样性/稳定性关联分析

# 关键R包/Python库
R: vegan, phyloseq, microbiome
Python: scikit-bio, cobra, numpy/scipy
专用工具: coda4microbiome (2024)