环境宏基因组分析（eDNA）¶
一句话概述¶
环境DNA（eDNA）技术通过从土壤、水体等环境样品中提取DNA，设计特异性引物进行扩增子测序或宏基因组测序，利用OBITools/DADA2等工具处理数据实现物种鉴定，是生态监测和生物多样性调查的革命性方法。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| eDNA概念 | 环境中游离的DNA片段用于物种检测 | - |
| 引物设计 | 通用引物/物种特异引物选择 | ecoPrimers, PrimerMiner |
| 文库构建 | 双索引策略、标签跳跃控制 | - |
| 序列处理 | 质控→去引物→去噪/聚类 | OBITools, Cutadapt, DADA2 |
| 分类注释 | ASV/OTU的物种鉴定 | BLAST, ecotag, RDP classifier |
| 假阳性控制 | 去污染、阴性对照检查 | microDecon, decontam |
| 定量分析 | 读数与生物量的关系 | - |
| 群落分析 | 多样性、组成、占域模型 | vegan, phyloseq, occupancy |
各步骤详解¶
第一步：理解eDNA技术原理¶
白话解释： 生物在环境中生活时会不断释放DNA到周围环境中（通过脱落皮肤细胞、排泄物、分泌物等）。这些漂浮在水中或存在于土壤中的DNA碎片就是"环境DNA"。通过采集水/土样本、提取其中的DNA、用特定引物扩增目标基因片段并测序，就能知道这个环境中有哪些生物存在——不需要直接看到或捕获它们。
技术细节： - eDNA半衰期：水中通常1-14天（温度、UV、pH影响） - 常用标记基因：COI(动物)、12S/16S rRNA(鱼类/脊椎动物)、ITS(真菌)、rbcL/matK(植物) - metabarcoding vs species-specific qPCR - 主要挑战：PCR偏好、引物覆盖度、参考数据库不完整、定量不精确
代码示例：
# eDNA工作流程概览
# 1. 野外采样(水过滤/土壤取样)
# 2. DNA提取(DNeasy PowerWater/PowerSoil)
# 3. PCR扩增(目标基因片段)
# 4. 文库构建(双索引)
# 5. 测序(Illumina MiSeq/NovaSeq)
# 6. 生信分析(本教程重点)
# 常用引物组合
# 鱼类12S: MiFish-U (Miya et al. 2015) ~170bp
# 脊椎动物12S: teleo (Valentini et al. 2016) ~60bp
# 动物COI: mlCOIintF/jgHCO2198 ~313bp
# 真菌ITS: ITS1F/ITS2 ~250bp
# 细菌16S: 515F/806R ~250bp
第二步：引物设计与评估¶
白话解释： 引物设计是eDNA研究成败的关键。理想引物需要：1）能扩增目标类群的所有物种（通用性好）；2）不扩增非目标类群（特异性好）；3）扩增片段能区分物种（分辨率高）；4）片段不能太长（eDNA已降解）。
代码示例：
# 使用ecoPrimers设计引物（OBITools套件）
# 需要先准备参考数据库
obiconvert --fasta reference_sequences.fasta \
--ecopcrdb-output=reference_db
# 运行ecoPrimers
ecoPrimers -d reference_db \
-l 50 -L 200 \ # 扩增片段长度范围
-e 3 \ # 允许3个错配
-t Vertebrata \ # 目标分类群
--taxonomy taxonomy_db \
-o candidate_primers.txt
# 使用PrimerMiner评估引物覆盖度
# R包
Rscript << 'EOF'
library(PrimerMiner)
# 从GenBank下载参考序列
batch_download(
taxon = "Actinopterygii",
gene = "12S",
output = "fish_12S_refs.fasta")
# 评估引物覆盖度
evaluate_primers(
primer_F = "GTCGGTAAAACTCGTGCCAGC", # MiFish-U-F
primer_R = "CATAGTGGGGTATCTAATCCCAGTTTG", # MiFish-U-R
ref_sequences = "fish_12S_refs.fasta",
mismatch_threshold = 3)
EOF
# 使用ecoPCR模拟in silico PCR
ecoPCR -d reference_db \
-e 3 \ # 允许3个错配
-l 100 -L 300 \ # 片段长度
GTCGGTAAAACTCGTGCCAGC CATAGTGGGGTATCTAATCCCAGTTTG \
> ecopcr_results.txt
# 统计覆盖物种数
awk -F'\t' '{print $2}' ecopcr_results.txt | sort -u | wc -l
第三步：OBITools数据处理流程¶
白话解释： OBITools是专为eDNA metabarcoding设计的数据处理套件。它从原始reads开始：配对合并→按样品分拣→去重→去噪→分类注释。相比其他流程，OBITools更注重对短片段和降解DNA的处理。
代码示例：
# 安装OBITools3（当前版本3.0.1b26）
# OBITools3 不在 bioconda 中（bioconda只有旧版OBITools1），需通过 pip 安装:
pip install OBITools3
# 或从源码安装：git clone https://git.metabarcoding.org/obitools/obitools3.git
# === OBITools3 完整流程 ===
# 1. 导入数据
obi import --fastq-input sample_R1.fastq.gz reads/reads1
obi import --fastq-input sample_R2.fastq.gz reads/reads2
# 2. 配对合并
obi alignpairedend -R reads/reads2 reads/reads1 reads/aligned
# 3. 去除未成功合并的reads（只保留mode=alignment的序列）
obi grep -a mode:alignment reads/aligned reads/goodreads
# 4. 按样品拆分（根据tag/index）
obi ngsfilter -t ngsfilter.txt -u reads/unidentified \
reads/goodreads reads/identified
# ngsfilter.txt格式:
# experiment sample forward_tag reverse_tag forward_primer reverse_primer
# exp1 sample1 ACGT TGCA GTCGGTAAAACTCGTGCCAGC CATAGTGGGGTATCTAATCCCAGTTTG
# 5. 去重复（合并相同序列并计数）
obi uniq -m sample reads/identified reads/dereplicated
# 6. 去噪（去除PCR/测序错误）
obi annotate -k COUNT -k MERGED_sample reads/dereplicated reads/annotated
obi grep -p "sequence['COUNT'] >= 10" reads/annotated reads/filtered
# 7. 去除单片段和短序列
obi grep -p "len(sequence) >= 80" reads/filtered reads/length_filtered
# 8. 清理PCR错误（obiclean）
obi clean -s MERGED_sample -r 0.05 reads/length_filtered reads/cleaned
# 9. 分类注释
# 先导入参考序列和分类信息
obi import --fasta-input EMBL_12S.fasta refs/db
# 建立参考数据库（-t 为相似度阈值）
obi build_ref_db -t 0.97 --taxonomy refs/taxonomy refs/db refs/db_indexed
# 使用ecotag进行分类（-m 为最小匹配度，-R 为参考数据库）
obi ecotag -m 0.97 --taxonomy refs/taxonomy -R refs/db_indexed reads/cleaned reads/assigned
# 10. 导出结果
obi export --tab-output reads/assigned > final_results.tsv
# 结果包含: 序列、计数、分类(科/属/种)、相似度、样品分布
第四步：DADA2处理流程（替代方案）¶
白话解释： DADA2是另一个主流的eDNA数据处理方案，生成精确的ASV（扩增子序列变体）而非传统OTU。ASV有单碱基分辨率，且在不同数据集间可直接比较。
代码示例：
library(dada2)
library(Biostrings)
# 设置路径
path <- "raw_fastq/"
fnFs <- sort(list.files(path, pattern="_R1.fastq.gz", full.names=TRUE))
fnRs <- sort(list.files(path, pattern="_R2.fastq.gz", full.names=TRUE))
sample.names <- gsub("_R1.fastq.gz", "", basename(fnFs))
# 1. 质量检查
plotQualityProfile(fnFs[1:4])
plotQualityProfile(fnRs[1:4])
# 2. 去引物（使用Cutadapt）
system2("cutadapt", args=c(
"-g", "GTCGGTAAAACTCGTGCCAGC", # Forward primer
"-G", "CATAGTGGGGTATCTAATCCCAGTTTG", # Reverse primer
"--discard-untrimmed",
"-o", "trimmed/sample_R1.fastq.gz",
"-p", "trimmed/sample_R2.fastq.gz",
"raw_fastq/sample_R1.fastq.gz",
"raw_fastq/sample_R2.fastq.gz"))
# 3. 过滤和修剪
filtFs <- file.path("filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path("filtered", paste0(sample.names, "_R_filt.fastq.gz"))
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
truncLen = c(140, 140), # 根据质量图确定
maxN = 0,
maxEE = c(2, 2),
truncQ = 2,
rm.phix = TRUE,
compress = TRUE,
multithread = TRUE)
# 4. 学习错误率
errF <- learnErrors(filtFs, multithread=TRUE, randomize=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE, randomize=TRUE)
plotErrors(errF, nominalQ=TRUE)
# 5. 去噪（核心步骤）
dadaFs <- dada(filtFs, err=errF, multithread=TRUE, pool="pseudo")
dadaRs <- dada(filtRs, err=errR, multithread=TRUE, pool="pseudo")
# 6. 合并配对
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)
# 7. 构建ASV表
seqtab <- makeSequenceTable(mergers)
dim(seqtab) # samples × ASVs
# 8. 去嵌合体
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus",
multithread=TRUE, verbose=TRUE)
cat("保留比例:", sum(seqtab.nochim)/sum(seqtab), "\n")
# 9. 分类注释
# 使用自定义参考数据库（eDNA需要针对目标类群的数据库）
taxa <- assignTaxonomy(seqtab.nochim,
"reference_db/fish_12S_taxonomy.fa.gz",
multithread=TRUE, minBoot=80)
# 补充种水平注释
taxa <- addSpecies(taxa, "reference_db/fish_12S_species.fa.gz")
# 10. 导出
write.csv(seqtab.nochim, "asv_table.csv")
write.csv(taxa, "taxonomy_table.csv")
第五步：去污染与阴性对照处理¶
白话解释： eDNA研究中污染是最大的威胁。阴性对照（field blank/extraction blank/PCR blank）中不应该有任何检测，如果有，说明存在污染。需要用统计方法判断哪些检测是真阳性，哪些是污染。
代码示例：
# 使用decontam包去污染
library(decontam)
library(phyloseq)
# 创建phyloseq对象
ps <- phyloseq(
otu_table(seqtab.nochim, taxa_are_rows=FALSE),
tax_table(taxa),
sample_data(metadata))
# 标记阴性对照
sample_data(ps)$is_neg <- sample_data(ps)$sample_type == "negative_control"
# 方法1: 频率法（需要DNA浓度信息）
contamdf.freq <- isContaminant(ps, method="frequency",
conc="quant_reading")
# 方法2: 流行率法（基于阴性对照中的检出率）
# threshold默认0.1；0.5更宽松，会标记更多为污染
contamdf.prev <- isContaminant(ps, method="prevalence",
neg="is_neg", threshold=0.5)
# 查看污染ASV
table(contamdf.prev$contaminant)
contam_asvs <- rownames(contamdf.prev[contamdf.prev$contaminant == TRUE, ])
# 去除污染ASV
ps_clean <- prune_taxa(!taxa_names(ps) %in% contam_asvs, ps)
# 去除阴性对照样本
ps_clean <- subset_samples(ps_clean, sample_type != "negative_control")
# 额外过滤：设置最小读数阈值
# eDNA研究中通常设置per-sample最小读数
otu <- as.data.frame(otu_table(ps_clean))
otu[otu < 5] <- 0 # 每个样本中<5 reads的检测视为不可靠
otu_table(ps_clean) <- otu_table(otu, taxa_are_rows=FALSE)
第六步：物种鉴定与生态分析¶
白话解释： 最终目标是知道环境中有哪些物种。将ASV序列与参考数据库比对，获得物种鉴定结果后，进行多样性分析、群落组成比较、占域模型等生态学分析。
代码示例：
# 物种检测结果总结
library(phyloseq)
library(vegan)
library(ggplot2)
# 按分类等级汇总
species_list <- unique(tax_table(ps_clean)[, "Species"])
species_list <- species_list[!is.na(species_list)]
cat("检测到物种数:", length(species_list), "\n")
# Alpha多样性
alpha_div <- estimate_richness(ps_clean, measures=c("Observed","Shannon","Simpson"))
# 按站点/时间比较
p_alpha <- plot_richness(ps_clean, x="site", measures=c("Observed","Shannon")) +
geom_boxplot() + theme_bw()
# Beta多样性
ord <- ordinate(ps_clean, method="NMDS", distance="bray")
plot_ordination(ps_clean, ord, color="site") + stat_ellipse()
# 物种占域模型（考虑检测不完美性）
library(unmarked)
# 准备检测历史数据（每个站点多次采样）
# y矩阵: sites × visits (1=detected, 0=not detected)
y_detection <- matrix(c(1,1,0, 0,1,1, 1,0,0, ...), ncol=3, byrow=TRUE)
# 单季节占域模型
umf <- unmarkedFrameOccu(y=y_detection,
siteCovs=site_covariates,
obsCovs=list(effort=effort_matrix))
fm <- occu(~effort ~habitat + elevation, data=umf)
summary(fm)
# 估计真实占域率（考虑检测概率）
predicted_occupancy <- predict(fm, type="state")
实战命令¶
#!/bin/bash
# === eDNA metabarcoding完整流程 ===
set -e
PRIMERS_F="GTCGGTAAAACTCGTGCCAGC"
PRIMERS_R="CATAGTGGGGTATCTAATCCCAGTTTG"
REF_DB="reference/fish_12S_db"
THREADS=16
# 1. 去引物
for f in raw/*_R1.fastq.gz; do
sample=$(basename $f _R1.fastq.gz)
cutadapt -g $PRIMERS_F -G $PRIMERS_R \
--discard-untrimmed --minimum-length 50 \
-o trimmed/${sample}_R1.fastq.gz -p trimmed/${sample}_R2.fastq.gz \
raw/${sample}_R1.fastq.gz raw/${sample}_R2.fastq.gz
done
# 2. DADA2 (在R中运行)
Rscript run_dada2.R
# 3. BLAST分类注释
blastn -query asv_sequences.fasta -db $REF_DB \
-out blast_results.txt -outfmt "6 qseqid sseqid pident length qcovs staxids" \
-max_target_seqs 10 -evalue 1e-10 -num_threads $THREADS
# 4. LCA分类(取最低共同祖先)
python3 lca_assignment.py blast_results.txt > taxonomy_lca.tsv
面试常问点¶
Q1: eDNA相比传统生物调查方法的优势和局限？ A: 优势：非侵入性（不伤害生物）、高灵敏度（能检测稀有种）、可覆盖大范围分类群、采样快速标准化。局限：不能确定活体（可能检测到死体DNA）、不能确定个体数量（半定量）、引物偏好导致部分物种漏检、参考数据库不完整、环境DNA降解导致时间分辨率低。
Q2: 如何控制eDNA研究中的假阳性？ A: 1）严格的阴性对照体系（采样空白、提取空白、PCR空白）；2）使用decontam统计去污染；3）设置最小读数阈值（如<5 reads视为噪声）；4）要求多次重复检测才确认（如3次采样至少2次检出）；5）实验室分区防止交叉污染。
Q3: OTU和ASV在eDNA分析中分别适用什么场景？ A: ASV有单碱基分辨率、可跨数据集比较、不丢失真实变异信息，是当前推荐。OTU (97%聚类)可能合并近缘种（如COI中97%一致性可能跨种）。但对于eDNA短片段（如teleo的60bp），即使ASV也可能无法区分所有物种。
Q4: eDNA能做定量分析吗？ A: eDNA读数与生物量的相关性有限。影响因素：PCR偏好（不同物种扩增效率不同）、DNA释放率差异（鱼类>>两栖类）、DNA降解速率、过滤效率。相对丰度可做趋势参考，但不能直接推断绝对生物量。qPCR配合标准曲线对单物种定量较好。
易错点¶
引物选择不当：通用引物可能对某些物种有严重扩增偏好（甚至不扩增），导致系统性漏检。需要用in silico PCR评估覆盖度。
tag jumping未处理：Illumina平台的index hopping/tag switching会导致reads被分配到错误样本（假阳性率可达0.1-2%）。需要使用独特双索引(UDI)和设置最小读数阈值。
参考数据库不完整：许多地方性物种（尤其是无脊椎动物）未被参考数据库覆盖，导致无法注释到种水平。
忽略检测概率：eDNA检测并非100%灵敏，未检出不等于物种不存在。应使用占域模型(occupancy model)估计真实存在概率。
片段长度与eDNA降解不匹配：环境中的eDNA多为短片段（<200bp），设计过长的扩增子会降低灵敏度。
多次PCR的reads直接合并：同一样本的多次PCR重复应该先独立处理再合并，直接合并可能放大PCR bias。
补充知识¶
eDNA标记基因选择指南¶
| 目标类群 | 推荐标记 | 引物 | 片段长度 | 分辨率 |
|---|---|---|---|---|
| 鱼类 | 12S rRNA | MiFish-U | ~170bp | 种级 |
| 哺乳动物 | 12S/16S | MiMammal | ~170bp | 种级 |
| 两栖类 | 12S/16S | AmphibiaDNA | ~100bp | 属/种 |
| 昆虫 | COI | mlCOIintF | ~313bp | 种级 |
| 植物 | trnL-P6 | gh/c | ~10-143bp | 科/属 |
| 真菌 | ITS | ITS1F/ITS2 | ~250bp | 种级 |