622_微生物溯源与传播链¶

# === 基因组比较和分型工具 ===
conda install -c bioconda snippy      # SNP分析核心工具
conda install -c bioconda gubbins     # 重组检测
conda install -c bioconda snp-dists   # SNP距离矩阵
conda install -c bioconda mlst        # MLST分型
conda install -c bioconda chewbbaca   # cgMLST/wgMLST
conda install -c bioconda iqtree      # 建树

# === 传播链推断工具 ===
pip install transphylo  # 传播树推断（Python接口）

# === 用Snippy找SNP ===
# 第1步：每个样本与参考基因组比较
for sample in samples/*.fastq.gz; do
    name=$(basename "$sample" _R1.fastq.gz)
    snippy \
        --cpus 8 \                         # 线程数
        --outdir snippy_${name} \          # 输出目录
        --ref reference.fasta \            # 参考基因组
        --R1 ${sample} \                   # 正向reads
        --R2 ${sample/_R1/_R2}             # 反向reads
done

# 第2步：生成核心SNP比对
snippy-core \
    --ref reference.fasta \                # 参考基因组
    snippy_*/                              # 所有Snippy结果目录
# 输出：core.aln（核心SNP比对）, core.full.aln（全长度比对）

# 第3步：只提取变异位点
snippy-clean_full_aln core.full.aln > clean.full.aln  # 清理比对
snp-sites -c clean.full.aln > core_snps.aln  # 只保留SNP位点

# Gubbins去除重组产生的SNP（对细菌很重要！）
run_gubbins.py \
    clean.full.aln \           # 输入全长度比对
    --threads 8 \              # 线程数
    --tree-builder iqtree \    # 用IQ-TREE建树
    --prefix gubbins_output    # 输出前缀

# 输出：
# gubbins_output.filtered_polymorphic_sites.fasta — 去除重组后的SNP
# gubbins_output.final_tree.tre — 最终树
# gubbins_output.recombination_predictions.gff — 重组区域

# 计算成对SNP距离
snp-dists core_snps.aln > snp_distance_matrix.tsv

# 查看距离矩阵
cat snp_distance_matrix.tsv

# R语言可视化SNP距离热图
library(pheatmap)  # 加载pheatmap

# 读取SNP距离矩阵
dist <- read.table("snp_distance_matrix.tsv", 
                    header = TRUE, row.names = 1, sep = "\t")

# 画热图
pheatmap(
    as.matrix(dist),
    display_numbers = TRUE,     # 显示数值
    number_format = "%.0f",     # 整数格式
    fontsize_number = 8,        # 数字大小
    color = colorRampPalette(c("white", "yellow", "red"))(100),
    main = "Pairwise SNP Distances"  # 标题
)

# 用IQ-TREE建树
iqtree2 \
    -s gubbins_output.filtered_polymorphic_sites.fasta \
    -m GTR+G4 \               # GTR模型（或用MFP自动选）
    -bb 1000 \                # UFBoot
    -alrt 1000 \              # SH-aLRT
    -nt AUTO \
    -pre transmission_tree

import pandas as pd  # 导入pandas
import networkx as nx  # 导入networkx画网络图
import matplotlib.pyplot as plt  # 导入matplotlib

# 读取SNP距离矩阵
dist = pd.read_csv("snp_distance_matrix.tsv", sep="\t", index_col=0)

# 读取流行病学数据（采样时间、病房等）
epi = pd.read_csv("epidemiological_data.csv")

# === 基于SNP阈值判断传播关系 ===
# 常用阈值（取决于病原体）：
# 金黄色葡萄球菌(MRSA): <40 SNPs（同一传播链）
# 结核分枝杆菌(MTB): <5-12 SNPs（近期传播）
# 大肠杆菌: <20 SNPs

SNP_THRESHOLD = 15  # 设置SNP阈值

# 构建传播网络
G = nx.Graph()
samples = dist.columns.tolist()
G.add_nodes_from(samples)

for i in range(len(samples)):
    for j in range(i+1, len(samples)):
        snp_dist = dist.iloc[i, j]
        if snp_dist <= SNP_THRESHOLD:  # SNP距离在阈值内
            G.add_edge(samples[i], samples[j], weight=snp_dist)

# 画传播网络图
pos = nx.spring_layout(G)  # 布局
nx.draw(G, pos, with_labels=True, node_size=500,
        node_color='lightblue', font_size=8)
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=6)
plt.title(f"Transmission Network (SNP threshold={SNP_THRESHOLD})")
plt.savefig("transmission_network.pdf")
plt.close()

# 输出传播簇
clusters = list(nx.connected_components(G))
print(f"发现 {len(clusters)} 个传播簇:")
for i, cluster in enumerate(clusters):
    print(f"  簇{i+1}: {', '.join(cluster)}")

# === 结合时间和空间信息 ===
library(ape)      # 系统发育分析
library(ggtree)   # 画树

# 读取树和元数据
tree <- read.tree("transmission_tree.treefile")
meta <- read.csv("epidemiological_data.csv")

# 画带元数据的系统发育树
p <- ggtree(tree) %<+% meta +  # 合并元数据
    geom_tiplab(aes(color = Ward), size = 3) +  # 按病房着色
    geom_tippoint(aes(shape = Source), size = 2) +  # 按来源形状
    theme_tree2() +
    labs(title = "Transmission Phylogeny")

# 添加时间线
p <- p + geom_text(
    aes(label = Collection_Date),
    size = 2, hjust = -0.1
)

ggsave("transmission_phylogeny.pdf", p, width = 12, height = 8)