612_泛基因组分析Roary¶

# === conda安装（推荐）===
conda create -n pangenome python=3.8  # 创建环境
conda activate pangenome  # 激活环境

# 安装Roary和Prokka（Roary需要Prokka生成的GFF文件）
conda install -c bioconda roary prokka  # 一起安装

# 验证安装
roary --version  # 查看版本
prokka --version  # 查看版本

# === Docker安装 ===
docker pull sangerpathogens/roary  # Docker镜像

# Phandango在线可视化（推荐）
# 访问 https://jameshadfield.github.io/phandango/

# 安装roary_plots.py的依赖
pip install matplotlib numpy biopython  # Python画图依赖
conda install -c bioconda fasttree  # 快速建树工具

# === 批量用Prokka注释所有基因组 ===
# Roary的输入必须是Prokka生成的GFF3文件

for genome in genomes/*.fasta; do
    name=$(basename "$genome" .fasta)  # 提取文件名

    prokka \
        --kingdom Bacteria \       # 细菌基因组
        --outdir prokka_${name} \  # 输出目录
        --prefix ${name} \         # 输出文件前缀
        --locustag ${name} \       # 基因座位标签
        --cpus 8 \                 # CPU数
        --force \                  # 覆盖已有结果
        ${genome}                  # 输入基因组

    echo "完成注释: ${name}"
done

# 收集所有GFF文件到一个目录
mkdir -p gff_files  # 创建GFF目录
cp prokka_*/*.gff gff_files/  # 复制所有GFF文件

# === 基本运行 ===
roary \
    -p 16 \                    # 16个CPU
    -f roary_output \          # 输出目录
    -e \                       # 生成核心基因多序列比对
    -n \                       # 用MAFFT做比对（推荐）
    -v \                       # 详细输出
    gff_files/*.gff            # 所有GFF文件

# === 高级参数 ===
roary \
    -p 16 \                    # CPU数
    -f roary_advanced \        # 输出目录
    -e -n \                    # 核心基因比对用MAFFT
    -i 95 \                    # BLASTp序列相似性阈值（默认95%）
    -cd 99 \                   # 核心基因定义：99%菌株共有
    -s \                       # 不拆分旁系同源基因
    -g 60000 \                 # 最大基因簇数（大数据集需增大）
    -v \                       # 详细输出
    gff_files/*.gff            # 输入文件

# Roary输出的关键文件
ls roary_output/

# gene_presence_absence.csv    — 最重要！基因在每个菌株的有/无
# gene_presence_absence.Rtab   — R语言友好的二元矩阵
# core_gene_alignment.aln      — 核心基因的多序列比对
# summary_statistics.txt       — 统计摘要
# pan_genome_reference.fa      — 泛基因组参考序列
# number_of_conserved_genes.Rtab  — 累积曲线数据
# number_of_new_genes.Rtab       — 新基因曲线数据
# number_of_genes_in_pan_genome.Rtab — 泛基因组大小曲线
# accessory_binary_genes.fa.newick — 附属基因组的树

# 查看统计摘要
cat roary_output/summary_statistics.txt

# === 方法1：Roary自带的Python脚本 ===
# 先用FastTree建树
FastTree -nt -gtr roary_output/core_gene_alignment.aln \
    > core_gene_tree.nwk  # 基于核心基因比对建NJ树

# 画泛基因组矩阵图
python roary_plots.py \
    core_gene_tree.nwk \                       # 系统发育树
    roary_output/gene_presence_absence.csv     # 基因存在/缺失矩阵
# 输出：pangenome_matrix.png, pangenome_pie.png, pangenome_frequency.png

# === 方法2：Phandango在线可视化 ===
# 1. 访问 https://jameshadfield.github.io/phandango/
# 2. 拖入 core_gene_tree.nwk
# 3. 拖入 gene_presence_absence.csv
# 交互式查看基因有无模式

# === 泛基因组统计分析 ===
library(ggplot2)  # 加载画图包

# 读取泛基因组大小累积数据
pan_data <- read.table(
    "roary_output/number_of_genes_in_pan_genome.Rtab",
    header = FALSE
)

# 读取核心基因组累积数据
core_data <- read.table(
    "roary_output/number_of_conserved_genes.Rtab",
    header = FALSE
)

# === 画泛基因组开放/封闭曲线 ===
n_genomes <- 1:ncol(pan_data)  # 基因组数量

# 计算每个基因组数量下的平均值
pan_means <- colMeans(pan_data)    # 泛基因组平均大小
core_means <- colMeans(core_data)  # 核心基因组平均大小

plot_df <- data.frame(
    genomes = rep(n_genomes, 2),  # 基因组数
    genes = c(pan_means, core_means),  # 基因数
    type = rep(c("Pan-genome", "Core genome"), each = length(n_genomes))
)

ggplot(plot_df, aes(x = genomes, y = genes, color = type)) +
    geom_line(linewidth = 1.2) +  # 画曲线
    geom_point(size = 2) +  # 画点
    labs(
        title = "Pan-genome and Core genome curves",
        x = "Number of genomes",  # x轴
        y = "Number of genes",    # y轴
        color = "Type"            # 图例标题
    ) +
    theme_minimal() +  # 简洁主题
    scale_color_manual(values = c("Core genome" = "blue", 
                                   "Pan-genome" = "red"))

ggsave("pangenome_curve.pdf", width = 10, height = 6)  # 保存

# === 分析基因频率分布 ===
# 读取基因存在/缺失矩阵
gene_pa <- read.csv(
    "roary_output/gene_presence_absence.Rtab",
    sep = "\t", row.names = 1
)

# 计算每个基因在多少个菌株中存在
gene_freq <- rowSums(gene_pa > 0)  # 每个基因的存在频率

# 分类统计
n_strains <- ncol(gene_pa)  # 总菌株数
core <- sum(gene_freq == n_strains)       # 核心基因
soft_core <- sum(gene_freq >= 0.95 * n_strains & gene_freq < n_strains)
shell <- sum(gene_freq >= 0.15 * n_strains & gene_freq < 0.95 * n_strains)
cloud <- sum(gene_freq < 0.15 * n_strains)

cat("核心基因(100%):", core, "\n")
cat("软核心基因(95-99%):", soft_core, "\n")
cat("Shell基因(15-95%):", shell, "\n")
cat("Cloud基因(<15%):", cloud, "\n")
cat("泛基因组总大小:", nrow(gene_pa), "\n")