药物预测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]
不同的数据库资源作为函数的训练集，得到的结果必然是不一样的哦！而且函数也可以调整很多参数。