679 微生物组基因丰度定量¶

# RPKM vs FPKM vs TPM vs Coverage

# RPKM (Reads Per Kilobase per Million mapped reads)
RPKM_i = (reads_i × 10^9) / (total_reads × gene_length_i)
# 先除以总reads（百万），再除以基因长度（千碱基）
# 问题：不同样本的RPKM总和不同，样本间不可直接比较

# FPKM (Fragments Per Kilobase per Million)
# = RPKM的双端测序版本（一对reads算一个fragment）

# TPM (Transcripts Per Million)
# 步骤1: RPK_i = reads_i / gene_length_i (千碱基)
# 步骤2: TPM_i = RPK_i / sum(RPK) × 10^6
# 优点：每个样本的TPM总和=10^6，样本间可直接比较

# Coverage (覆盖度/深度)
Coverage_i = (reads_i × read_length) / gene_length_i
# 平均每个碱基被测了多少次

# ============ 完整的基因丰度定量流程 ============

# 1. 基因预测（从contigs预测ORF）
prodigal \
  -i contigs.fasta \             # 输入组装contigs
  -o genes.gff \                 # GFF格式输出
  -a proteins.faa \              # 蛋白质序列
  -d genes.fna \                 # 核酸序列
  -p meta \                      # 宏基因组模式
  -f gff                        # 输出格式

# 2. 构建非冗余基因目录（去冗余）
cd-hit-est \
  -i all_genes.fna \             # 所有样本的基因合并
  -o gene_catalog.fna \          # 非冗余基因目录
  -c 0.95 \                      # 95%一致性阈值
  -aS 0.9 \                      # 90%短序列覆盖度
  -G 0 \                         # 全局比对
  -M 0 \                         # 不限内存
  -T 16 \                        # 线程数
  -d 0                          # 不截断基因名

echo "基因目录大小:"
grep -c ">" gene_catalog.fna     # 统计基因数量

# 3. 构建比对索引
# 方法1: BWA-MEM2（推荐）
bwa-mem2 index gene_catalog.fna  # 建BWA索引

# 方法2: Bowtie2
bowtie2-build \
  --threads 16 \                 # 线程数
  gene_catalog.fna \             # 基因目录
  gene_catalog_bt2               # 索引前缀

# 方法3: Minimap2（更快）
minimap2 -d gene_catalog.mmi gene_catalog.fna  # 建索引

# 4. Reads比对到基因目录
# 使用BWA-MEM2
bwa-mem2 mem \
  -t 16 \                       # 线程数
  gene_catalog.fna \             # 基因目录
  sample_R1.fastq.gz \           # 正向reads
  sample_R2.fastq.gz \           # 反向reads
  | samtools sort -@ 4 -o sample_sorted.bam  # 排序BAM

samtools index sample_sorted.bam  # 建索引

# 5. 计算基因丰度
# 方法1: featureCounts（推荐，快速准确）
featureCounts \
  -a gene_catalog.gff \          # 基因注释文件
  -o gene_counts.txt \           # 输出计数
  -T 16 \                       # 线程数
  -p \                           # 双端测序
  --countReadPairs \             # 按片段计数
  -B \                           # 两端都比对
  -Q 10 \                       # 最低比对质量
  sample_sorted.bam              # 输入BAM文件

# 方法2: CoverM（覆盖度计算）
coverm genome \
  -b sample_sorted.bam \         # BAM文件
  -m trimmed_mean \              # 方法：截尾均值
  --min-covered-fraction 0.1 \   # 最低覆盖比例
  -o coverage_results.tsv        # 输出

# 方法3: samtools + bedtools
# 计算每个基因的reads数
bedtools coverage \
  -a gene_catalog.bed \          # 基因坐标
  -b sample_sorted.bam \         # BAM文件
  -counts \                      # 只输出计数
  > gene_raw_counts.txt          # 原始计数

# Python实现基因丰度标准化
import numpy as np   # 数值计算
import pandas as pd  # 数据处理

# 1. 读取featureCounts输出
def read_featurecounts(filepath):
    """读取featureCounts输出文件"""
    df = pd.read_csv(filepath, sep='\t', comment='#')  # 读取
    # featureCounts输出列：Geneid, Chr, Start, End, Strand, Length, counts...
    gene_ids = df["Geneid"]              # 基因ID
    lengths = df["Length"]               # 基因长度
    count_cols = [c for c in df.columns if c.endswith(".bam")]  # 计数列
    counts = df[count_cols]              # 提取计数矩阵
    # 清理列名
    counts.columns = [c.replace("_sorted.bam", "") for c in counts.columns]
    counts.index = gene_ids              # 设基因ID为索引
    return counts, lengths

counts, gene_lengths = read_featurecounts("gene_counts.txt")
print(f"基因数: {len(counts)}")
print(f"样本数: {len(counts.columns)}")

# 2. 计算RPKM
def calculate_rpkm(counts, gene_lengths):
    """计算RPKM值"""
    # RPKM = (reads × 10^9) / (total_reads × gene_length)
    total_reads = counts.sum(axis=0)     # 每个样本的总reads
    rpkm = counts.copy().astype(float)   # 初始化
    for sample in counts.columns:
        rpkm[sample] = (counts[sample] * 1e9) / \
                        (total_reads[sample] * gene_lengths)  # RPKM公式
    return rpkm

rpkm = calculate_rpkm(counts, gene_lengths)
print("\nRPKM范围:")
print(rpkm.describe())

# 3. 计算TPM（推荐）
def calculate_tpm(counts, gene_lengths):
    """计算TPM值"""
    # 步骤1: 计算RPK (reads per kilobase)
    rpk = counts.divide(gene_lengths / 1000, axis=0)  # 除以基因长度(kb)
    # 步骤2: 缩放到百万
    tpm = rpk.divide(rpk.sum(axis=0), axis=1) * 1e6   # 每个样本总和=10^6
    return tpm

tpm = calculate_tpm(counts, gene_lengths)
print("\nTPM范围:")
print(tpm.describe())
# 验证：每个样本的TPM总和应该=10^6
print("\n各样本TPM总和（应为1e6）:")
print(tpm.sum(axis=0))

# 4. 计算Coverage
def calculate_coverage(counts, gene_lengths, read_length=150):
    """计算平均覆盖度"""
    # Coverage = (reads × read_length) / gene_length
    coverage = counts.multiply(read_length, axis=0).divide(
        gene_lengths, axis=0)            # 覆盖度公式
    return coverage

cov = calculate_coverage(counts, gene_lengths)

# 5. 过滤低丰度基因
def filter_low_abundance(tpm, min_tpm=1, min_samples=0.1):
    """过滤低丰度基因"""
    n_samples = len(tpm.columns)         # 样本数
    min_n = int(n_samples * min_samples) # 最少出现的样本数
    # 在至少min_samples比例的样本中TPM≥min_tpm
    keep = (tpm >= min_tpm).sum(axis=1) >= min_n  # 过滤条件
    filtered = tpm[keep]                 # 保留的基因
    print(f"过滤前: {len(tpm)} 基因")
    print(f"过滤后: {len(filtered)} 基因")
    print(f"去除: {len(tpm) - len(filtered)} 个低丰度基因")
    return filtered

tpm_filtered = filter_low_abundance(tpm, min_tpm=1, min_samples=0.1)

# 6. 保存结果
tpm_filtered.to_csv("gene_tpm_filtered.tsv", sep='\t')  # 保存TPM表
print("TPM表已保存: gene_tpm_filtered.tsv")

# 将基因丰度汇总到功能注释级别
import pandas as pd  # 数据处理

# 1. 读取功能注释（eggNOG-mapper/KEGG/COG）
def summarize_to_function(tpm, annotation_file, level="KO"):
    """将基因TPM汇总到功能分类级别"""
    # 读取功能注释
    annot = pd.read_csv(annotation_file, sep='\t',
                        index_col=0)     # 基因→功能映射
    # 合并TPM和注释
    tpm_annot = tpm.join(annot[[level]], how="inner")  # 内连接
    # 按功能分组求和
    func_tpm = tpm_annot.groupby(level).sum()  # 汇总到功能级别
    print(f"基因级别: {len(tpm)} 条 → 功能级别: {len(func_tpm)} 条")
    return func_tpm

# KEGG pathway汇总
# ko_tpm = summarize_to_function(tpm, "eggnog_annotations.tsv", level="KEGG_ko")

# 2. KEGG module completeness
def kegg_module_completeness(ko_tpm, module_definitions):
    """计算KEGG模块的完整度"""
    results = {}
    for module, kos in module_definitions.items():
        present_kos = [ko for ko in kos if ko in ko_tpm.index]
        completeness = len(present_kos) / len(kos)  # 完整度
        if completeness > 0:
            avg_abundance = ko_tpm.loc[present_kos].mean().mean()  # 平均丰度
            results[module] = {
                "completeness": completeness,
                "n_present": len(present_kos),
                "n_total": len(kos),
                "avg_tpm": avg_abundance
            }
    return pd.DataFrame(results).T

# 3. 差异功能分析
def differential_function(func_tpm, metadata, group_col,
                           group1, group2):
    """两组间的差异功能分析"""
    from scipy.stats import mannwhitneyu  # Mann-Whitney U检验

    g1_samples = metadata[metadata[group_col] == group1].index
    g2_samples = metadata[metadata[group_col] == group2].index

    results = []
    for func in func_tpm.index:
        vals1 = func_tpm.loc[func, g1_samples]  # 组1值
        vals2 = func_tpm.loc[func, g2_samples]  # 组2值
        stat, pval = mannwhitneyu(vals1, vals2,
                                  alternative='two-sided')  # U检验
        fc = vals1.mean() / (vals2.mean() + 1e-10)  # fold change
        results.append({
            "function": func,
            "mean_g1": vals1.mean(),
            "mean_g2": vals2.mean(),
            "fold_change": fc,
            "log2fc": np.log2(fc + 1e-10),
            "pvalue": pval
        })

    result_df = pd.DataFrame(results)
    result_df["padj"] = result_df["pvalue"] * len(result_df)  # Bonferroni
    result_df["padj"] = result_df["padj"].clip(upper=1.0)
    return result_df.sort_values("pvalue")

# 基因丰度定量流程
1. 基因预测: Prodigal (meta模式)
2. 构建基因目录: CD-HIT-EST (95%一致性)
3. Reads比对: BWA-MEM2 / Minimap2
4. 计数: featureCounts / bedtools
5. 标准化: TPM（推荐）
6. 过滤: TPM≥1 在≥10%样本中
7. 功能汇总: 按KO/COG/Pathway聚合

# 标准化选择
样本间比较: TPM（总和恒定=10^6）
统计模型输入: Raw counts（DESeq2自行标准化）
展示丰度: TPM 或 相对丰度

# 公共基因目录
IGC2: 人肠道非冗余基因集（1504万基因）
UHGG: 统一人肠道基因组（28万基因组+13亿基因）
OM-RGC: 海洋参考基因目录

# 去冗余工具对比
CD-HIT-EST: 经典，适合<10M基因
MMseqs2 linclust: 更快，适合>10M基因
VSEARCH: CD-HIT替代品

# 比对工具对比
BWA-MEM2: 最准确，适合精确定量
Minimap2: 最快，适合大数据集
Bowtie2: 经典，参数丰富