环境DNA (eDNA) Metabarcoding分析¶

eDNA实验流程：
1. 野外采样 → 采集水样/土壤样本（注意防污染！）
2. 过滤/沉淀 → 用0.45μm滤膜过滤水样，截留DNA
3. DNA提取 → 从滤膜上提取总DNA
4. PCR扩增 → 用通用引物扩增目标marker（如12S、COI）
5. 文库构建 → 加接头和索引，准备上机
6. 测序 → Illumina MiSeq 2×250bp（最常用）
7. 生信分析 → 质控→去噪→物种鉴定→生态分析

常用Marker选择指南：
  鱼类   → 12S rRNA（MiFish引物）
  两栖类 → 12S或16S
  哺乳类 → 12S
  无脊椎 → COI
  真菌   → ITS
  植物   → rbcL / trnL

library(dada2)                           # 加载DADA2

# ===== 设置文件路径 =====
path <- "/data/edna_fastq/"              # FASTQ文件目录
fnFs <- sort(list.files(path, pattern = "_R1_001.fastq.gz", full.names = TRUE))
fnRs <- sort(list.files(path, pattern = "_R2_001.fastq.gz", full.names = TRUE))
sample_names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)  # 提取样本名

# ===== 质量图检查 =====
plotQualityProfile(fnFs[1:2])            # 查看正向reads质量
plotQualityProfile(fnRs[1:2])            # 查看反向reads质量

# ===== 过滤与修剪 =====
filtFs <- file.path(path, "filtered", paste0(sample_names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample_names, "_R_filt.fastq.gz"))

out <- filterAndTrim(
  fnFs, filtFs,                          # 正向：原始→过滤后
  fnRs, filtRs,                          # 反向：原始→过滤后
  truncLen = c(200, 160),                # 截断长度（根据质量图调整）
  maxEE = c(2, 2),                       # 最大期望错误数
  truncQ = 2,                            # 质量值低于2则截断
  maxN = 0,                              # 不允许N碱基
  rm.phix = TRUE,                        # 去除PhiX spike-in
  compress = TRUE,                       # 压缩输出
  multithread = TRUE                     # 多线程
)
# 注意：对于ITS等变长marker，不要使用truncLen截断！

# ===== 学习错误率 =====
errF <- learnErrors(filtFs, multithread = TRUE)  # 学习正向错误模型
errR <- learnErrors(filtRs, multithread = TRUE)  # 学习反向错误模型
plotErrors(errF, nominalQ = TRUE)        # 检查错误模型拟合

# ===== 去噪（核心步骤）=====
dadaFs <- dada(filtFs, err = errF, multithread = TRUE)  # 正向去噪
dadaRs <- dada(filtRs, err = errR, multithread = TRUE)  # 反向去噪

# ===== 合并配对reads =====
mergers <- mergePairs(
  dadaFs, filtFs,                        # 正向结果和文件
  dadaRs, filtRs,                        # 反向结果和文件
  verbose = TRUE                         # 显示详情
)

# ===== 构建ASV表 =====
seqtab <- makeSequenceTable(mergers)     # 序列×样本矩阵
dim(seqtab)                              # 查看维度（样本数 × ASV数）
table(nchar(getSequences(seqtab)))       # 查看ASV长度分布

# ===== 去除嵌合体 =====
seqtab.nochim <- removeBimeraDenovo(
  seqtab,                               # 输入ASV表
  method = "consensus",                  # 共识方法
  multithread = TRUE,                    # 多线程
  verbose = TRUE                         # 显示进度
)
sum(seqtab.nochim) / sum(seqtab)         # 嵌合体去除后保留比例

# ===== 方法1：DADA2内置的assignTaxonomy（适合16S/18S）=====
taxa <- assignTaxonomy(
  seqtab.nochim,                         # ASV序列表
  "/db/MIDORI2_NB_12S.fasta.gz",         # 参考数据库（如MIDORI2的12S）
  multithread = TRUE                     # 多线程
)

# ===== 方法2：BLAST比对（eDNA推荐）=====
# 先导出ASV序列
library(Biostrings)                      # 序列处理包
asv_seqs <- DNAStringSet(getSequences(seqtab.nochim))  # 获取ASV序列
names(asv_seqs) <- paste0("ASV_", seq_along(asv_seqs))  # 命名
writeXStringSet(asv_seqs, "ASV_sequences.fasta")  # 写入FASTA文件

# BLAST比对
blastn \                                 # 核酸BLAST
  -query ASV_sequences.fasta \           # 输入ASV序列
  -db /db/nt \                           # NCBI nt数据库
  -out blast_results.txt \               # 输出结果
  -outfmt "6 qseqid sseqid pident length evalue stitle" \  # 格式化输出
  -max_target_seqs 5 \                   # 每个ASV最多5个匹配
  -evalue 1e-10 \                        # E值阈值
  -num_threads 16                        # 16线程

# 也可以用BOLD数据库（动物条形码专用）
# 在线：https://www.boldsystems.org/index.php/IDS_OpenIdEngine

library(phyloseq)                        # 微生物/eDNA生态分析
library(vegan)                           # 群落生态学
library(ggplot2)                         # 可视化

# 创建phyloseq对象
ps <- phyloseq(
  otu_table(seqtab.nochim, taxa_are_rows = FALSE),  # ASV丰度表
  tax_table(taxa),                       # 分类信息
  sample_data(metadata)                  # 样本元数据（站点、环境因子等）
)

# Alpha多样性
plot_richness(
  ps,                                    # phyloseq对象
  x = "Site",                            # x轴：采样站点
  measures = c("Observed", "Shannon"),   # 观测到的物种数和Shannon指数
  color = "Season"                       # 按季节着色
)

# Beta多样性（物种组成差异）
ord <- ordinate(ps, method = "NMDS", distance = "bray")  # NMDS排序
plot_ordination(
  ps, ord,                               # 数据和排序结果
  color = "Site",                        # 按站点着色
  shape = "Season"                       # 按季节标记
) + geom_point(size = 3)                 # 散点大小

# PERMANOVA：检验组间差异是否显著
adonis2(
  distance(ps, method = "bray") ~        # Bray-Curtis距离
    Site * Season,                       # 站点×季节交互效应
  data = as(sample_data(ps), "data.frame"),  # 元数据
  permutations = 999                     # 置换次数
)

# 物种检出热图
plot_heatmap(
  ps,                                    # phyloseq对象
  taxa.label = "Species",                # 物种标签
  sample.label = "Site"                  # 样本标签
)

# ========== eDNA Metabarcoding速查 ==========

# DADA2核心流程（R）
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(200,160))
errF <- learnErrors(filtFs, multithread=TRUE)
dadaFs <- dada(filtFs, err=errF, multithread=TRUE)
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs)
seqtab <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab)
taxa <- assignTaxonomy(seqtab.nochim, ref_db)

# 物种鉴定（命令行）
blastn -query ASV.fasta -db nt -out results.txt -outfmt 6 -evalue 1e-10

# 常用参考数据库
# 12S(鱼): MIDORI2, MitoFish
# COI(动物): BOLD, MIDORI2
# ITS(真菌): UNITE
# 18S(真核): SILVA, PR2