405_蛋白基因组学Proteogenomics¶

# ========== 蛋白基因组学分析流程示例 ==========

# ========== 1. 从VCF生成突变肽段数据库 ==========
def vcf_to_mutant_peptides(vcf_file, reference_fasta, upstream=10, downstream=10):
    """
    从VCF文件生成突变肽段序列
    原理：取突变位点前后各10个氨基酸，生成突变版本和野生型版本
    """
    # 简化示例，实际使用 pyvcf 或 cyvcf2
    mutations = [
        {"chrom": "chr1", "pos": 1000, "ref": "G", "alt": "A", 
         "gene": "TP53", "codon_change": "R175H"},
        {"chrom": "chr7", "pos": 2000, "ref": "C", "alt": "T", 
         "gene": "KRAS", "codon_change": "G12D"},
    ]

    # 氨基酸单字母代码对照表
    aa_table = {
        "Arg": "R", "His": "H", "Lys": "K", "Asp": "D", "Glu": "E",
        "Gly": "G", "Ala": "A", "Val": "V", "Leu": "L", "Ile": "I",
        "Pro": "P", "Phe": "F", "Trp": "W", "Met": "M", "Ser": "S",
        "Thr": "T", "Cys": "C", "Tyr": "Y", "Asn": "N", "Gln": "Q"
    }

    peptides = []
    for mut in mutations:
        # 解析氨基酸变化（例如 "R175H" = 第175位R突变为H）
        ref_aa = mut["codon_change"][0]      # 野生型氨基酸
        mut_aa = mut["codon_change"][-1]     # 突变氨基酸
        position = int(mut["codon_change"][1:-1])  # 位置

        # 模拟蛋白质序列（实际从FASTA文件提取）
        mock_protein = "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP"

        # 确保位置在蛋白范围内
        if position <= len(mock_protein):
            start = max(0, position - upstream - 1)
            end = min(len(mock_protein), position + downstream)

            # 提取突变肽段（包含突变位点）
            mutant_peptide = mock_protein[start:position-1] + mut_aa + mock_protein[position:end]
            wild_peptide = mock_protein[start:end]  # 野生型对照

            peptides.append({
                "gene": mut["gene"],
                "mutation": mut["codon_change"],
                "mutant_peptide": mutant_peptide,
                "wild_type_peptide": wild_peptide
            })

    return peptides

mutations = vcf_to_mutant_peptides("tumor.vcf", "hg38.fa")
for m in mutations:
    print(f"基因: {m['gene']}, 突变: {m['mutation']}")
    print(f"  突变肽段: {m['mutant_peptide']}")
    print(f"  野生型: {m['wild_type_peptide']}")

# ========== 2. 生成自定义搜库数据库（FASTA格式） ==========
def create_custom_database(mutations, output_fasta="custom_db.fasta"):
    """
    将突变肽段写入FASTA格式数据库，用于质谱搜库
    格式：>蛋白ID 描述信息
          氨基酸序列
    """
    with open(output_fasta, "w") as f:
        for i, mut in enumerate(mutations):
            # 突变肽段条目
            f.write(f">MUT_{mut['gene']}_{mut['mutation']} mutant_peptide\n")
            f.write(f"{mut['mutant_peptide']}\n")
            # 野生型对照条目
            f.write(f">WT_{mut['gene']}_{mut['mutation']} wild_type_peptide\n")
            f.write(f"{mut['wild_type_peptide']}\n")
    print(f"数据库已写入: {output_fasta}")

create_custom_database(mutations)

# ========== 3. 分析质谱搜库结果（PSM过滤） ==========
import pandas as pd
import numpy as np

# 模拟质谱搜库输出（PSM = Peptide-Spectrum Match，肽段-谱图匹配）
psm_results = pd.DataFrame({
    "peptide_sequence": ["SVEPPLSQETF", "LMLSPDDIEQW", "MEEPQSDPSVEP", "KLLPENNVLS"],
    "protein_id": ["MUT_TP53_R175H", "WT_TP53_R175H", "MUT_KRAS_G12D", "WT_KRAS_G12D"],
    "score": [85.2, 72.1, 91.5, 68.3],       # 搜库打分
    "pep": [0.001, 0.01, 0.0005, 0.02],       # Posterior Error Probability（越小越好）
    "unique_peptide": [True, True, True, True]  # 是否为蛋白特异肽段
})

# ========== 4. 过滤高置信度突变肽段 ==========
# FDR < 1%：PEP < 0.01 且分数足够高
confident_psm = psm_results[
    (psm_results["pep"] < 0.01) &          # 假阳性率 < 1%
    (psm_results["score"] > 75) &           # 打分阈值
    (psm_results["protein_id"].str.startswith("MUT_"))  # 只保留突变肽段
]

print("\n高置信度突变肽段:")
print(confident_psm[["peptide_sequence", "protein_id", "score", "pep"]].to_string(index=False))

# ========== 5. 统计各类别蛋白的分布 ==========
psm_results["category"] = psm_results["protein_id"].apply(
    lambda x: "Mutant" if x.startswith("MUT_") else "Wild-type"
)
category_counts = psm_results.groupby("category").size()
print("\n各类别PSM数量:")
print(category_counts)