373_食品微生物组分析¶

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

# 步骤1：导入食品样本的双端测序数据
qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path food_manifest.tsv \   # manifest文件：样本名-文件路径对应表
  --output-path food_sequences.qza \
  --input-format PairedEndFastqManifestPhred33V2

# 步骤2：DADA2 去噪（食品样本建议截短更短，因为基质干扰大）
qiime dada2 denoise-paired \
  --i-demultiplexed-seqs food_sequences.qza \
  --p-trim-left-f 20 \             # 去掉前20bp引物区
  --p-trim-left-r 20 \
  --p-trunc-len-f 220 \            # 正向截短到220bp（质量下降位点）
  --p-trunc-len-r 180 \            # 反向截短到180bp
  --p-n-threads 8 \                # 8线程并行
  --o-table food_table.qza \       # 输出：ASV丰度表
  --o-representative-sequences food_rep_seqs.qza \
  --o-denoising-stats food_stats.qza

# 步骤3：分类注释（使用 SILVA 数据库）
qiime feature-classifier classify-sklearn \
  --i-classifier silva-138-99-515-806-nb-classifier.qza \  # SILVA 138 分类器
  --i-reads food_rep_seqs.qza \
  --o-classification food_taxonomy.qza \
  --p-n-jobs 8

# 步骤4：导出丰度表供下游分析
qiime tools export \
  --input-path food_table.qza \
  --output-path food_exported/
biom convert -i food_exported/feature-table.biom \
  -o food_exported/feature-table.tsv --to-tsv   # 转换为可读tsv格式

import subprocess
import os

# 用 Kraken2 快速鉴定致病菌
def run_kraken2_pathogen_detection(input_r1, input_r2, output_dir, sample_name):
    """
    快速检测食品样本中的致病菌
    input_r1, input_r2: 双端fastq文件路径
    output_dir: 输出目录
    sample_name: 样本名称
    """
    os.makedirs(output_dir, exist_ok=True)  # 创建输出目录

    # 运行 Kraken2 比对到标准数据库
    cmd = [
        "kraken2",
        "--db", "/data/kraken2_standard_db",  # Kraken2标准数据库路径
        "--paired",                            # 双端模式
        "--threads", "8",
        "--output", f"{output_dir}/{sample_name}_kraken.txt",
        "--report", f"{output_dir}/{sample_name}_report.txt",
        "--classified-out", f"{output_dir}/{sample_name}_classified#.fastq",
        input_r1, input_r2
    ]
    subprocess.run(cmd, check=True)

    # 解析报告，提取致病菌
    pathogen_list = [
        "Listeria monocytogenes",
        "Salmonella enterica",
        "Escherichia coli",
        "Staphylococcus aureus"
    ]

    report_file = f"{output_dir}/{sample_name}_report.txt"
    with open(report_file) as f:
        for line in f:
            parts = line.strip().split("\t")
            if len(parts) >= 6:
                percent = float(parts[0])   # 比对百分比
                taxon_name = parts[5].strip()   # 物种名
                # 检查是否为致病菌且占比超过0.01%
                if any(p in taxon_name for p in pathogen_list) and percent > 0.01:
                    print(f"[警告] 检测到致病菌: {taxon_name} ({percent:.3f}%)")

# 示例调用
run_kraken2_pathogen_detection(
    "sample1_R1.fastq.gz",
    "sample1_R2.fastq.gz",
    "pathogen_results",
    "food_sample1"
)