药物预测R包之oncoPredict

有了前面的教程：药物预测之认识表达量矩阵和药物IC50 的背景知识铺垫，认识了Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 两个数据库资源。

现在我们可以尝试一下使用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同一个人，相当于是炒冷饭吧！

使用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


差不多是最简单的最常用的机器学习统计学方法了，如果你对其它机器学习方法感兴趣，可以看生信菜鸟团鲍志炜的专栏分享：

• #开新坑#一起来学习机器学习吧~
• 机器学习实战-k-邻近算法
• 机器学习实战-决策树
• 机器学习实战-朴素贝叶斯
• 机器学习实战-逻辑回归
• 机器学习实战-支持向量机（一）
• 机器学习实战-支持向量机·sklearn参数详解
• 机器学习实战-Adaboost
• 机器学习实战-数据预处理
• 机器学习实战-机器学习性能指标
• 机器学习实战-特征选择
• 机器学习实战-分类算法学习目录
• bioBakeryWorkflows-宏基因组分析流程

其实在生物信息学领域，可能仍然是R语言应用更为广泛，简单的一个R包或者函数即可。所以并不是一定要学习python哈。

解读药物预测结果

前面的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]


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