药物预测之认识表达量矩阵和药物IC50

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

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

那么，比较权威的资源一般就是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个是废物药物啦。