374_土壤微生物组分析¶

# 激活环境
conda activate qiime2-2024.10

# 步骤1：导入土壤样本（单端测序，部分土壤研究用单端）
qiime tools import \
  --type 'SampleData[SequencesWithQuality]' \
  --input-path soil_manifest.tsv \    # 样本-文件对应表
  --output-path soil_seqs.qza \
  --input-format SingleEndFastqManifestPhred33V2

# 步骤2：去噪（土壤样本质量稍低，截短保守一点）
qiime dada2 denoise-single \
  --i-demultiplexed-seqs soil_seqs.qza \
  --p-trim-left 19 \      # 去掉引物区（515F引物长度）
  --p-trunc-len 230 \     # 截短到230bp
  --p-n-threads 16 \
  --o-table soil_table.qza \
  --o-representative-sequences soil_rep_seqs.qza \
  --o-denoising-stats soil_stats.qza

# 步骤3：系统发育树构建（计算UniFrac必须有树）
qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences soil_rep_seqs.qza \
  --o-alignment soil_aligned_seqs.qza \
  --o-masked-alignment soil_masked_aligned_seqs.qza \
  --o-tree soil_unrooted_tree.qza \
  --o-rooted-tree soil_rooted_tree.qza \
  --p-n-threads 8

# 步骤4：Alpha 多样性（含Faith's PD）
qiime diversity alpha-rarefaction \
  --i-table soil_table.qza \
  --i-phylogeny soil_rooted_tree.qza \
  --p-max-depth 10000 \    # 稀疏化到最小样本量
  --m-metadata-file soil_metadata.tsv \
  --o-visualization soil_alpha_rarefaction.qzv

# 步骤5：Beta 多样性（weighted UniFrac）
qiime diversity core-metrics-phylogenetic \
  --i-phylogeny soil_rooted_tree.qza \
  --i-table soil_table.qza \
  --p-sampling-depth 5000 \    # 稀疏化深度（取最小样本reads数）
  --m-metadata-file soil_metadata.tsv \
  --output-dir soil_core_metrics/

library(vegan)     # 生态学分析包
library(ggplot2)
library(dplyr)

# 加载菌群丰度和土壤理化指标
otu_table  <- read.csv("soil_otu_table.csv", row.names = 1)
soil_chem  <- read.csv("soil_chemistry.csv", row.names = 1)
# soil_chem 包含：pH、有机碳(SOC)、全氮(TN)、含水量等

# 确保样本顺序一致
common_samples <- intersect(rownames(otu_table), rownames(soil_chem))
otu_sub  <- otu_table[common_samples, ]    # 取公共样本
chem_sub <- soil_chem[common_samples, ]

# Bray-Curtis 距离矩阵
bray_dist <- vegdist(otu_sub, method = "bray")

# PERMANOVA：各理化指标解释多少菌群差异？
permanova_result <- adonis2(
  bray_dist ~ pH + SOC + TN + Moisture,   # 同时检验多个环境变量
  data = chem_sub,
  permutations = 999,
  by = "margin"    # 边际检验（控制其他变量后检验每个变量）
)
print(permanova_result)

# 偏Mantel检验：控制地理距离后，pH与菌群相关吗？
geo_dist <- dist(chem_sub[, c("Longitude", "Latitude")])  # 地理距离
ph_dist  <- dist(chem_sub$pH)                              # pH差异距离

partial_mantel <- mantel.partial(
  bray_dist,   # 菌群距离
  ph_dist,     # pH距离
  geo_dist,    # 控制的地理距离
  method = "spearman",
  permutations = 999
)
print(partial_mantel)   # 看控制地理效应后pH是否还显著