微生物菌株水平分析¶

# 菌株分析工具对比
# | 工具 | 原理 | 输入 | 优势 |
# |------|------|------|------|
# | StrainPhlAn | marker gene SNP | MetaPhlAn输出 | 快速，无需组装 |
# | inStrain | reads→genome mapping | BAM+genome | 定量精确，种群遗传 |
# | DESMAN | 共变异推断菌株 | contig variants | 混合菌株分解 |
# | StrainGE | reference-guided | FASTQ | 直接菌株检测 |

# === StrainPhlAn 4 流程 ===

# 前提：已运行MetaPhlAn4获得比对输出文件
# 注意：MetaPhlAn 4.2起参数名改为 --mapout（旧版用 --bowtie2out）
metaphlan sample.fq.gz --input_type fastq \
    --bowtie2out sample.bowtie2.bz2 \
    -o sample_profile.txt --nproc 16

# Step1: 从MetaPhlAn比对中提取样本标记基因
# 注意：StrainPhlAn 4.1起输出从.pkl改为.json格式
sample2markers.py -i sample.bowtie2.bz2 \
    -o consensus_markers/ \
    --nprocs 8

# Step2: 提取目标物种的参考标记基因
extract_markers.py -c s__Escherichia_coli \
    -o ref_markers/

# Step3: 运行StrainPhlAn（4.0用*.pkl，4.1+用*.json）
strainphlan -s consensus_markers/*.pkl \
    -m ref_markers/s__Escherichia_coli.fna \
    -r ref_genomes/*.fna \
    -o strainphlan_output/ \
    -n 8 \
    -c s__Escherichia_coli \
    --mutation_rates \
    --phylophlan_mode accurate

# 输出：
# - RAxML系统发育树
# - 多序列比对
# - 突变率统计

# StrainPhlAn结果解析
import pandas as pd
from ete3 import Tree

# 读取系统发育树
tree = Tree("strainphlan_output/RAxML_bestTree.s__Escherichia_coli.tre")

# 计算样本间遗传距离
# 距离近的样本=可能同一菌株传播
for leaf1 in tree.get_leaves():
    for leaf2 in tree.get_leaves():
        if leaf1.name < leaf2.name:
            dist = tree.get_distance(leaf1, leaf2)
            if dist < 0.001:  # 非常近
                print(f"Possible same strain: {leaf1.name} - {leaf2.name} (dist={dist:.6f})")

# === inStrain 流程 ===

# Step1: 将reads比对到参考基因组
bowtie2 --very-sensitive -p 16 \
    -x reference_genome \
    -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
    | samtools sort -@ 8 -o sample_vs_ref.sorted.bam

samtools index sample_vs_ref.sorted.bam

# Step2: 运行inStrain profile
inStrain profile sample_vs_ref.sorted.bam reference_genome.fasta \
    -o inStrain_output/ \
    -p 16 \
    --min_cov 5 \
    --min_freq 0.05 \
    --database_mode

# Step3: 多样本比较
inStrain compare \
    -i sample1_inStrain/ sample2_inStrain/ sample3_inStrain/ \
    -o compare_output/ \
    -p 16 \
    --database_mode

# 关键输出指标：
# - popANI: 种群平均核苷酸一致性（>99.999%=同一菌株）
# - conANI: 共识序列ANI
# - coverage: 基因组覆盖度
# - breadth: 基因组覆盖宽度
# - nucl_diversity: 核苷酸多样性(π)

# inStrain结果分析
import inStrain
import inStrain.SNVprofile
import pandas as pd

# 读取profile结果
IS = inStrain.SNVprofile.SNVprofile("inStrain_output/")

# 获取基因组级别统计（读取输出目录下的TSV文件）
genome_info = pd.read_csv("inStrain_output/output/genome_info.tsv", sep="\t")
print(genome_info[['genome', 'coverage', 'breadth', 'nucl_diversity']])

# 读取compare结果
compare_df = pd.read_csv("compare_output/comparisonsTable.tsv", sep="\t")

# 判断同一菌株
# popANI > 99.999% 且 共享位点 > 覆盖的50%
same_strain = compare_df[
    (compare_df['popANI'] > 0.99999) & 
    (compare_df['percent_genome_compared'] > 0.5)
]
print("Likely same strain pairs:")
print(same_strain[['name1', 'name2', 'popANI']])

# === DESMAN 流程 ===

# 前提：已有contig binning结果，选定目标物种的bin

# Step1: 从BAM中提取变异位点频率
python Variant_Filter.py contig_coverage.csv \
    -o filtered_variants/ \
    -m 0.01 -f 0.05

# Step2: 运行DESMAN（尝试不同菌株数G）
for G in 2 3 4 5; do
    desman filtered_variants/freq_var.csv \
        -e filtered_variants/freq_df.csv \
        -o desman_G${G}/ \
        -g ${G} \
        -s 100 \
        -i 100 \
        --seed 42
done

# Step3: 选择最优菌株数（基于ELBO/posterior deviance）
python resolvenhap.py desman_G2/ desman_G3/ desman_G4/ desman_G5/

# 输出：
# - 各菌株的基因型（在每个变异位点的等位基因）
# - 各样本中各菌株的相对频率
# - 模型拟合度（用于选择G）

# === 母婴菌株传递分析 ===

# 使用inStrain compare的结果
compare_results <- read.csv("compare_output/comparisonsTable.tsv", sep = "\t")

# 筛选母婴配对的比较结果
mother_infant_pairs <- compare_results %>%
  filter(grepl("Mother", name1) & grepl("Infant", name2))

# 判断传递：popANI > 99.999%
transmitted_strains <- mother_infant_pairs %>%
  filter(popANI > 0.99999 & percent_genome_compared > 0.5)

cat(sprintf("Transmitted strains: %d out of %d species compared\n",
            nrow(transmitted_strains), nrow(mother_infant_pairs)))

# 传递率统计（按物种）
transmission_rate <- mother_infant_pairs %>%
  group_by(genome) %>%
  summarise(
    n_pairs = n(),
    n_transmitted = sum(popANI > 0.99999, na.rm = TRUE),
    rate = n_transmitted / n_pairs
  ) %>%
  arrange(desc(rate))

# 可视化菌株传递网络
library(igraph)
edges <- transmitted_strains[, c("name1", "name2")]
g <- graph_from_data_frame(edges, directed = TRUE)
plot(g, vertex.size = 8, edge.arrow.size = 0.5)

# === 菌株功能差异 ===

# 使用inStrain的基因水平分析
# 新版inStrain推荐在profile阶段直接用-g指定基因文件，gene_profile已废弃
inStrain profile sample_vs_ref.sorted.bam reference.fasta \
    -g reference_genes.fna \
    -o gene_profile_output/ \
    -p 16

# 基因存在/缺失变异（accessory genome）
# 通过覆盖度判断基因是否存在
# coverage > 0.5 * genome_mean = present
# coverage < 0.1 * genome_mean = absent

# 菌株水平的基因含量差异
import pandas as pd

# 读取基因覆盖度
gene_coverage = pd.read_csv("gene_profile_output/gene_info.tsv", sep="\t")

# 标记基因存在/缺失
genome_mean_cov = gene_coverage['coverage'].mean()
gene_coverage['present'] = gene_coverage['coverage'] > 0.5 * genome_mean_cov

# 统计可变基因（在部分样本存在）
gene_prevalence = gene_coverage.groupby('gene')['present'].mean()
accessory_genes = gene_prevalence[(gene_prevalence > 0.1) & (gene_prevalence < 0.9)]
print(f"Accessory genes: {len(accessory_genes)}")

# 耐药基因/毒力因子的菌株特异性
# 检查哪些样本的菌株携带特定功能基因

# StrainPhlAn
sample2markers.py -i sample.bowtie2.bz2 -o markers/
strainphlan -s markers/*.pkl -m ref_markers.fna -o output/ -c species_name

# inStrain
bowtie2 --very-sensitive -x ref -1 R1.fq.gz -2 R2.fq.gz | samtools sort -o aligned.bam
inStrain profile aligned.bam ref.fa -o profile_out/ -p 16
inStrain compare -i profile1/ profile2/ -o compare_out/