20260910 113033 宏基因组微生物分类鉴定全流程实现详解

fastp -i R1.fastq.gz -I R2.fastq.gz \
    -o clean_R1.fastq.gz -O clean_R2.fastq.gz \
    --qualified_quality_phred 20 --length_required 50 \
    --detect_adapter_for_pe --thread 8

bowtie2 -x human_genome_index -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \
    --un-conc-gz host_removed_R%.fastq.gz --very-sensitive

kraken2 --db "$KRAKEN_DB" \
    --paired --gzip-compressed --threads 8 \
    --confidence 0.1 \            # 从默认0提高到0.1-0.4，抑制单k-mer假阳性
    --minimum-hit-groups 2 \      # 要求>=2个独立hit区域（临床建议3）
    --report out.kreport \
    --output out.kraken \
    R1.fq.gz R2.fq.gz

kraken2 --db "$KRAKEN_DB" --paired --confidence 0.1 \
    --report-minimizer-data \     # 插入2列：总 minimizer + 独特 minimizer
    --report out.krep --output /dev/null \
    R1.fq.gz R2.fq.gz
# 高reads但低独特minimizer的物种 = 假阳性（单一保守区域反复命中）

kraken2-build --download-taxonomy --db custom_db
for lib in bacteria archaea viral human UniVec_Core; do
    kraken2-build --download-library "$lib" --db custom_db
done
kraken2-build --build --db custom_db --threads 16
kraken2-build --clean --db custom_db   # 删除 library/ + taxonomy/ 减小体积
bracken-build -d custom_db -t 16 -k 35 -l 150   # -k 必须等于 Kraken2 的 k(35)；-l = 读长

bracken -d "$KRAKEN_DB" -i out.kreport -o out.bracken -w out.bracken.kreport \
    -r 150 \   # 必须匹配 built databaseRLENmers.kmer_distrib 且等于真实读长
    -l S -t 10 # species level；-t 是重分配下限，不是置信度

# paired-end 作为 ONE 逗号分隔参数（MetaPhlAn 按单端处理，不用配对信息）
metaphlan reads_R1.fastq.gz,reads_R2.fastq.gz \
    --input_type fastq \
    --index mpa_vJun23_CHOCOPhlAnSGB_202403 \   # 钉住 index，数据库版本是批次变量
    --nproc 8 \
    --mapout sample.map.bz2 \                    # 缓存 read->marker 映射（4.2 前是 --bowtie2out）
    --output_file profile.txt

metaphlan sample.map.bz2 --input_type mapout \   # 前4.2是 bowtie2out
    --tax_lev s \            # k,p,c,o,f,g,s,t
    --stat_q 0.2 \          # 截断均值：去掉上下20%，平均中间60%
    --output_file profile_species.txt

# 预QC：接头/质量trim + 去宿主（例如 KneadData）。宿主读段膨胀 UNMAPPED 且浪费 DIAMOND 时间
cat sample_R1.fq.gz sample_R2.fq.gz > sample.fq.gz   # 配对读段没有原生配对，拼接成单文件
humann --input sample.fq.gz --output out/ \
    --taxonomic-profile sample_metaphlan.tsv \   # 复用 MetaPhlAn profile；不要 --remove-temp-output
    --threads 8   # 默认：prescreen 0.01, translated-id 80, gap-fill on, minpath on

# 归一化：跨样本统计前先 PER SAMPLE 归一化（RPK 是深度依赖的）
humann_renorm_table -i out/sample_genefamilies.tsv -o out/sample_cpm.tsv -u cpm
humann_join_tables -i out -o merged_pathabundance.tsv --file_name pathabundance
humann_regroup_table -i merged_pathabundance.tsv -g uniref90_ko -o merged_ko.tsv   # 会加 UNGROUPED 行
humann_split_stratified_table -i merged_pathabundance.tsv -o .   # 统计跑在 UNSTRATIFIED 文件上

import pandas as pd

df = pd.read_csv('merged_pathabundance_unstratified.tsv', sep='\t', index_col=0)
meta = pd.read_csv('metadata.tsv', sep='\t', index_col=0)
unmapped = df.loc[['UNMAPPED', 'UNINTEGRATED']]     # 分母，检查但不要丢弃
g1 = meta.index[meta['condition'] == 'healthy']
g2 = meta.index[meta['condition'] == 'disease']
unmapped_shift = unmapped[g2].mean(axis=1) - unmapped[g1].mean(axis=1)
print(unmapped_shift.round(1))
# 然后交给 MaAsLin2/ANCOM-BC（对组分与零膨胀建模），而不是裸 t-test/Mann-Whitney

qiime feature-table summarize --i-table table.qza --o-visualization table.qzv  # 深度信息在这
qiime diversity alpha-rarefaction \
    --i-table table.qza --i-phylogeny rooted-tree.qza \
    --p-max-depth 20000 \
    --m-metadata-file metadata.tsv --o-visualization alpha-rarefaction.qza
qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza --i-table table.qza \
    --p-sampling-depth 10000 \    # 在 observed-features 平台上；会静默丢弃低深度样本
    --m-metadata-file metadata.tsv --output-dir core-metrics-results

qiime fragment-insertion sepp \
    --i-representative-sequences rep-seqs.qda \
    --i-reference-database sepp-refs-gg-13-8.qda \
    --p-threads 4 \
    --o-tree insertion-tree.qda --o-placements insertion-placements.qda
qiime fragment-insertion filter-features \
    --i-table table.qda --i-tree insertion-tree.qda \
    --o-filtered-table table-sepp.qda --o-removed-table removed-table.qda

library(phyloseq); library(vegan)
ps_rare <- rarefy_even_depth(ps, sample.size = chosen_depth, rngseed = 42, replace = FALSE)
alpha <- estimate_richness(ps_rare, measures = c('Observed', 'Shannon', 'InvSimpson'))
alpha$Group <- sample_data(ps_rare)$Group
alpha$Shannon_eff <- exp(alpha$Shannon)   # 有效种数（Hill q=1），跨工具可比
kruskal.test(Shannon ~ Group, data = alpha)

wu  <- UniFrac(ps_rare, weighted = TRUE)    # 优势谱系视角
uwu <- UniFrac(ps_rare, weighted = FALSE)   # 稀有谱系 + 拓扑视角
gu  <- as.dist(GUniFrac::GUniFrac(t(as(otu_table(ps_rare), 'matrix')), phy_tree(ps_rare), alpha = 0.5)$unifracs[, , 'd_0.5'])

meta <- data.frame(sample_data(ps_rare))
adonis2(wu ~ Group, data = meta, permutations = 999)   # 显著 = 质心位置差异
permutest(betadisper(wu, meta$Group))                  # 强制：是离散度差异还是位置差异？

krepdeport2krona.py -r kraken.kreport -o krona.txt && ktImportText krona.txt -o krona.html

alpha <- estimate_richness(ps_counts, measures = c('Shannon', 'InvSimpson'))
hill_q1 <- exp(alpha$Shannon)   # 有效种数 q=1
hill_q2 <- alpha$InvSimpson     # q=2 优势度加权

library(vegan)
dist_bc <- vegdist(otu_matrix, method = 'bray')
pm <- adonis2(dist_bc ~ group, permutations = 999, by = 'terms')
bd <- betadisper(dist_bc, group); permutest(bd, permutations = 999)  # 必须配对

import pandas as pd
from skbio.stats.composition import clr, multi_replace
from sklearn.decomposition import PCA

ab = pd.read_csv('merged_abundance.txt', sep='\t', index_col=0)
ab = ab[ab.index.str.contains(r'\|s__') & ~ab.index.str.contains(r'\|t__')]  # 物种行
proportions = (ab.T.values / ab.T.values.sum(axis=1, keepdims=True))
clr_mat = clr(multi_replace(proportions))      # 零替换 + CLR（不是 StandardScaler）
pca = PCA(n_components=2).fit(clr_mat)
coords = pca.transform(clr_mat)

原始 shotgun reads（R1+R2）
  │
  ├─ Step 1: fastp 质控（去接头、去低质量）→ Bowtie2 宿主去除
  │
  ├─ Step 2: Kraken2 分类（--confidence 0.1 --minimum-hit-groups 2）
  │          → Bracken 丰度重估计（-r 读长匹配）
  │          → 假阳性控制（--report-minimizer-data / KrakenUniq）
  │
  ├─ Step 3: MetaPhlAn 4 标记基因剖面（细胞分数，uSGB 量化未知 taxa）
  │          （注意：与 Kraken 分数不可合并）
  │
  ├─ Step 4: HUMAnN 3 功能谱（MetaPhl 预筛 → Bowtie2 泛基因组 → DIAMOND）
  │          → 保留 UNMAPPED/UNINTEGRATED → 归一化到 CPM
  │
  ├─ Step 5: 多样性（α: Hill q=0/1/2；β: adonis2 + betadisper 配对）
  │
  └─ Step 6: 可视化（Krona / 堆叠柱 / CLR-PCA / 差异丰度共识）