47. 微生物组与宿主互作分析（Microbiome-Host Interaction Analysis）¶

# NetCoMi 网络构建示例
library(NetCoMi)  # 加载 NetCoMi 包

# 构建网络：使用 SParCC 计算相关性
net_result <- netConstruct(
  data = otu_table,           # 输入 OTU/ASV 丰度表
  measure = "sparcc",         # 用 SParCC 方法计算相关（适合组成型数据）
  measurePar = list(          # SParCC 参数
    iter = 20,                # 迭代次数
    inner_iter = 10           # 内部迭代次数
  ),
  filtTax = "highestFreq",    # 过滤策略：保留高频taxa
  filtTaxPar = list(highestFreq = 50),  # 保留前50个最常见的taxa
  sparsMethod = "threshold",  # 稀疏化方法：阈值法
  thresh = 0.3                # 相关系数阈值：绝对值>0.3的保留
)

# 分析网络属性
net_props <- netAnalyze(
  net_result,                 # 网络构建结果
  centrLCC = TRUE,            # 在最大连通分量上计算中心性
  clustMethod = "cluster_fast_greedy"  # 社区检测算法
)

# 可视化网络
plot(net_props,
     nodeColor = "cluster",   # 节点颜色按社区着色
     nodeSize = "eigenvector", # 节点大小按特征向量中心性
     labelScale = FALSE,      # 标签不缩放
     cexLabels = 0.8          # 标签字体大小
)

人类基因组 SNP → 影响宿主表型（如免疫因子、黏液分泌、pH）
                    ↓
               塑造肠道微环境
                    ↓
             决定哪些菌能定植、丰度多少

# mbQTL 分析简化流程（需要大样本量，通常 >1000 人）

# 1. 准备宿主基因型数据（PLINK 格式）
plink --bfile host_genotype \   # 输入基因型文件
  --maf 0.05 \                  # 最小等位基因频率过滤
  --geno 0.02 \                 # 基因型缺失率过滤
  --mind 0.02 \                 # 个体缺失率过滤
  --hwe 1e-6 \                  # Hardy-Weinberg 平衡检验
  --make-bed \                  # 输出 bed 格式
  --out host_qc                 # 输出文件名

# 2. 准备菌群丰度数据（CLR 转换后）
# 在 R 中用 MaAsLin2 或自定义脚本做 CLR 转换

# 3. 运行关联分析（通常用线性混合模型）
# 常用工具：PLINK2、BOLT-LMM、GEMMA
plink2 --bfile host_qc \        # QC后的基因型
  --pheno microbiome_clr.txt \  # CLR转换后的菌群丰度（每个菌属一列）
  --glm \                       # 广义线性模型
  --covar covariates.txt \      # 协变量（年龄、性别、BMI、前10个PC）
  --out mbqtl_results           # 输出结果

# 4. 多重检验校正（Bonferroni 或 FDR）
# SNP数 × 菌属数 = 巨大的检验次数，需严格校正

# 宏转录组分析流程

# 1. 质控（去接头、去低质量）
fastp -i raw_R1.fq.gz \         # 输入原始正向reads
  -I raw_R2.fq.gz \             # 输入原始反向reads
  -o clean_R1.fq.gz \           # 输出质控后正向reads
  -O clean_R2.fq.gz             # 输出质控后反向reads

# 2. 去除宿主 RNA（非常关键！宿主RNA通常占 >90%）
bowtie2 -x human_genome \      # 人类参考基因组索引
  -1 clean_R1.fq.gz \          # 质控后的reads
  -2 clean_R2.fq.gz \
  --un-conc-gz microbial.fq.gz \  # 未比对上的reads = 微生物reads
  -S host_aligned.sam           # 比对上的reads = 宿主reads

# 3. 去除 rRNA（微生物 rRNA 占大量，不是我们要的 mRNA）
sortmerna --ref rRNA_database \  # rRNA参考数据库
  --reads microbial.fq.gz \     # 微生物reads
  --aligned rRNA_reads \        # rRNA reads（丢弃）
  --other mRNA_reads \          # mRNA reads（保留！这才是目标）
  --fastx                       # 输出fastq格式

# 4. 功能注释（用 HUMAnN 做通路分析）
humann --input mRNA_reads.fq.gz \    # 输入mRNA reads
  --output metatranscriptome_out \   # 输出目录
  --nucleotide-database chocophlan \ # 核酸数据库
  --protein-database uniref         # 蛋白数据库

# 5. 计算 RNA/DNA 比值 = 基因活性指标
# RNA丰度 / DNA丰度 = 转录活性
# 比值 > 1 说明该基因/通路被"积极激活"

膳食纤维 → 肠道菌群发酵 → SCFA → 多重效应：
  ├─ 丁酸 → 供能给肠上皮细胞 → 维持肠道屏障
  ├─ 丙酸 → 激活 GPR41/43 受体 → 刺激 GLP-1/PYY 分泌 → 调节血糖
  └─ 乙酸 → 进入血液 → 影响肝脏脂肪合成和外周代谢

菌群的 BSH 酶（胆盐水解酶）→ 解偶联胆汁酸
                ↓
      次级胆汁酸（DCA、LCA）
                ↓
    激活 FXR 受体 → 调节肝脏糖异生
    激活 TGR5 受体 → 促进 GLP-1 分泌 → 改善胰岛素分泌

# MaAsLin2 完整分析示例
library(Maaslin2)  # 加载包（Bioconductor安装：BiocManager::install("Maaslin2")）

# 输入数据准备
# features：行=样本，列=微生物特征（物种/通路等）
# metadata：行=样本，列=临床变量

result <- Maaslin2(
  input_data = "species_abundance.tsv",  # 微生物丰度表（TSV格式）
  input_metadata = "metadata.tsv",       # 样本元数据表
  output = "maaslin2_output",            # 输出目录

  # 模型参数
  fixed_effects = c("T2D_status",        # 固定效应：T2D状态（主要研究变量）
                     "age", "BMI",        # 需要校正的混杂因素
                     "sex"),
  random_effects = c("subject_id"),       # 随机效应：受试者ID（用于纵向数据）

  # 统计方法
  analysis_method = "LM",                # 分析方法：LM=线性模型（也可选CPLM, NEGBIN, ZINB）
  normalization = "TSS",                 # 标准化：TSS=总和缩放（相对丰度）
  transform = "LOG",                     # 转换：对数转换

  # 过滤
  min_abundance = 0.0001,               # 最小丰度阈值（去掉极低丰度的菌）
  min_prevalence = 0.1,                  # 最小流行度（至少10%样本中出现）

  # 多重检验校正
  correction = "BH",                     # Benjamini-Hochberg FDR 校正
  significance_threshold = 0.25,         # q-value 阈值（MaAsLin2 默认 0.25）

  # 参考水平
  reference = c("T2D_status,control")    # 参考组：健康对照
)

# 输出文件说明：
# significant_results.tsv → 显著关联的特征
# all_results.tsv → 所有检验结果
# figures/ → 散点图（每个显著关联一张）

# SParCC Python 使用示例（fastspar 是 C++ 加速版）
# 安装：conda install -c bioconda fastspar

# 命令行运行 FastSpar（SParCC 的快速实现）
# 1. 计算相关系数
fastspar \
  --otu_table otu_table.tsv \       # 输入OTU表（tab分隔）
  --correlation cor_matrix.tsv \    # 输出相关系数矩阵
  --covariance cov_matrix.tsv \     # 输出协方差矩阵
  --iterations 50 \                 # 迭代次数（越多越准，默认50）
  --threads 4                       # 线程数

# 2. Bootstrap 计算 p 值
mkdir bootstrap_counts              # 创建bootstrap目录
fastspar_bootstrap \
  --otu_table otu_table.tsv \       # 原始OTU表
  --number 1000 \                   # bootstrap次数（推荐>=1000）
  --prefix bootstrap_counts/boot    # 输出前缀

# 3. 对每个bootstrap样本计算相关系数
parallel fastspar \
  --otu_table {} \
  --correlation bootstrap_cor/cor_{/} \
  --covariance bootstrap_cor/cov_{/} \
  --iterations 10 \
  ::: bootstrap_counts/boot_*.tsv

# 4. 计算 p 值
fastspar_pvalues \
  --otu_table otu_table.tsv \        # 原始OTU表
  --correlation cor_matrix.tsv \     # 原始相关系数
  --prefix bootstrap_cor/cor_ \      # bootstrap相关系数前缀
  --permutations 1000 \              # 排列次数
  --outfile pvalues.tsv              # 输出p值矩阵

# NetCoMi 两组网络比较示例（T2D vs 健康人）
library(NetCoMi)

# 构建两组网络
net_compare <- netConstruct(
  data = otu_t2d,               # T2D组的OTU表
  data2 = otu_healthy,          # 健康组的OTU表
  measure = "sparcc",           # 使用SParCC计算相关
  filtTax = "highestFreq",      # 过滤：保留高频taxa
  filtTaxPar = list(highestFreq = 30),  # 保留前30个
  sparsMethod = "threshold",    # 阈值稀疏化
  thresh = 0.3                  # 相关阈值
)

# 分析网络差异
net_anal <- netAnalyze(
  net_compare,
  centrLCC = TRUE,              # 在最大连通分量上计算
  clustMethod = "cluster_fast_greedy"
)

# 比较两个网络的差异（置换检验）
net_comp <- netCompare(
  net_anal,
  permTest = TRUE,              # 做置换检验
  nPerm = 1000,                 # 置换次数
  seed = 42                     # 随机种子
)

summary(net_comp)               # 查看差异显著性
# 输出：哪些网络属性（连通性、模块化、Hub节点）在两组间显著不同

# MMUPHin 批次校正 + 荟萃分析示例
library(MMUPHin)  # 安装：BiocManager::install("MMUPHin")

# 1. 批次效应校正（不同研究的系统差异）
fit_adjust <- adjust_batch(
  feature_abd = abundance_matrix,  # 微生物丰度矩阵（行=特征，列=样本）
  batch = "study_id",              # 批次变量：研究来源
  covariates = "disease_status",   # 需要保留的生物学差异
  data = sample_metadata           # 样本元数据
)
# fit_adjust$feature_abd_adj 就是校正后的丰度矩阵

# 2. 荟萃分析：跨队列发现与T2D相关的菌
fit_meta <- lm_meta(
  feature_abd = abundance_matrix,  # 丰度矩阵
  batch = "study_id",              # 批次变量
  exposure = "T2D_status",         # 暴露变量（T2D/对照）
  covariates = c("age", "sex", "BMI"),  # 协变量
  data = sample_metadata           # 元数据
)

# 查看结果：哪些菌在多个队列中一致地与T2D相关
fit_meta$meta_fits  # 荟萃分析结果，包含效应量和p值

# CIBERSORT 反卷积 + 菌群关联分析流程

# 第一步：上传 bulk RNA-seq 数据到 CIBERSORTx 网站
# 网址：https://cibersortx.stanford.edu/
# 选择 LM22 签名矩阵（22种免疫细胞）
# 下载结果：每个样本的22种免疫细胞比例

# 第二步：在 R 中关联免疫细胞比例与菌群丰度
library(Maaslin2)

# 将免疫细胞比例作为 metadata 的一部分
immune_microbe_result <- Maaslin2(
  input_data = "species_abundance.tsv",  # 菌群丰度
  input_metadata = "immune_metadata.tsv", # 包含免疫细胞比例的元数据
  output = "immune_microbe_output",
  fixed_effects = c("M1_Macrophages",    # 关注的免疫细胞类型
                     "CD8_T_cells",
                     "Tregs",
                     "age", "sex"),        # 校正变量
  normalization = "TSS",
  transform = "LOG"
)

# 或者用 Spearman 相关做快速探索
# 菌群丰度矩阵 vs 免疫细胞比例矩阵
cor_result <- cor(
  microbe_abundance,    # 菌群丰度矩阵
  immune_fractions,     # 免疫细胞比例矩阵
  method = "spearman"   # Spearman相关（适合非正态数据）
)

# 可视化：热图
library(pheatmap)
pheatmap(cor_result,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "Microbiome-Immune Cell Correlation")

T2D 患者肠道菌群失调（Dysbiosis）
    │
    ├── SCFA产生菌减少 → 丁酸↓ → 肠道屏障受损 → 细菌 LPS 入血
    │                                              │
    │                          ┌────────────────────┘
    │                          ↓
    ├── LPS 激活 TLR4 → NF-κB 通路 → 慢性低度炎症（TNF-α, IL-6↑）
    │                                              │
    │                                              ↓
    ├── 炎症因子干扰胰岛素信号通路 → IRS-1 丝氨酸磷酸化 → 胰岛素抵抗
    │
    ├── 胆汁酸代谢紊乱 → FXR/TGR5 信号异常 → GLP-1 分泌减少
    │
    └── 支链氨基酸（BCAA）代谢改变 → mTOR 过度激活 → β细胞功能障碍

你要研究什么？
│
├── 菌群与临床指标的关联 → MaAsLin2
│     └── 有混杂因素？ → 必须用 MaAsLin2（不能用 LEfSe）
│
├── 菌群之间的相互关系 → SParCC + NetCoMi
│     └── 要比较两组网络？ → NetCoMi（netCompare）
│
├── 多个队列数据整合 → MMUPHin
│     └── 有批次效应？ → 先用 adjust_batch 校正
│
├── 菌群与免疫细胞的关系 → CIBERSORTx + MaAsLin2
│
├── 因果推断 → 孟德尔随机化（TwoSampleMR）
│
└── 多组学整合 → MOFA2 / DIABLO

第1天：理解互作概念 + SCFA/胆汁酸通路
    ↓
第2天：跑通 MaAsLin2（菌群-表型关联）
    ↓
第3天：学会 SParCC/FastSpar（相关网络）
    ↓
第4天：用 NetCoMi 构建和比较网络
    ↓
第5天：理解 CIBERSORTx 反卷积 + 免疫-菌群关联
    ↓
第6天：MMUPHin 多队列整合
    ↓
第7天：整合到该 T2D 项目叙事中