Posted on

药物预测需要训练集,一般来说推荐使用权威资源作为训练集建好模型,这样就可以去预测你自己的数据。

权威的药物预测训练集资源

那么,比较权威的资源一般就是Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC)

Cancer Therapeutics Response Portal (CTRP)

目前主要是CTRP v2,官网是:http://portals.broadinstitute.org/ctrp.v2.1/

  • 481 compounds X 860 CCLs
  • correlations to copy-number and gene-expression data
  • mutation data integrate CCLE and Sanger/MGH calls
  • correlation and enrichment analysis on-the-fly
  • box-whisker visualization in addition to enrichment heatmaps
  • drill-down to scatter plots and concentration-response curves
  • flter by lineage/subtype, growth mode

Genomics of Drug Sensitivity in Cancer (GDSC)

官网是:https://www.cancerrxgene.org/

如果是v2的版本,有809 Cell lines 以及 198 Compounds

如果是看v1版本,987 Cell lines 和 367 Compounds

资源都被整理好了

我们这里直接使用R包oncoPredict整理好的这两个数据库的rdata文件,下载链接是:https://osf.io/c6tfx/

oncoPredict
Contributors: Danielle Maeser Robert Gruener
Date created: 2021-03-26 01:39 PM | Last Updated: 2021-08-15 10:44 PM

下载约700M,重要的文件 如下所示;

 171M Aug 14 17:10 CTRP2_Expr (RPKM, log2(x+1) Transformed).rds
  177M Apr  3  2021 CTRP2_Expr (TPM, log2(x+1) Transformed).rds
  1.1M Apr  3  2021 CTRP2_Res.rds
  119M Apr  3  2021 GDSC1_Expr (RMA Normalized and Log Transformed).rds
  2.0M Apr  3  2021 GDSC1_Res.rds
  100M Apr  3  2021 GDSC2_Expr (RMA Normalized and Log Transformed).rds
  906K Apr  3  2021 GDSC2_Res.rds

可以看到 Cancer Therapeutics Response Portal (CTRP) 数据库里面的细胞系表达量矩阵是来自于转录组测序, 所以 提供了 FPKM和TPM两个版本供用户选择。

然后呢 Genomics of Drug Sensitivity in Cancer (GDSC) 数据库里面的细胞系表达量矩阵应该是芯片,因为它使用了  RMA Normalized and Log Transformed ,标准的芯片数据处理方法。

代码探索 (GDSC) 数据库

直接看 v2的版本,有809 Cell lines 以及 198 Compounds

主要是八百多个细胞系的约2万个基因的表达量矩阵,以及对应八百多细胞系的约200个药物的IC50值。

library(reshape2)
library(ggpubr)
th=theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
dir='./DataFiles/Training Data/'
GDSC2_Expr = readRDS(file=file.path(dir,'GDSC2_Expr (RMA Normalized and Log Transformed).rds'))
dim(GDSC2_Expr)  
GDSC2_Expr[1:4, 1:4]
boxplot(GDSC2_Expr[,1:4])
df=melt(GDSC2_Expr[,1:4])
head(df)
p1=ggboxplot(df, "Var2", "value") +th

# Read GDSC2 response data. rownames() are samples, colnames() are drugs. 
dir
GDSC2_Res = readRDS(file = file.path(dir,"GDSC2_Res.rds"))
dim(GDSC2_Res)  # 805 198
GDSC2_Res[1:4, 1:4]
p2=ggboxplot(melt(GDSC2_Res[ , 1:4]), "Var2", "value") +th ; p2
# IMPORTANT note: here I do e^IC50 since the IC50s are actual ln values/log transformed already, and the calcPhenotype function Paul #has will do a power transformation (I assumed it would be better to not have both transformations)
GDSC2_Res <- exp(GDSC2_Res)  
p3=ggboxplot(melt(GDSC2_Res[ , 1:4]), "Var2", "value") +th ; p3

library(patchwork)
p1+p2+p3

如下所示 :

图片

表达量矩阵被归一化很好,就是(RMA Normalized and Log Transformed),跟绝大部分芯片数据分析一样的,介于 0到15之间。

药物的ic50值,最开始的rds文件里面,也就是说从 (GDSC) 数据库下载得到的是被log转换的,所以又重新使用幂函数转回来。

其中半抑制浓度,或称半抑制率,即IC50,其在间接竞争ELISA标准曲线中是一个非常重要的数据。这样就出现了好几个专有名词啦,不过我们可以简单一点理解,就是能杀死癌细胞的药物浓度的一半。

  • 比如设计 0.01,0.1,1,10,100 这样的药物梯度
  • 发现从浓度1开始可以杀死癌细胞啦,所以我们就认为IC50是  0.5(一般来说不太可能出现浓度更高反而杀不死癌细胞的反常情况)
  • 可以看到,浓度 梯度设计决定了 IC50的分辨率,就是一个区间值就足够啦。
  • 如果IC50超级高,比如解决好几百大几千说明这是一个废物药物,约等于没有疗效。
  • 如果IC50超级低,比如无限接近于0,说明这个药物就是传说中的灵丹妙药!

前面的箱线图,我们展现的是某个药物的八百多细胞系的IC50,这样可以看得出来有一些药物在很多癌症细胞系的表现就是废材,比如 “Cisplatin_1005”    和 “Cytarabine_1006″,当然了,因为我仅仅是展现了4个药物,所以说它们是废材仅仅是相当于 “Camptothecin_1003” 和”Vinblastine_1004″来说。

还有另外一个展现方式,就是看针对具体的细胞系来说,那些药物有奇效那些药物是打酱油。代码如下所示:

   ggboxplot(melt(GDSC2_Res[ 1:4 ,]), "Var1", "value") +th 
图片

因为每个细胞系的箱线图里面都是约200个药物,所以这样的可视化看不出来具体 的药物表现,并没有太大的意义。我们应该是直接看top药物即可:

round(apply(GDSC2_Res[ 1:4 ,], 1, function(x){
  return(c(
    head(sort(x)),
    tail(sort(x))
  ))
}),2)

     COSMIC_906826 COSMIC_687983 COSMIC_910927 COSMIC_1240138
 [1,]          0.00          0.00          0.00           0.05
 [2,]          0.00          0.00          0.00           0.07
 [3,]          0.01          0.01          0.01           0.09
 [4,]          0.04          0.01          0.01           0.21
 [5,]          0.05          0.01          0.01           0.95
 [6,]          0.05          0.01          0.01           0.98
 [7,]       2174.67        310.59        286.42         922.01
 [8,]       2285.10        405.44        388.50         925.43
 [9,]       2859.17        413.01        471.66         939.02
[10,]       3736.69        436.98        489.76         989.15
[11,]       5118.44        626.42        623.64        1105.89
[12,]      15431.05        973.87        803.89        1457.11

可以看到, 每个细胞系都是有自己的特异性药物和废物药物,IC50接近于0的就是神药,那些大几千的就是辣鸡药物啦。

但是,我们可能是更想看到的是药物名字啦!

apply(GDSC2_Res[ 1:4 ,], 1, function(x){ 
  names(x)=gsub('_[0-9]*','',colnames(GDSC2_Res))
  return(c(
    names(head(sort(x))),
    names(tail(sort(x)))
  ))
})

    COSMIC_906826          COSMIC_687983          COSMIC_910927          COSMIC_1240138
 [1,] "Sepantronium bromide" "Daporinad"            "Dactinomycin"         "Luminespib"  
 [2,] "Bortezomib"           "Sepantronium bromide" "Sepantronium bromide" "CDK9"        
 [3,] "Dactinomycin"         "Dactinomycin"         "Bortezomib"           "Dinaciclib"  
 [4,] "Rapamycin"            "Bortezomib"           "Vinblastine"          "Dactinomycin"
 [5,] "Dactolisib"           "Docetaxel"            "Docetaxel"            "CDK9"        
 [6,] "Docetaxel"            "Vinblastine"          "Luminespib"           "Sabutoclax"  
 
 [7,] "KU-55933"             "EPZ5676"              "AZD1208"              "IAP"         
 [8,] "EPZ5676"              "Carmustine"           "Nelarabine"           "Mirin"       
 [9,] "Acetalax"             "Selumetinib"          "Doramapimod"          "PFI3"        
[10,] "Nelarabine"           "Nelarabine"           "SB216763"             "ML323"       
[11,] "Doramapimod"          "Acetalax"             "Carmustine"           "AZD5991"     
[12,] "Temozolomide"         "5-Fluorouracil"       "Temozolomide"         "Carmustine"

前面的6个药物是各自细胞系的神药,后面的6个是废物药物啦。


现在我们可以尝试一下使用R包之oncoPredict对你的表达量矩阵进行药物反应预测啦!

发表oncoPredict这个包的文献非常新:《oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data》,这个通讯作者就是2014年r包pRRophetic同一个人,相当于是炒冷饭吧!

https://ask.qcloudimg.com/http-save/yehe-1154415/a170f5be17dc7b5a57f85aca931ad946.png?imageView2/2/w/1620

使用oncoPredict之前先安装,代码如下:

   install.packages("oncoPredict") 

如果遇到版本问题,请看:https://mp.weixin.qq.com/s/HGfePIQP4yP_nvhjiWdpAQ

使用方法超级简单

首先需要读入训练集的表达量矩阵和药物处理信息,参考前面的教程:药物预测之认识表达量矩阵和药物IC50

rm(list = ls())  ## 魔幻操作,一键清空~
options(stringsAsFactors = F)
library(oncoPredict)
library(data.table)
library(gtools)
library(reshape2)
library(ggpubr)
th=theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
dir='./DataFiles/Training Data/'
GDSC2_Expr = readRDS(file=file.path(dir,'GDSC2_Expr (RMA Normalized and Log Transformed).rds'))
GDSC2_Res = readRDS(file = file.path(dir,"GDSC2_Res.rds"))
GDSC2_Res <- exp(GDSC2_Res) 

这里仍然是以Genomics of Drug Sensitivity in Cancer (GDSC) 的v2作为例子,有了训练集的表达量矩阵和药物处理信息,还需要读入你需要做预测的表达量矩阵。

因为我们这个是教程,所以我就不读取自己的表达量矩阵了,直截了当的从Genomics of Drug Sensitivity in Cancer (GDSC) 的v2里面随机挑选10个细胞系作为要预测的矩阵。

testExpr<- GDSC2_Expr[,sample(1:ncol(GDSC2_Expr),10)]
testExpr[1:4,1:4]  
colnames(testExpr)=paste0('test',colnames(testExpr))
dim(testExpr) 

了训练集的表达量矩阵和药物处理信息,然后也有了待预测的表达量矩阵,接下来就是一个函数的事情啦!这个函数calcPhenotype就是R包 oncoPredict的核心,超级方便!

calcPhenotype(trainingExprData = GDSC2_Expr,
              trainingPtype = GDSC2_Res,
              testExprData = testExpr,
              batchCorrect = 'eb',  #   "eb" for ComBat  
              powerTransformPhenotype = TRUE,
              removeLowVaryingGenes = 0.2,
              minNumSamples = 10, 
              printOutput = TRUE, 
              removeLowVaringGenesFrom = 'rawData' )

这个函数运行取决于你的计算资源,需要半个小时左右。好像也没有多线程的可能性,所以只能是慢慢等了,喝一杯咖啡吧,如果可以的话希望你在咱们《生信技能树》公众号任意教程末尾打赏一杯咖啡也行,我们一起慢慢喝,慢慢等!

从函数运行的log日志来看,本质上就是一个岭回归:

17419  gene identifiers overlap between the supplied expression matrices... 
 
Found2batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data 
 4650 low variabilty genes filtered. 
Fitting Ridge Regression model 
Calculating predicted phenotype... 
Done making prediction for drug 1 of 198 
Fitting Ridge Regression model... 
Calculating predicted phenotype... 
Done making prediction for drug 2 of 198

解读药物预测结果

前面的R包 oncoPredict的核心函数calcPhenotype运行完毕后,会在当前工作目录下面输出 calcPhenotype_Output 文件夹,里面有一个 DrugPredictions.csv的文件,这个都是函数calcPhenotype写死了的。

library(data.table)
testPtype <- fread('./calcPhenotype_Output/DrugPredictions.csv', data.table = F)
testPtype[1:4, 1:4]

不同的数据库资源作为函数的训练集,得到的结果必然是不一样的哦!而且函数也可以调整很多参数。

发表评论

邮箱地址不会被公开。 必填项已用*标注