GDSC(Genomics of Drug Sensitivity in Cancer)數(shù)據(jù)庫由英國桑格研究院開發(fā)，收集腫瘤細(xì)胞對(duì)藥物的敏感度和反應(yīng)。是目前腫瘤細(xì)胞藥物敏感性的最大公共資源，收錄1000多種癌細(xì)胞系中的621個(gè)抗癌化合物的作用廣泛，涉及24種通路。2023年10月更新到了8.5 version。
網(wǎng)址：https://www.cancerrxgene.org/
1. 數(shù)據(jù)集介紹
? GDSC內(nèi)部由和
兩個(gè)數(shù)據(jù)庫組成：
1）主要收錄的是2010-2015年間的測(cè)序和試驗(yàn)結(jié)果。最新版本的GDSC1包含有970個(gè)細(xì)胞系，403個(gè)抗癌化合物，以及333292個(gè)藥物反應(yīng)IC50值和AUC值。（GDSC1 的許多實(shí)驗(yàn)已在 GDSC2 中重復(fù)，如果存在重復(fù)的 IC50，官方建議使用 GDSC2 的結(jié)果）
2）主要收錄的是2015至今的測(cè)序和試驗(yàn)的結(jié)果。最新版本的GDSC2包含有969個(gè)細(xì)胞系，297個(gè)抗癌化合物，以及243466個(gè)藥物反應(yīng)IC50值。
? IC50（half maximal inhibitory concentration），半數(shù)抑制濃度， 指在一定條件下，化合物或藥物能夠抑制生物過程或活性的濃度，使得生物過程或活性被抑制50%。通常，IC50用于評(píng)估藥物的抗生物活性。IC50值越低，說明藥物越強(qiáng)效，因?yàn)樗梢栽诟偷臐舛认乱种颇繕?biāo)生物分子的活性。
? AUC（area under concentration-time curve），血藥濃度-時(shí)間曲線下面積。代表藥物的生物利用度（藥物活性成分從制劑釋放吸收進(jìn)入全身循環(huán)的程度和速度），AUC大則生物利用度高，反之則低。AUC0→∞指藥物從零時(shí)間至所有原形藥物全部消除這一段時(shí)間的藥時(shí)曲線下總面積，反映藥物進(jìn)入血循環(huán)的總量?？偟膩碚f就是較低的AUC值表明細(xì)胞對(duì)治療的敏感性增加。
? EC50（half maximal effective concentration），半數(shù)有效濃度，指在給定實(shí)驗(yàn)條件下，藥物或化合物能夠產(chǎn)生特定效應(yīng)（例如細(xì)胞激活、酶活性增加等）50%的濃度。通常，EC50用于評(píng)估藥物對(duì)疾病相關(guān)分子或細(xì)胞的激活或促進(jìn)作用。較低的EC50值表示藥物對(duì)目標(biāo)生物過程的激活效果更好。
★ IC50和EC50的主要區(qū)別在于它們描述的效應(yīng)類型不同。IC50描述了抑制效應(yīng)的濃度，而EC50描述了激活效應(yīng)的濃度。
2. 數(shù)據(jù)下載
若想研究基因表達(dá)與藥物敏感性的相關(guān)性，需要獲取基因在不同細(xì)胞系中的表達(dá)譜和每個(gè)細(xì)胞系中各種藥物的IC50或者AUC值。
1）表達(dá)譜信息：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")
第一列和第二列為基因名，后續(xù)為各個(gè)細(xì)胞系中對(duì)應(yīng)的基因表達(dá)譜數(shù)據(jù)
2）細(xì)胞系信息：Cell-line-annotation
3）藥物信息：Compounds-annotation
4）每種藥物在每種細(xì)胞中的IC50/AUC值：GDSC1-dataset；GDSC2-dataset
3. oncoPredict包分析
1）載入訓(xùn)練數(shù)據(jù)：
oncoPredict包安裝后需要配套的訓(xùn)練數(shù)據(jù)：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"
# 讀取訓(xùn)練集表達(dá)數(shù)據(jù)
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
# 讀取訓(xùn)練集藥物反應(yīng)數(shù)據(jù)
GDSC2_drug = readRDS(file = file.path(dir,"GDSC2_Res.rds"))
GDSC2_drug <- exp(GDSC2_drug) #下載的數(shù)據(jù)是被log轉(zhuǎn)換過的，逆轉(zhuǎn)回去
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 （二者的樣本名稱對(duì)應(yīng)）
? GDSC2_exp是標(biāo)準(zhǔn)的表達(dá)矩陣格式，行是基因，列是細(xì)胞系，一共17419個(gè)基因，805個(gè)細(xì)胞系
? GDSC2_drug是每個(gè)細(xì)胞系對(duì)每個(gè)藥物的IC50值，行是細(xì)胞系，列是藥物，一共805種細(xì)胞系，198個(gè)藥物
2）加載自己待預(yù)測(cè)的表達(dá)數(shù)據(jù)
此處數(shù)據(jù)通過UCSC Xena（http://xenabrowser.net）網(wǎng)站進(jìn)行下載
# 讀取臨床信息
phe = fread("/data/shumin/GDSC/BRCA/TCGA-BRCA.GDC_phenotype.tsv.gz", header = T, sep = '\t',data.table = F)
# 讀取表達(dá)數(shù)據(jù)
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)
# 提取前兩列，進(jìn)行轉(zhuǎn)換
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函數(shù)對(duì)gene列中的相同基因進(jìn)行合并（一個(gè)探針可能對(duì)應(yīng)多個(gè)gene）
BRCA.fkpm.pro = distinct(BRCA.fkpm.pro,gene,.keep_all = T)
dim(BRCA.fkpm.pro)
## [1] 58387 1219
#把表達(dá)矩陣中g(shù)ene列變成行名
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）預(yù)測(cè)藥物反應(yīng)
calcPhenotype(trainingExprData = GDSC2_exp,
trainingPtype = GDSC2_drug,
testExprData = as.matrix(BRCA.fkpm.pro), # 需要matrix
batchCorrect = 'eb', # 校正參數(shù)，"eb" for array, "standardize" for rnaseq
#IC50是對(duì)數(shù)轉(zhuǎn)換的，所以表達(dá)矩陣也用對(duì)數(shù)轉(zhuǎn)換過的
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
# 取臨床數(shù)據(jù)一列，分組示例
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
# 合并數(shù)據(jù)
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, # 箱型圖邊線的粗細(xì)
font.x = 15, font.y = 13, # x、y軸標(biāo)簽的大小
font.xtickslab = 13, font.ytickslab = 13, #x、y軸坐標(biāo)的大小
# 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個(gè)藥物繪圖
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 # 設(shè)置凹槽及凹槽寬度
) +
# 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()
