606 真菌ITS测序分析¶

一句话概述： ITS（内转录间隔区）是真菌分类的"通用条码"，通过ITS1/ITS2扩增子测序+QIIME2/DADA2分析流程，可以揭示样本中真菌群落的组成和多样性。

核心知识点速查表¶

知识点	白话解释
ITS	Internal Transcribed Spacer，核糖体RNA基因之间的间隔区，变异大，适合真菌分类
ITS1	18S和5.8S之间的间隔区，长约200-600bp
ITS2	5.8S和28S之间的间隔区，长约200-400bp
UNITE	真菌ITS分类的"金标准"参考数据库（相当于细菌的SILVA）
ITSxpress	从测序数据中精确提取ITS1或ITS2区域的工具
DADA2	去噪算法，将原始reads校正为精确的ASV（扩增子序列变体）
ASV	Amplicon Sequence Variant，精确到单碱基分辨率的序列单元
OTU	Operational Taxonomic Unit，按97%相似性聚类的序列单元（旧方法）
真菌组(mycobiome)	某个环境中所有真菌的集合
cutadapt	去除引物序列的工具（ITS分析中必用）

一、ITS vs 16S：为什么真菌不用16S？¶

白话解释：

细菌用16S rRNA基因做分类，因为16S在不同细菌之间"刚好够不同"。但真菌的16S/18S rRNA基因太保守了，不同真菌之间几乎一样，分不开。

ITS区域就像真菌的"指纹"：

rRNA基因组结构：
18S ─── ITS1 ─── 5.8S ─── ITS2 ─── 28S
(保守)  (高变)   (保守)   (高变)   (保守)
         ↑                  ↑
   用这两个区域做真菌分类

对比项	16S rRNA（细菌）	ITS（真菌）
靶标	16S核糖体RNA基因	核糖体间隔区
长度	约1542bp，V3V4约450bp	ITS1约200-600bp，ITS2约200-400bp
长度变异	小（不同细菌长度接近）	大（不同真菌ITS长度差异很大！）
参考数据库	SILVA / Greengenes2	UNITE
DADA2处理	可以设固定截断长度	不能设固定截断长度（长度变异大）

ITS1 vs ITS2 怎么选？ - ITS1：变异更大，分辨率更高，但引物覆盖度稍差 - ITS2：长度更均一，DADA2处理更友好，2024年多数研究推荐 - 建议：没有特殊需求就选ITS2

二、分析流程详解¶

2.1 数据导入与质控¶

# ========== 步骤1：导入数据到QIIME2 ==========

# 激活QIIME2环境
conda activate qiime2-amplicon-2024.10  # 使用最新版QIIME2

# 导入双端数据（Casava格式）
qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \  # 数据类型：双端带质量
    --input-path raw_data/ \                                # 输入：原始数据目录
    --input-format CasavaOneEightSingleLanePerSampleDirFmt \  # Casava格式
    --output-path demux.qza                                 # 输出：QIIME2格式

# 查看数据质量
qiime demux summarize \
    --i-data demux.qza \            # 输入：导入的数据
    --o-visualization demux.qzv     # 输出：可视化文件（用QIIME2 View查看）
# 在 https://view.qiime2.org/ 上传qzv文件查看质量分布

2.2 去除引物序列（关键步骤）¶

# ========== 步骤2：去除引物序列 ==========
# ITS引物可能部分进入测序读段，必须去除
# 否则引物序列会被当作生物序列参与下游分析

# 常用ITS引物：
# ITS1F: CTTGGTCATTTAGAGGAAGTAA
# ITS2:  GCTGCGTTCTTCATCGATGC
# ITS3:  GCATCGATGAAGAACGCAGC
# ITS4:  TCCTCCGCTTATTGATATGC

# 用cutadapt去引物
qiime cutadapt trim-paired \
    --i-demultiplexed-sequences demux.qza \          # 输入
    --p-front-f CTTGGTCATTTAGAGGAAGTAA \             # 正向引物（ITS1F）
    --p-front-r GCTGCGTTCTTCATCGATGC \               # 反向引物（ITS2）
    --p-discard-untrimmed \                           # 丢弃没有引物的reads
    --p-match-read-wildcards \                        # 允许引物中的简并碱基
    --p-match-adapter-wildcards \                     # 允许接头中的简并碱基
    --p-minimum-length 50 \                           # 去引物后最短50bp
    --o-trimmed-sequences trimmed.qza \               # 输出：去引物后的数据
    --verbose                                         # 显示详细信息

# 检查去引物效果
qiime demux summarize \
    --i-data trimmed.qza \             # 输入：去引物后的数据
    --o-visualization trimmed.qzv      # 输出：可视化

2.3 ITSxpress提取ITS区域¶

# ========== 步骤3：用ITSxpress精确提取ITS区域 ==========
# 为什么要这一步？
# ITS两侧有保守的rRNA序列（18S和5.8S的片段）
# 这些保守序列会干扰ASV的推断，必须去掉
# ITSxpress v2（2024更新）保留质量分数，适合后续DADA2

# 安装ITSxpress QIIME2插件
pip install q2-itsxpress  # 安装插件

# 提取ITS2区域（示例）
qiime itsxpress trim-pair-output-unmerged \
    --i-per-sample-sequences trimmed.qza \       # 输入：去引物后的数据
    --p-region ITS2 \                            # 提取ITS2区域
    --p-taxa F \                                 # 目标类群：F=真菌
    --p-threads 16 \                             # 16个线程
    --o-trimmed its2_extracted.qza               # 输出：提取ITS2后的数据

# 查看提取效果
qiime demux summarize \
    --i-data its2_extracted.qza \
    --o-visualization its2_extracted.qzv

2.4 DADA2去噪¶

# ========== 步骤4：DADA2去噪 ==========
# ITS分析的DADA2参数与16S不同！
# 核心区别：ITS长度变异大，不能设固定的trunc-len

# 方法1：不截断（推荐用于ITS）
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs its2_extracted.qza \  # 输入：提取ITS后的数据
    --p-trunc-len-f 0 \                          # 正向不截断（0=不截断）
    --p-trunc-len-r 0 \                          # 反向不截断
    --p-max-ee-f 2.0 \                           # 正向最大期望错误数
    --p-max-ee-r 2.0 \                           # 反向最大期望错误数
    --p-trunc-q 2 \                              # 遇到质量值2就截断
    --p-min-fold-parent-over-abundance 1.0 \     # 嵌合体检测参数
    --p-n-threads 16 \                           # 16个线程
    --o-table table.qza \                        # 输出：特征表（ASV丰度矩阵）
    --o-representative-sequences rep-seqs.qza \  # 输出：代表序列
    --o-denoising-stats stats.qza                # 输出：去噪统计

# 查看去噪统计
qiime metadata tabulate \
    --m-input-file stats.qza \        # 输入：去噪统计
    --o-visualization stats.qzv       # 输出：可视化

# 查看特征表摘要
qiime feature-table summarize \
    --i-table table.qza \             # 输入：特征表
    --o-visualization table.qzv \     # 输出：可视化
    --m-sample-metadata-file metadata.tsv  # 样本元数据

2.5 物种分类注释¶

# ========== 步骤5：物种分类（基于UNITE数据库） ==========

# 下载UNITE数据库（2025年2月最新版 v10）
# 从 https://unite.ut.ee/repository.php 下载
# 选择 "QIIME release" → "All eukaryotes" → "dynamic"

# 训练分类器（首次使用需要，比较耗时）
qiime feature-classifier fit-classifier-naive-bayes \
    --i-reference-reads unite_refs.qza \              # UNITE参考序列
    --i-reference-taxonomy unite_taxonomy.qza \       # UNITE分类信息
    --o-classifier unite_classifier.qza               # 输出：训练好的分类器

# 分类注释
qiime feature-classifier classify-sklearn \
    --i-classifier unite_classifier.qza \      # 输入：训练好的分类器
    --i-reads rep-seqs.qza \                   # 输入：ASV代表序列
    --p-n-jobs 16 \                            # 16个线程
    --p-confidence 0.7 \                       # 置信度阈值（0.7是默认值）
    --o-classification taxonomy.qza            # 输出：分类结果

# 查看分类结果
qiime metadata tabulate \
    --m-input-file taxonomy.qza \
    --o-visualization taxonomy.qzv

# 生成分类柱状图
qiime taxa barplot \
    --i-table table.qza \                # 特征表
    --i-taxonomy taxonomy.qza \          # 分类结果
    --m-sample-metadata-file metadata.tsv \  # 样本元数据
    --o-visualization taxa-bar-plots.qzv     # 输出：分类柱状图

# 补充：用BLASTn做分类（对某些真菌更准确）
qiime feature-classifier classify-consensus-blast \
    --i-query rep-seqs.qza \              # 查询序列
    --i-reference-reads unite_refs.qza \  # 参考序列
    --i-reference-taxonomy unite_taxonomy.qza \  # 参考分类
    --p-perc-identity 0.97 \              # 最低一致性97%
    --o-classification taxonomy_blast.qza  # 输出分类结果

2.6 多样性分析¶

# ========== 步骤6：Alpha和Beta多样性分析 ==========

# 构建系统发育树（用于UniFrac等需要进化信息的指标）
qiime phylogeny align-to-tree-mafft-fasttree \
    --i-sequences rep-seqs.qza \          # 输入：代表序列
    --o-alignment aligned-rep-seqs.qza \  # 输出：比对结果
    --o-masked-alignment masked-aligned.qza \  # 输出：mask后的比对
    --o-tree unrooted-tree.qza \          # 输出：无根树
    --o-rooted-tree rooted-tree.qza       # 输出：有根树

# 核心多样性分析（一步生成所有常用指标）
qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza \        # 系统发育树
    --i-table table.qza \                  # 特征表
    --p-sampling-depth 5000 \              # 抽平深度（看table.qzv确定）
    --m-metadata-file metadata.tsv \       # 样本元数据
    --output-dir diversity_results/        # 输出目录

# Alpha多样性组间比较
qiime diversity alpha-group-significance \
    --i-alpha-diversity diversity_results/shannon_vector.qza \  # Shannon指数
    --m-metadata-file metadata.tsv \      # 元数据
    --o-visualization shannon_significance.qzv  # 输出可视化

# Beta多样性组间比较（PERMANOVA）
qiime diversity beta-group-significance \
    --i-distance-matrix diversity_results/bray_curtis_distance_matrix.qza \  # Bray-Curtis距离
    --m-metadata-file metadata.tsv \       # 元数据
    --m-metadata-column group \            # 分组列
    --p-pairwise \                         # 两两比较
    --o-visualization bray_curtis_significance.qzv  # 输出

2.7 差异丰度分析¶

# ========== 步骤7：差异丰度分析（R语言） ==========
library(phyloseq)     # 微生物组分析框架
library(DESeq2)       # 差异分析
library(ggplot2)      # 可视化

# 从QIIME2导入数据到R
library(qiime2R)  # QIIME2-R桥接包

# 读取QIIME2输出
physeq <- qza_to_phyloseq(
  features = "table.qza",       # 特征表
  tree = "rooted-tree.qza",     # 系统发育树
  taxonomy = "taxonomy.qza",    # 分类注释
  metadata = "metadata.tsv"     # 样本元数据
)

# 过滤低丰度ASV（至少在5%的样本中出现）
physeq_filt <- filter_taxa(physeq,               # 输入phyloseq对象
                           function(x) sum(x > 0) > (0.05 * length(x)),  # 过滤条件
                           TRUE)                  # 保留通过的

# 用DESeq2做差异分析
dds <- phyloseq_to_deseq2(physeq_filt,    # 转换为DESeq2对象
                           ~ group)         # 按group分组

dds <- DESeq(dds, test = "Wald",           # Wald检验
             fitType = "parametric")        # 参数拟合

# 提取结果
res <- results(dds, contrast = c("group", "Disease", "Healthy"))  # 疾病vs健康
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)     # 显著差异

# 查看显著差异的真菌
cat("显著差异ASV数:", nrow(res_sig), "\n")  # 打印数量

# 将结果与分类信息关联
tax_table_df <- as.data.frame(tax_table(physeq_filt))  # 提取分类表
sig_taxa <- tax_table_df[rownames(res_sig), ]           # 显著差异ASV的分类信息
print(sig_taxa[, c("Phylum", "Class", "Order", "Family", "Genus")])  # 打印分类

三、ITS分析的特殊注意事项¶

========== ITS分析 vs 16S分析的关键区别 ==========

1. DADA2截断长度：
   16S：设固定截断长度（如trunc-len-f=240, trunc-len-r=200）
   ITS：设为0（不截断），因为ITS长度变异太大

2. 引物去除：
   16S：可以用DADA2的trimLeft参数
   ITS：必须用cutadapt，因为引物可能包含简并碱基

3. ITS区域提取：
   16S：不需要
   ITS：需要用ITSxpress去除保守的rRNA侧翼序列

4. 参考数据库：
   16S：SILVA / Greengenes2
   ITS：UNITE（唯一选择）

5. 系统发育树：
   16S：可以构建（序列可以可靠比对）
   ITS：比对质量差（长度变异大），UniFrac可能不可靠
   建议ITS用Bray-Curtis而非UniFrac做Beta多样性

6. 分类器选择：
   16S：Naive Bayes效果好
   ITS：有时BLASTn更准确（特别是低丰度未知真菌）

常见报错与解决¶

报错	原因	解决方案
`DADA2: zero reads passed filter`	截断长度设太长，ITS短序列全被丢弃	设trunc-len为0（不截断）
`ITSxpress: No ITS region found`	引物没去干净，或序列方向反了	先用cutadapt去引物，检查序列方向
`cutadapt: 0 reads with adapters`	引物序列写错了	检查引物序列是否与测序文库匹配
`UNITE分类器训练失败`	内存不足（UNITE数据库很大）	至少16G内存，或用预训练的分类器
`大部分ASV分类为Unidentified`	UNITE数据库不完整或置信度太高	降低confidence（如0.6），或用BLASTn
`ITS2长度分布异常（偏短）`	未去除引物就提取ITS	先cutadapt去引物，再ITSxpress
`Alpha多样性组间无差异`	抽平深度不合适	尝试不同的抽平深度，检查样本量是否足够

速查表¶

========== ITS分析标准流程 ==========
原始数据 → 导入QIIME2 → cutadapt去引物 → ITSxpress提取ITS
    → DADA2去噪(trunc-len=0) → UNITE分类注释
    → 多样性分析 → 差异丰度分析 → 可视化

========== 核心工具 ==========
框架平台 ：QIIME2 / mothur
去噪     ：DADA2（推荐）/ Deblur
引物去除 ：cutadapt
ITS提取  ：ITSxpress v2（2024更新，保留质量分数）
分类注释 ：Naive Bayes + UNITE / BLASTn + UNITE
多样性   ：QIIME2 diversity plugin
差异分析 ：DESeq2 / ANCOM-BC / ALDEx2
可视化   ：R (phyloseq + ggplot2) / QIIME2 View

========== 常用引物对 ==========
ITS1区域：ITS1F + ITS2    → ~250-600bp
ITS2区域：ITS3 + ITS4     → ~200-400bp
ITS全长 ：ITS1F + ITS4    → ~400-900bp（不推荐Illumina短读段）

========== UNITE数据库版本 ==========
最新版：v10（2025年2月发布）
选择：All eukaryotes → dynamic → developer format（QIIME2用）
注意：每年更新，建议使用最新版

========== DADA2参数对比 ==========
                  16S分析         ITS分析
trunc-len-f    ：240             0（不截断）
trunc-len-r    ：200             0（不截断）
max-ee         ：2.0             2.0
trunc-q        ：2               2
引物去除       ：trimLeft        cutadapt（必须）

面试高频问题¶

Q1：为什么真菌不用16S rRNA做分类？¶

答： 16S rRNA在真菌中（实际是18S rRNA）太保守了，不同真菌属之间18S的序列差异很小，分辨率不够。而ITS区域（内转录间隔区）位于18S和28S之间，进化速度快，不同物种之间差异大，非常适合做真菌分类。ITS被国际条码生命项目（iBOL）确定为真菌DNA条码标准区域。

Q2：ITS1和ITS2怎么选？有什么区别？¶

答： - ITS1：位于18S和5.8S之间，变异更大，分辨率更高，但不同真菌间长度差异很大（200-600bp），DADA2处理时可能出问题 - ITS2：位于5.8S和28S之间，长度更均一（200-400bp），更适合DADA2处理 - 2024年的研究表明ITS2在DADA2处理时产生的假阳性更少 - 建议：肠道真菌研究多用ITS2，土壤真菌研究有时用ITS1

Q3：ITS分析中DADA2的参数设置与16S有何不同？¶

答：最关键的区别是截断长度（trunc-len）设为0。

原因：16S的V3V4区域长度基本一致（约450bp），可以设固定截断长度去掉3'端低质量碱基。但ITS的长度变异很大（200-600bp），如果设固定截断长度，短ITS序列会被全部丢弃，导致某些真菌类群完全消失。

设trunc-len=0表示不截断，让DADA2根据质量分数（trunc-q=2）自动决定截断位置。

Q4：UNITE数据库是什么？怎么选择版本？¶

答： UNITE是真菌ITS分类的唯一主流参考数据库，由爱沙尼亚塔尔图大学维护： - dynamic版本：按物种假说（Species Hypothesis）动态聚类，推荐使用 - static版本：按97%相似性固定聚类 - All eukaryotes版本：包含所有真核生物（推荐，可以检测到非真菌的真核生物） - developer format：适合QIIME2使用的格式

最新版本：v10（2025年2月发布），建议始终使用最新版。

Q5：ITS分析结果中大量"Unidentified"怎么办？¶

答：真菌ITS分析中高比例的"未鉴定"是正常的，原因： 1. 真菌多样性极高，UNITE数据库覆盖不完整 2. 很多环境真菌还没有被培养和鉴定

解决方法： 1. 降低分类器的置信度阈值（从0.7降到0.6或0.5） 2. 用BLASTn代替Naive Bayes做分类 3. 对"Unidentified"的ASV做BLAST手动查询 4. 接受一定比例的Unidentified，在结果中报告为"未鉴定真菌"