733. 微生物组生长速率估计GRiD（Growth Rate InDex）¶

慢速生长的细菌：
oriC ████████████████████████████████ ter
覆盖度很均匀，PTR ≈ 1

快速生长的细菌：
oriC ████████████████████████████████ ter
████████████                          
████████                              
覆盖度呈梯度分布，PTR >> 1

# 方法一：conda安装（推荐）
conda create -n growth python=3.8  # 创建专用环境
conda activate growth  # 激活环境
conda install -c bioconda grid  # 从bioconda安装GRiD

# 方法二：GitHub安装
git clone https://github.com/ohlab/GRiD.git  # 克隆仓库
cd GRiD  # 进入目录
bash setup.sh  # 运行安装脚本

# 依赖检查
which bowtie2  # 检查bowtie2是否安装（比对工具）
which samtools  # 检查samtools是否安装（BAM处理）
which bedtools  # 检查bedtools是否安装（区间操作）
Rscript -e "library(dplyr)"  # 检查R包是否可用

# ===== 模式一：使用参考基因组数据库 =====
# 适合有已知菌种的分析

# 准备参考基因组目录
mkdir -p grid_database/  # 创建数据库目录

# 将参考基因组fasta文件放入目录
cp reference_genomes/*.fna grid_database/  # 复制参考基因组文件

# 构建GRiD数据库
update_database.sh -d grid_database/ \  # 输入：参考基因组目录
                   -g grid_db  # 输出：GRiD数据库前缀名

# ===== 模式二：使用MAG =====
# 适合没有参考基因组的未知菌种

# 直接用MAG作为参考
# MAG文件放在一个目录里即可
ls mags_directory/  # 确认MAG文件都在这个目录
# bin.1.fa  bin.2.fa  bin.3.fa ...

# ===== 单样本模式 =====
grid single \
  -r reads_R1.fastq.gz \  # 输入：正向reads（质控后的）
  -e fastq \  # 输入文件格式
  -d grid_db \  # 参考数据库路径
  -o output_single/ \  # 输出目录
  -n 8  # 使用8个线程加速

# ===== 多样本模式（推荐） =====
grid multiplex \
  -r reads_directory/ \  # 输入：包含多个样本reads的目录
  -e fastq \  # 输入文件格式
  -d grid_db \  # 参考数据库路径
  -o output_multi/ \  # 输出目录
  -n 16 \  # 使用16个线程
  -p  # 输出每个样本的PDF图

# ===== 使用MAG的模式 =====
grid single \
  -r reads_R1.fastq.gz \  # 输入：原始reads
  -e fastq \  # 格式
  -g mags_directory/ \  # 使用MAG目录替代参考数据库
  -o output_mag/ \  # 输出目录
  -n 8  # 线程数

# GRiD输出文件结构
# output/
# ├── sample1.GRiD.txt           # 主结果文件
# ├── sample1.GRiD.pdf           # 覆盖度分布图
# └── sample1.GRiD.filtered.txt  # 过滤后的可靠结果

# 查看结果
cat output_single/sample1.GRiD.txt
# Genome             GRiD   Coverage   Species_heterogeneity
# Ecoli_K12          1.82   15.3       0.12
# Bfragilis          1.15   8.7        0.05
# Lactobacillus_r    2.34   22.1       0.08

# ===== 用Python可视化GRiD结果 =====
import pandas as pd  # 导入数据处理包
import matplotlib.pyplot as plt  # 导入绑图包
import seaborn as sns  # 导入增强绑图包
import numpy as np  # 导入数值计算包

# 读取多个样本的GRiD结果
import glob  # 导入文件匹配模块

results = []  # 存储所有结果
for f in glob.glob("output_multi/*.GRiD.txt"):  # 遍历所有结果文件
    sample = f.split("/")[-1].replace(".GRiD.txt", "")  # 提取样本名
    df = pd.read_csv(f, sep="\t")  # 读取结果表
    df["Sample"] = sample  # 添加样本名列
    results.append(df)  # 添加到结果列表

all_results = pd.concat(results)  # 合并所有结果

# 过滤：只保留覆盖度>=1x且GRiD值合理的结果
filtered = all_results[
    (all_results["Coverage"] >= 1) &  # 覆盖度至少1x才可靠
    (all_results["GRiD"] >= 1) &  # GRiD值不能小于1（物理意义上不可能）
    (all_results["GRiD"] <= 10) &  # GRiD值太大通常不可靠
    (all_results["Species_heterogeneity"] < 0.3)  # 种内异质性太高说明不是纯菌
]

# 绘制生长速率热图
pivot = filtered.pivot_table(
    index="Genome",  # 行：菌种名
    columns="Sample",  # 列：样本名
    values="GRiD",  # 值：生长速率
    aggfunc="mean"  # 聚合：取平均
)

plt.figure(figsize=(14, 8))  # 设置画布
sns.heatmap(pivot, cmap="RdYlBu_r", center=1.5,  # 绘制热图，红色=快速生长
            annot=True, fmt=".2f",  # 显示数值
            vmin=1, vmax=4)  # 设置色标范围
plt.title("Bacterial Growth Rates (GRiD) Across Samples")  # 标题
plt.xlabel("Sample")  # x轴标签
plt.ylabel("Bacterial Species")  # y轴标签
plt.tight_layout()  # 自动调整
plt.savefig("growth_rate_heatmap.png", dpi=300)  # 保存图片
plt.show()  # 显示

# 绘制分组箱线图（疾病vs健康）
# 假设样本名包含分组信息
filtered["Group"] = filtered["Sample"].apply(
    lambda x: "T2D" if "T2D" in x else "Healthy"  # 根据样本名判断分组
)

plt.figure(figsize=(10, 6))
top_species = filtered.groupby("Genome")["GRiD"].mean().nlargest(10).index  # 取Top10菌种
plot_data = filtered[filtered["Genome"].isin(top_species)]  # 筛选数据

sns.boxplot(data=plot_data, x="Genome", y="GRiD", hue="Group")  # 绘制箱线图
plt.xticks(rotation=45, ha="right")  # 旋转x轴标签
plt.ylabel("Growth Rate (GRiD)")  # y轴标签
plt.title("Growth Rate Comparison: T2D vs Healthy")  # 标题
plt.tight_layout()
plt.savefig("growth_rate_comparison.png", dpi=300)  # 保存
plt.show()

# 安装CoPTR
pip install coptr  # 安装CoPTR

# 第一步：建立参考索引
coptr index \
  reference_genomes/ \  # 输入：参考基因组目录
  coptr_ref  # 输出：索引前缀

# 第二步：比对reads到参考基因组
coptr extract \
  sample.bam \  # 输入：比对好的BAM文件
  coptr_ref \  # 参考索引
  sample.pkl  # 输出：提取的覆盖度信息

# 第三步：估算生长速率
coptr estimate \
  extract_output/ \  # 输入：所有样本的pkl文件目录
  coptr_results.csv  # 输出：生长速率结果表

# ===== GRiD速查 =====

# 安装
conda install -c bioconda grid

# 建库
update_database.sh -d genomes/ -g grid_db

# 运行（单样本）
grid single -r reads.fq.gz -e fastq -d grid_db -o output/ -n 8

# 运行（多样本）
grid multiplex -r reads_dir/ -e fastq -d grid_db -o output/ -n 16

# 结果过滤标准
# Coverage >= 1x
# GRiD >= 1 且 <= 10
# Species_heterogeneity < 0.3

# ===== CoPTR速查 =====
pip install coptr
coptr index ref_genomes/ coptr_ref
coptr extract sample.bam coptr_ref sample.pkl
coptr estimate extract_dir/ results.csv

# ===== 关键文献 =====
# GRiD: Emiola & Oh, 2018, Nature Communications
# iRep: Brown et al., 2016, Nature Biotechnology
# CoPTR: Kundu et al., 2022, Genome Biology