药物基因组学分析¶

一句话概述¶

药物基因组学(Pharmacogenomics)研究基因变异如何影响药物反应——白话说就是：同一种药，有人吃了有效，有人吃了没用甚至中毒，原因就在基因里。药物基因组学帮你找到"哪些基因决定了药物效果"。

核心知识点表格¶

知识点	说明
GDSC	Genomics of Drug Sensitivity in Cancer，最大的癌症药物敏感性数据库
CCLE	Cancer Cell Line Encyclopedia，Broad研究所的细胞系百科全书
CellMinerCDB	NCI跨数据库药物基因组学分析平台（v2.2，2025年更新）
IC50	药物抑制50%细胞生长所需浓度，越低说明越敏感
AUC	剂量-反应曲线下面积，衡量总体药物敏感性
药物靶点预测	通过基因组数据推测药物作用靶点
耐药机制	基因突变、表达变化等导致的药物抗性
生物标志物	预测药物反应的基因/蛋白标记
oncoPredict	R包，用细胞系数据预测临床样本药物反应
CellHit	2025年ML方法，从GDSC/PRISM数据预测药物敏感性

各步骤详解¶

第一步：理解药物敏感性数据¶

白话解释： 科学家在上千种癌细胞系上测试了几百种药物，记录每种"细胞+药物"组合的效果（IC50/AUC）。同时对这些细胞系做了全基因组测序、RNA-seq等。这样就能找到"基因A突变 → 对药物X敏感"的关系。

主要数据库对比：
GDSC:  ~1000个细胞系 × ~500种药物 × 基因组/表达/甲基化
CCLE:  ~1500个细胞系 × 基因组/表达/蛋白组
CTRP:  ~860个细胞系 × ~500种化合物
PRISM: ~900个细胞系 × ~4000种药物（最大药物库）

第二步：从GDSC获取药物敏感性数据¶

白话解释： GDSC网站可以直接下载药物IC50数据和对应的基因组数据。

# ===== 方法1：使用R包PharmacoGx =====
# install.packages("BiocManager")
# BiocManager::install("PharmacoGx")

library(PharmacoGx)                      # 加载PharmacoGx

# 查看可用数据集
availablePSets()                          # 列出所有可用的药理学数据集

# 下载GDSC数据（可能需要几分钟）
gdsc <- downloadPSet("GDSC2")           # 下载GDSC2数据集

# 查看基本信息
show(gdsc)                              # 数据集概览
drugNames(gdsc)[1:10]                    # 前10种药物名
cellNames(gdsc)[1:10]                    # 前10个细胞系名

# 提取药物敏感性数据
drug_sensitivity <- summarizeSensitivityProfiles(
  gdsc,                                  # 数据集
  sensitivity.measure = "ic50_published", # 使用发表的IC50
  summary.stat = "median"                # 中位数汇总
)

# 提取基因表达数据
gene_expression <- summarizeMolecularProfiles(
  gdsc,                                  # 数据集
  mDataType = "rna",                     # RNA表达数据
  summary.stat = "median"                # 中位数汇总
)

# ===== 方法2：直接从GDSC网站下载 =====
# https://www.cancerrxgene.org/downloads/bulk_download

# 下载药物敏感性数据
# wget https://www.cancerrxgene.org/downloads/drug_data_all.csv
# 下载基因表达数据
# wget https://www.cancerrxgene.org/downloads/gene_expression_all.csv

第三步：药物敏感性预测（oncoPredict）¶

白话解释： oncoPredict用细胞系的药物反应数据训练模型，然后预测临床肿瘤样本（如TCGA）对各种药物的敏感性——就像用"已知答案"训练AI，再让它预测"未知答案"。

# install.packages("oncoPredict")
library(oncoPredict)                     # 加载oncoPredict

# 准备训练数据（GDSC细胞系）
# 1. 药物敏感性矩阵：行=细胞系，列=药物，值=IC50
# 2. 基因表达矩阵：行=基因，列=细胞系
# 这些数据可以从GDSC下载或用PharmacoGx提取

# 准备测试数据（临床样本，如TCGA）
# 基因表达矩阵：行=基因，列=样本

# 运行药物敏感性预测
calcPhenotype(
  trainingExprData = gdsc_expression,    # 训练集表达矩阵（GDSC）
  trainingPtype = gdsc_ic50,             # 训练集药物敏感性（IC50）
  testExprData = tcga_expression,        # 测试集表达矩阵（TCGA）
  batchCorrect = "eb",                   # 批次校正方法（经验贝叶斯）
  powerTransformPhenotype = TRUE,        # 对IC50做幂变换
  removeLowVaryingGenes = 0.2,           # 移除低变异基因（20%）
  minNumSamples = 10,                    # 每种药物最少10个细胞系
  selection = -1,                        # 使用所有公共基因
  printOutput = TRUE,                    # 打印进度
  pcr = FALSE,                           # 不使用PCR
  report_pc = FALSE                      # 不报告PC
)
# 输出：DrugPredictions.csv，每个TCGA样本的预测IC50

第四步：药物靶点关联分析¶

白话解释： 找出"基因X的表达/突变与药物Y的敏感性显著相关"——如果某基因高表达的细胞系对某药物特别敏感，这个基因可能就是药物的靶点或biomarker。

library(ggplot2)                         # 可视化

# 读取药物敏感性和基因表达数据
# drug_data: 行=细胞系, 列=药物, 值=IC50
# expr_data: 行=基因, 列=细胞系

# 计算每个基因与药物敏感性的相关性
drug_name <- "Erlotinib"                 # 目标药物
drug_ic50 <- drug_data[, drug_name]      # 提取该药物的IC50

# 计算所有基因与IC50的Spearman相关
correlations <- apply(expr_data, 1, function(gene_expr) {
  # 取有效配对
  valid <- !is.na(gene_expr) & !is.na(drug_ic50)
  if (sum(valid) < 20) return(c(cor = NA, pval = NA))  # 样本太少跳过

  test <- cor.test(                      # 相关性检验
    gene_expr[valid],                    # 基因表达
    drug_ic50[valid],                    # 药物IC50
    method = "spearman"                  # Spearman秩相关
  )
  c(cor = test$estimate, pval = test$p.value)  # 返回相关系数和p值
})

# 整理结果
cor_results <- data.frame(
  gene = rownames(expr_data),            # 基因名
  correlation = correlations["cor", ],   # 相关系数
  pvalue = correlations["pval", ]        # p值
)
cor_results$fdr <- p.adjust(cor_results$pvalue, method = "BH")  # FDR校正

# 筛选显著相关的基因
significant <- cor_results[cor_results$fdr < 0.05, ]  # FDR < 0.05
significant <- significant[order(abs(significant$correlation), decreasing = TRUE), ]
head(significant, 20)                    # 查看top 20关联基因

# 火山图可视化
ggplot(cor_results, aes(x = correlation, y = -log10(pvalue))) +
  geom_point(alpha = 0.3) +              # 散点图
  geom_point(data = significant, color = "red") +  # 显著基因标红
  theme_minimal() +                      # 简洁主题
  labs(title = paste("Gene associations with", drug_name))  # 标题

第五步：CellMinerCDB跨数据库分析¶

白话解释： CellMinerCDB让你同时查询GDSC、CCLE、CTRP等多个数据库，比较药物-基因关系的一致性。

# CellMinerCDB主要通过Web界面使用
# 网址：https://discover.nci.nih.gov/cellminercdb/

# 也可以用R API访问NCI-60数据
# install.packages("rcellminer")
library(rcellminer)                      # 加载rcellminer

# 获取NCI-60药物活性数据
drugAct <- getAct(rcellminerData::drugData)  # 药物活性矩阵

# 获取NCI-60基因表达数据
geneExpr <- getAllFeatureData(rcellminerData::molData)[["exp"]]  # 表达矩阵

# 计算特定药物与基因的相关性
drug_id <- "94600"                       # 药物NSC编号（如5-FU）
gene_name <- "TYMS"                      # 基因名（胸苷酸合成酶）

# 提取数据并计算相关
drug_vals <- drugAct[drug_id, ]          # 药物活性值
gene_vals <- geneExpr[gene_name, ]       # 基因表达值

cor.test(drug_vals, gene_vals, method = "pearson", use = "complete.obs")

常见报错与解决¶

报错	原因	解决方案
`downloadPSet: timeout`	数据集太大下载超时	设置options(timeout=600)，或用断点续传
`oncoPredict: no common genes`	训练集和测试集基因名格式不一致	统一基因名格式（Ensembl ID vs Gene Symbol）
`calcPhenotype: singular matrix`	基因表达矩阵有共线性	增加removeLowVaryingGenes比例
`PharmacoGx: assay not found`	数据集版本不匹配	用availablePSets()确认可用数据集名称
`批次效应严重`	细胞系和临床样本差异大	使用batchCorrect="eb"或"combat"

速查表¶

# ========== 药物基因组学速查 ==========

# 1. 获取GDSC数据
library(PharmacoGx)
gdsc <- downloadPSet("GDSC2")
drug_data <- summarizeSensitivityProfiles(gdsc, "ic50_published")
expr_data <- summarizeMolecularProfiles(gdsc, "rna")

# 2. 预测临床样本药物敏感性
library(oncoPredict)
calcPhenotype(trainingExprData, trainingPtype, testExprData, batchCorrect="eb")

# 3. 药物-基因相关分析
cor.test(gene_expression, drug_ic50, method="spearman")

# 4. 关键数据库
# GDSC: https://www.cancerrxgene.org/
# CCLE: https://depmap.org/portal/
# CellMinerCDB: https://discover.nci.nih.gov/cellminercdb/
# PRISM: https://depmap.org/repurposing/

面试高频问题¶

Q1: 什么是IC50？AUC和IC50有什么区别？

IC50是抑制50%细胞生长的药物浓度，数值越低说明药物越有效。AUC是剂量-反应曲线下面积，综合反映所有浓度下的效果。AUC更稳健（不依赖曲线拟合的外推），但IC50更直观。GDSC同时提供两者。

Q2: 用细胞系数据预测临床药物反应可靠吗？

有一定参考价值但不完全可靠。局限性：(1) 细胞系缺少肿瘤微环境（免疫细胞、基质细胞）；(2) 长期培养导致基因组漂移；(3) 2D培养vs体内3D环境差异大。改进方向：用Celligner等工具对齐细胞系和临床样本的转录组，或用类器官(organoid)数据替代。

Q3: 药物基因组学有哪些主要分析方法？

(1) 关联分析：基因表达/突变与IC50的相关性；(2) 机器学习预测：随机森林、深度学习（CellHit、DrugS）从多组学数据预测药物反应；(3) 网络药理学：构建药物-靶点-疾病网络；(4) 孟德尔随机化：用遗传变异作为工具变量推断因果关系。

Q4: 什么是拷贝数对药物敏感性的影响？

基因拷贝数扩增可能导致靶蛋白过表达（如HER2扩增对曲妥珠单抗敏感），也可能导致耐药（如MYC扩增对多种化疗药耐药）。分析时需要整合拷贝数、表达、突变等多层组学数据。

Q5: 如何验证预测的药物靶点？

计算预测后需要实验验证：(1) 体外验证：siRNA/shRNA敲低候选靶基因，看是否影响药物敏感性；(2) CRISPR筛选：用MAGeCK等分析基因敲除与药物反应的关系；(3) PDX模型：在患者来源异种移植模型中验证；(4) 临床关联：对比临床试验中的基因型和药物反应数据。