GDSC(Genomics of Drug Sensitivity in Cancer)数据库由英国桑格研究院开发，收集肿瘤细胞对药物的敏感度和反应。是目前肿瘤细胞药物敏感性的最大公共资源，收录1000多种癌细胞系中的621个抗癌化合物的作用广泛，涉及24种通路。2023年10月更新到了8.5 version。
1. 数据集介绍
ꔷ GDSC内部由和
两个数据库组成：
1）主要收录的是2010-2015年间的测序和试验结果。最新版本的GDSC1包含有970个细胞系，403个抗癌化合物，以及333292个药物反应IC50值和AUC值。（GDSC1 的许多实验已在 GDSC2 中重复，如果存在重复的 IC50，官方建议使用 GDSC2 的结果）
2）主要收录的是2015至今的测序和试验的结果。最新版本的GDSC2包含有969个细胞系，297个抗癌化合物，以及243466个药物反应IC50值。
ꔷ IC50（half maximal inhibitory concentration），半数抑制浓度， 指在一定条件下，化合物或药物能够抑制生物过程或活性的浓度，使得生物过程或活性被抑制50%。通常，IC50用于评估药物的抗生物活性。IC50值越低，说明药物越强效，因为它可以在更低的浓度下抑制目标生物分子的活性。
ꔷ AUC（area under concentration-time curve），血药浓度-时间曲线下面积。代表药物的生物利用度（药物活性成分从制剂释放吸收进入全身循环的程度和速度），AUC大则生物利用度高，反之则低。AUC0→∞指药物从零时间至所有原形药物全部消除这一段时间的药时曲线下总面积，反映药物进入血循环的总量。总的来说就是较低的AUC值表明细胞对治疗的敏感性增加。
ꔷ EC50（half maximal effective concentration），半数有效浓度，指在给定实验条件下，药物或化合物能够产生特定效应（例如细胞激活、酶活性增加等）50%的浓度。通常，EC50用于评估药物对疾病相关分子或细胞的激活或促进作用。较低的EC50值表示药物对目标生物过程的激活效果更好。
★ IC50和EC50的主要区别在于它们描述的效应类型不同。IC50描述了抑制效应的浓度，而EC50描述了激活效应的浓度。
2. 数据下载
若想研究基因表达与药物敏感性的相关性，需要获取基因在不同细胞系中的表达谱和每个细胞系中各种药物的IC50或者AUC值。
1）表达谱信息：Download from GDSC1000 resourse > Dataset（"EXP"-"Preprocessed Cell-lines"-"RMA normalised expression data for cell-lines"-"RMA normalised basal expression profiles for all the cell-lines"-"Pathway Activity Scores")
第一列和第二列为基因名，后续为各个细胞系中对应的基因表达谱数据
2）细胞系信息：Cell-line-annotation
3）药物信息：Compounds-annotation
4）每种药物在每种细胞中的IC50/AUC值：GDSC1-dataset；GDSC2-dataset
3. oncoPredict包分析
1）载入训练数据：
oncoPredict包安装后需要配套的训练数据：https://osf.io/c6tfx/
# 安装并加载所需的R包
# install.packages("oncoPredict")
# BiocManager::install("TCGAbiolinks")
library(oncoPredict)
library(data.table)
library(gtools)
library(reshape2)
library(ggpubr)
library(TCGAbiolinks)
library(dplyr)
library(tidyverse)
dir='/data/shumin/GDSC/DataFiles/Training Data/'
dir(dir)
## [1] "CTRP2_Expr (RPKM, not log transformed).rds" "CTRP2_Expr (TPM, not log transformed).rds" "CTRP2_Res.rds"
## [4] "GDSC1_Expr (RMA Normalized and Log Transformed).rds" "GDSC1_Res.rds" "GDSC2_Expr (RMA Normalized and Log Transformed).rds"
## [7] "GDSC2_Res.rds"
# 读取训练集表达数据
GDSC2_exp = readRDS(file=file.path(dir,'GDSC2_Expr (RMA Normalized and Log Transformed).rds'))
dim(GDSC2_exp)
## [1] 17419 805
GDSC2_exp[1:4,1:4]
## COSMIC_906826 COSMIC_687983 COSMIC_910927 COSMIC_1240138
## TSPAN6 7.632023 7.548671 8.712338 7.797142
## TNMD 2.964585 2.777716 2.643508 2.817923
## DPM1 10.379553 11.807341 9.880733 9.883471
## SCYL3 3.614794 4.066887 3.956230 4.063701
# 读取训练集药物反应数据
GDSC2_drug = readRDS(file = file.path(dir,"GDSC2_Res.rds"))
GDSC2_drug <- exp(GDSC2_drug) #下载的数据是被log转换过的，逆转回去
dim(GDSC2_drug)
## [1] 805 198
GDSC2_drug[1:4,1:4]
## Camptothecin_1003 Vinblastine_1004 Cisplatin_1005 Cytarabine_1006
## COSMIC_906826 0.3158373 0.208843106 1116.05899 18.5038719
## COSMIC_687983 0.2827342 0.013664227 26.75839 16.2943594
## COSMIC_910927 0.0295671 0.006684071 12.09379 0.3387418
## COSMIC_1240138 7.2165789 NA NA NA
identical(rownames(GDSC2_drug),colnames(GDSC2_exp))
## [1] TRUE （二者的样本名称对应）
ꔷ GDSC2_exp是标准的表达矩阵格式，行是基因，列是细胞系，一共17419个基因，805个细胞系
ꔷ GDSC2_drug是每个细胞系对每个药物的IC50值，行是细胞系，列是药物，一共805种细胞系，198个药物
2）加载自己待预测的表达数据
此处数据通过UCSC Xena（http://xenabrowser.net）网站进行下载
# 读取临床信息
phe = fread("/data/shumin/GDSC/BRCA/TCGA-BRCA.GDC_phenotype.tsv.gz", header = T, sep = '\t',data.table = F)
# 读取表达数据
BRCA.fkpm = fread("/data/shumin/GDSC/BRCA/TCGA-BRCA.htseq_fpkm.tsv.gz",header = T, sep = '\t',data.table = F)
# 读取探针信息
BRCA.pro = fread("/data/shumin/GDSC/BRCA/gencode.v22.annotation.gene.probeMap",header = T, sep = '\t',data.table = F)
# 提取前两列，进行转换
BRCA.pro = BRCA.pro[, c(1,2)]
BRCA.fkpm.pro = merge(x = BRCA.pro, y = BRCA.fkpm, by.y = "Ensembl_ID", by.x = "id" )
dim(BRCA.fkpm.pro)
## [1] 60483 1219
# 使用aggregate函数对gene列中的相同基因进行合并（一个探针可能对应多个gene）
BRCA.fkpm.pro = distinct(BRCA.fkpm.pro,gene,.keep_all = T)
dim(BRCA.fkpm.pro)
## [1] 58387 1219
#把表达矩阵中gene列变成行名
BRCA.fkpm.pro <- column_to_rownames(BRCA.fkpm.pro, "gene")
#把Ensembl_ID列去掉
BRCA.fkpm.pro = BRCA.fkpm.pro[,-1]
BRCA.fkpm.pro[1:4, 1:4]
## TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A1EU-11A TCGA-A8-A06X-01A
## TSPAN6 1.707383 4.06802095 4.269338 4.26380546
## TNMD 0.000000 0.08710664 2.420087 0.04326277
## DPM1 5.393644 5.15791344 4.753093 4.73014799
## SCYL3 2.353152 1.73310574 2.509152 2.45011821
3）预测药物反应
calcPhenotype(trainingExprData = GDSC2_exp,
trainingPtype = GDSC2_drug,
testExprData = as.matrix(BRCA.fkpm.pro), # 需要matrix
batchCorrect = 'eb', # 校正参数，"eb" for array, "standardize" for rnaseq
#IC50是对数转换的，所以表达矩阵也用对数转换过的
powerTransformPhenotype = F,
minNumSamples = 20,
printOutput = T,
removeLowVaryingGenes = 0.2,
removeLowVaringGenesFrom = "homogenizeData"
)
result <- read.csv("/data/shumin/GDSC/BRCA/calcPhenotype_Output/DrugPredictions.csv", header = T , stringsAsFactors = F ,check.names = F)
names(result)[1] <- "submitter_id.samples"
result[1:4, 1:4]
## submitter_id.samples Camptothecin_1003 Vinblastine_1004 Cisplatin_1005
## 1 TCGA-E9-A1NI-01A 2.0801021 1.6722196 249.18844
## 2 TCGA-A1-A0SP-01A -1.8898854 -1.0342157 -195.56903
## 3 TCGA-BH-A1EU-11A -0.4090808 -0.1844808 50.07022
## 4 TCGA-A8-A06X-01A 1.5502492 1.0955549 268.94364
# 取临床数据一列，分组示例
phe_test <- phe[,c("submitter_id.samples", "age_at_initial_pathologic_diagnosis")]
phe_test$age <- ifelse(as.numeric(phe_test$age_at_initial_pathologic_diagnosis) >= 60, '>60',
ifelse(as.numeric(phe_test$age_at_initial_pathologic_diagnosis) >= 35, "35-60", "<35"))
phe_test$age <- factor(phe_test$age, levels = c('<35','35-60','>60'))
phe_test <- na.omit(phe_test)
phe_test[1:3,1:3]
## submitter_id.samples age_at_initial_pathologic_diagnosis age
## 1 TCGA-A2-A0CY-01A 63 >60
## 2 TCGA-B6-A40B-01A 76 >60
## 3 TCGA-AO-A0J8-01A 61 >60
# 合并数据
BRCA_FPKM_drug <- result %>% inner_join(phe_test, "submitter_id.samples")
# 生成两两之间的list
group = levels(factor(phe_test$age))
comp = combn(group,2)
comp
## [,1] [,2] [,3]
## [1,] "<35" "<35" ">60"
## [2,] ">60" "35-60" "35-60"
# 箱线图
p1 <- ggboxplot(BRCA_FPKM_drug, x = "age", y = "Vinblastine_1004",
color = "age", palette =c("#00AFBB", "#E7B800", "#FC4E07"), # 轮廓颜色
bxp.errorbar = T, # 显示误差条
bxp.errorbar.width = 0.5, # 误差条大小
size = 1, # 箱型图边线的粗细
font.x = 15, font.y = 13, # x、y轴标签的大小
font.xtickslab = 13, font.ytickslab = 13, #x、y轴坐标的大小
# ggtheme = theme_bw(), # 修改主题
legend = "right" # 图例放右边
)
my_comparisons = list()
for(i in 1:ncol(comp)){
my_comparisons[[i]] <- comp[,i]
}
p1 + stat_compare_means(comparisons = my_comparisons) +
stat_compare_means(label.x = 0.7, label.y = 4) # Add global p-value
# 选前6个药物绘图
library(ggsci)
result %>% select(1:7) %>%
inner_join(phe_test, "submitter_id.samples") %>%
pivot_longer(2:7, names_to = "drugs", values_to = "ic50") %>%
ggplot(., aes(age, ic50)) +
geom_boxplot(aes(fill = age),
notch = TRUE, notchwidth = 0.9 # 设置凹槽及凹槽宽度
) +
# scale_fill_manual(values = c("<35" = "#00AFBB", "35-60" = "#E7B800", ">60" = "#FC4E07")) + # 自定义配色
scale_fill_jama() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.title.x = element_blank(),
legend.position = "none")+
facet_wrap(vars(drugs), scales = "free_y", nrow = 2)+
stat_compare_means()
