药物预测之相关性为什么不行

生信技能树

发布于 2021-10-12 12:16:08

1.3K0

发布于 2021-10-12 12:16:08

文章被收录于专栏：生信技能树

有了前面的教程：药物预测之认识表达量矩阵和药物IC50 的背景知识铺垫，认识了Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 两个数据库资源。也介绍了2021年7月新鲜出炉的 药物预测R包之oncoPredict的用法以及一个超级老旧的R包的用法，见：药物预测R包之pRRophetic

现在可以尝试一下理解药物预测的原理啦，首先我们看看仅仅是相关性，能不能搞定药物预测！

首先认识需要预测的表达量矩阵

我们还是使用Bortezomib的表达量矩阵，来预测里面的全部样品对Bortezomib这个药物的反应情况！Bortezomib 是第一种蛋白酶体抑制剂，具有抗癌活性。在 pRRophetic 包里面有 bortezomibData数据集，主要是一个表达量矩阵，264个样品，以及这264个样品的结局（药物处理后）。

library(pRRophetic)
data(bortezomibData)
dim(exprDataBortezomib)
exprDataBortezomib[1:4,1:4]
boxplot(exprDataBortezomib[1:4,])
# 264个样品的表达量矩阵
table(studyResponse)
# 264个病人的结局 （药物处理后）
PGx_Responder = IE PGx_Responder = NR  PGx_Responder = R 
                25                126                113

具体是什么癌症什么病人就需要看该文章以及该数据集的来源文献啦，但是药物处理结局事件是很明显的，主要是区分R和NR，应该是有无响应的简单情况。

然后提取Bortezomib这个药物的数据资源

有了上面的已知信息，现在我们来看看Genomics of Drug Sensitivity in Cancer (GDSC) 这个数据库资源。

rm(list = ls())
library(pRRophetic)   
# exprMatCcle (ccleData)       
data(cgp2016ExprRma) 
dim(cgp2016ExprRma) 
data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016)
length(unique(drugData2016$Drug.name))
unique(drugData2016$Drug.name)

如下所示的表达量矩阵和药物信息：

> cgp2016ExprRma[1:4,1:4]
         CAL-120   DMS-114   CAL-51    H2869
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
> drugData2016[1:4,1:4]
  Drug.name Drug.Id Cell.line.name Cosmic.sample.Id
1 Erlotinib       1         MC-CAR           683665
2 Erlotinib       1            ES3           684055
3 Erlotinib       1            ES5           684057
4 Erlotinib       1            ES7           684059

现在我们提前Bortezomib这个药物的反应情况：

Bort_drugs=drugData2016[drugData2016$Drug.name=='Bortezomib',]
dim(Bort_drugs)
kp =  Bort_drugs$Cell.line.name %in% colnames(cgp2016ExprRma) 
table(kp)
Bort_drugs=Bort_drugs[kp,]
pos=match( Bort_drugs$Cell.line.name ,colnames(cgp2016ExprRma)  )
Bort_ExprRma=cgp2016ExprRma[,pos]

可以看到 Genomics of Drug Sensitivity in Cancer (GDSC) 这个数据库里面虽然是有1018个细胞系，但是有Bortezomib药物处理信息的只有389个，有意思的是有38个细胞系有Bortezomib药物处理信息但是没有表达量矩阵。

合并两个表达量矩阵

两个表达量矩阵各自的技术都不一样，所以肯定是存在批次效应的，

ids=intersect(rownames(Bort_ExprRma),rownames(exprDataBortezomib))
m=cor(cbind(exprDataBortezomib[ids,],
          Bort_ExprRma[ids,]))
library(pheatmap)
pheatmap(m)
ac=data.frame(batch=c(
  rep('test',ncol(exprDataBortezomib)),
  rep('GDSC',ncol(Bort_ExprRma))
))
rownames(ac)=colnames(m)
pheatmap(m,show_rownames = F,show_colnames = F,
         annotation_col = ac)

两个表达量矩阵的基因交集居然只有一万多个：

> dim(exprDataBortezomib)
[1] 22283   264
> dim(Bort_ExprRma)
[1] 17419   389
> length(ids)
[1] 11437

各自的表达量矩阵的相关性如下所示：

批次效应

可以看到，妥妥的批次效应，不过，我们仍然是可以尝试给每个细胞系一个相关性最高选项。

ch_m=m[1:ncol(exprDataBortezomib),
       (ncol(exprDataBortezomib)+1):ncol(m)]
pheatmap(ch_m,show_rownames = F,show_colnames = F)
ch_m[1:4,1:4]
colnames(Bort_drugs)
pred_drug = Bort_drugs[apply(ch_m, 1, which.max), "IC50"]
boxplot(pred_drug~studyResponse)

可以看到，基本上没有如何差异，明明是不同样品对Bortezomib药物处理应该是有R和NR的差异！

而且我把这样的简单的相关性预测结果，跟之前的R包预测结果（见：药物预测R包之pRRophetic ）对比发现，也没有一致性。

既然有批次效应，去除后会不会好一点呢？

最简单的当然是sva去除啦！

    combat_edata1 = ComBat(dat=as.matrix(m), 
                            batch=c(
                             rep('test',ncol(exp_test)),
                             rep('GDSC',ncol(exp_db))
                           ), mod=NULL, 
                           par.prior=T, prior.plots=F)
    combat_edata1[1:4,1:4]

有意思的是，sva确实是把两个表达量矩阵的差异效应给去除了，而且预测的药物效果稍微合理了一点，至少在R组（药物敏感）里面的IC50值显著的低，如下所示：

它虽然是跟之前的R包预测结果（见：药物预测R包之pRRophetic ）相关性并不是非常高，但是起码是有统计学显著的正相关，而且也不能说之前的R包预测结果（见：药物预测R包之pRRophetic ）就是金标准！

也仅仅是，我们简简单单的把需要预测的表达量矩阵，去跟数据库里面有药物作用信息的表达量矩阵，进行一个相关性计算（前提是去除了批次效应），然后给需要预测的样品每个都找到了一个最相关的数据库里面的样品，直接把药物作用值赋值给预测样品。居然都可以勉勉强强得到还算是合理的结果啊！

相关性预测药物效果的逻辑链

给需要预测的样品去数据库资源里面找到一个跟它表达量相关性最高的样品，直接把该样品的药物作用值赋值给预测样品，是勉勉强强得到还算是合理的结果。这个推理的逻辑基础就是，药物作用跟表达量是线性相关的。

一个很简单的例子就是，根据颜值来预测收入，我们有一个数据库记录了2万人的颜值和这2万人的收入，发现确实是线性相关，颜值越高收入越高，这里面的相关系数目前不知道，我们也不想去分析。这个时候有一个全新的团队里面的人的收入是未知的，但是颜值是可以判断的，我们对每个人都判断颜值，然后去2万人里面找到最匹配它的颜值，就把匹配到的人收入直接辅助给它作为预测值。

听起来似乎是很合理，但是它的假设前提就是2万人的颜值和这2万人的收入确实是线性相关，比如颜值越高收入越高。但是实际情况是部分颜值低到可怜的人居然收入反转了比如马云，也就是说这个相关性有可能是一个一元二次方程，如果你还记得一点点以前学过的数学，就可以想象得到同样的超高收入会对应两个极端的颜值（马云和李彦宏）。这个时候如果仅仅是看颜值相似性你就没办法判断它的收入了，非常的麻烦。

那，我们来探索一下， Genomics of Drug Sensitivity in Cancer (GDSC) 这个数据库资源里面的各个样品的表达量信息是否跟各个药物处理的IC50值有相关性！代码如下：

d=as.matrix(dist(cor(Bort_ExprRma[,order(Bort_drugs$IC50)])))
df=data.frame(dist=d[1,],
              ic50=sort(Bort_drugs$IC50))
library(ggpubr)
ggscatter(df, x = "dist", y = "ic50",
          color = "black", shape = 21, size = 3, # Points color, shape and size
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE, # Add confidence interval
          cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
          cor.coeff.args = list(method = "pearson",  label.sep = "\n")
)

略微有点尴尬，基本上没有相关性：