免疫治疗生物标志物分析¶

体细胞突变 → 氨基酸变化 → HLA分型 → MHC-肽段结合预测 → 新抗原候选

# === CIBERSORT 免疫细胞浸润分析 ===

# 安装和加载
# CIBERSORT 需要从官网下载 R 源码（需注册）
source("CIBERSORT.R")                  # 加载CIBERSORT函数

# 1. 准备表达矩阵（TPM 或 FPKM，行=基因，列=样本）
expr_matrix <- read.table(
  "expression_tpm.txt",               # TPM表达矩阵
  header = TRUE, row.names = 1,       # 有表头和行名
  sep = "\t"                           # 制表符分隔
)

# 2. 运行 CIBERSORT（LM22 参考签名）
results <- CIBERSORT(
  sig_matrix = "LM22.txt",            # 22种免疫细胞参考签名
  mixture_file = "expression_tpm.txt", # 表达矩阵
  perm = 100,                          # 置换次数（统计检验）
  QN = TRUE                            # 分位数归一化
)

# 3. 提取结果（前22列是细胞比例）
immune_fractions <- results[, 1:22]    # 22种免疫细胞比例
p_values <- results[, "P-value"]       # 统计显著性

# 4. 可视化
library(ggplot2)                       # 绑图库
library(tidyr)                         # 数据整理

# 转为长格式
long_df <- as.data.frame(immune_fractions)
long_df$sample <- rownames(long_df)    # 添加样本列
long_df <- pivot_longer(long_df,       # 转长格式
  cols = -sample,                      # 除样本列外
  names_to = "cell_type",             # 细胞类型
  values_to = "fraction"              # 比例
)

# 绘制堆积柱状图
ggplot(long_df, aes(x = sample, y = fraction, fill = cell_type)) +
  geom_bar(stat = "identity") +        # 堆积柱状图
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  # X轴旋转
  labs(title = "肿瘤免疫细胞浸润分析",
       y = "细胞比例", x = "样本")

# === 新抗原预测流程 ===

import subprocess                      # 系统命令

# 1. HLA 分型（从 WES BAM 推断）
def run_optitype(bam_path, output_dir):
    """用 OptiType 进行 HLA-I 分型"""
    cmd = (
        f"OptiTypePipeline.py "        # OptiType主程序
        f"-i {bam_path} "             # 输入BAM
        f"--dna "                      # DNA模式
        f"-o {output_dir} "           # 输出目录
        f"-v"                          # 详细输出
    )
    subprocess.run(cmd, shell=True)    # 执行

# 2. 从突变生成肽段序列
def generate_peptides(vcf_path, ref_fasta, peptide_length=9):
    """从 VCF 突变生成突变肽段（简化示例）"""
    # 实际使用 pvacseq 或 vaxrank 等专用工具
    cmd = (
        f"pvacseq run "               # pVACseq 新抗原预测
        f"{vcf_path} "                # 输入VCF（需VEP注释）
        f"sample_name "               # 样本名
        f"HLA-A*02:01,HLA-B*07:02 "  # HLA类型（从OptiType获取）
        f"NetMHCpan "                 # 结合预测方法
        f"pvacseq_output/ "           # 输出目录
        f"-e1 8,9,10,11 "            # MHC-I肽段长度
        f"--iedb-install-directory /opt/iedb "  # IEDB路径
        f"-t 8"                        # 线程数
    )
    subprocess.run(cmd, shell=True)    # 执行

# 3. 筛选高亲和力新抗原
import pandas as pd
neoantigen_df = pd.read_csv("pvacseq_output/MHC_Class_I/sample.filtered.tsv",
                            sep="\t")

# 筛选条件
strong_binders = neoantigen_df[
    (neoantigen_df["Median MT IC50"] < 500) &     # 突变肽结合力 <500nM
    (neoantigen_df["Median WT IC50"] > 500) &      # 野生型结合力 >500nM
    (neoantigen_df["Tumor DNA VAF"] > 0.1)         # 突变VAF >10%
]
print(f"强结合新抗原候选数: {len(strong_binders)}")