有了前面的教程：药物预测之认识表达量矩阵和药物IC50 的背景知识铺垫，认识了Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 两个数据库资源。也介绍了2021年7月新鲜出炉的 药物预测R包之oncoPredict
还可以尝试一下同一个团队早在2014年就出品的R包之 pRRophetic ，也可以对你的表达量矩阵进行药物反应预测啦！很有意思的是这个包虽然是2014就发表了，文章是：《pRRophetic: An R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels》
但是它一直没有把包提交到cran或者bioconductor，而是自己维护：
文章短小精悍，就5个参考文献，而且它并不使用我们介绍的Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 两个数据库资源，而是使用 Cancer Genome Project (CGP)计划里面的表达量矩阵和药物处理信息，该数据库有138 anticancer drugs against 727 cell lines
使用 pRRophetic 之前先安装，自己下载pRRophetic_0.5.tar.gz这个压缩包（大于500M哦），然后当前工作目录下面使用代码如下：
# 自己下载pRRophetic_0.5.tar.gz这个压缩包（大于500M哦）
# website (http://genemed.uchicago.edu/~pgeeleher/pRRophetic)
# GitHub (https://github.com/paulgeeleher/pRRophetic).
BiocManager::install(c('sva', 'car', 'genefilter', 'preprocessCore', 'ridge'))
# 假如提前你缺啥，就使用上面的代码安装缺失的包
install.packages("pRRophetic_0.5.tar.gz", repos = NULL, dependencies = TRUE)
安装成功之后，在你的R包pRRophetic的文件夹路径下面有如下所示文件 ：
77K Oct 4 10:46 Cell_Lines_Details-1.csv
292M Oct 4 10:46 Cell_line_RMA_proc_basalExp.txt
23M Oct 4 10:46 PANCANCER_IC_Tue Aug 9 15_28_57 2016.csv
6.4M Oct 4 10:46 PANCANCER_IC_Tue_Aug_9_15_28_57_2016.RData
25M Oct 4 10:46 bortezomibData.RData
122M Oct 4 10:46 ccleData.RData
126M Oct 4 10:46 cgp2016ExprRma.RData
205B Oct 4 10:46 datalist
72M Oct 4 10:46 drugAndPhenoCgp.RData
可以看到虽然这个包是2014的，但是里面有一些数据是2016的，应该是作者还在维护。不过它毕竟是没有提交到cran或者bioconductor，还是有点担心啊！
Bortezomib (PS-341) 是一种可逆性和选择性的蛋白酶体 (proteasome) 抑制剂，通过靶向苏氨酸残基有效抑制 20S 蛋白酶体 (Ki=0.6 nM)。Bortezomib 破坏细胞周期、诱导细胞凋亡以及抑制核因子 NF-κB。
简而言之：Bortezomib 是第一种蛋白酶体抑制剂，具有抗癌活性。
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，应该是有无响应的简单情况。
library(pRRophetic)
# exprMatCcle (ccleData)
data(ccleData)
dim(exprMatCcle)
exprMatCcle[1:4,1:4]
boxplot(exprMatCcle[,1:4])
# 1037个细胞系的18988个基因组成的 表达量矩阵
head(sensDataCcle)
tmp=as.data.frame(table(sensDataCcle[,2:3]) )
length(unique(sensDataCcle$Compound))
length(unique(sensDataCcle$Primary.Cell.Line.Name))
可以看到ccle数据库里面的表达量矩阵信息很丰富，有1037个细胞系的18988个基因表达。
但是药物处理信息只有24个药物，五百多个细胞系，如下所示：
> tmp=as.data.frame(table(sensDataCcle[,2:3]) )
> length(unique(sensDataCcle$Compound))
[1] 24
> length(unique(sensDataCcle$Primary.Cell.Line.Name))
[1] 504
> head(tmp)
Primary.Cell.Line.Name Compound Freq
1 1321N1 17-AAG 1
2 22Rv1 17-AAG 1
3 42-MG-BA 17-AAG 1
4 5637 17-AAG 1
5 639-V 17-AAG 1
6 697 17-AAG 1
> dim(tmp)
[1] 12096 3
全部的24种药物是：
> unique(sensDataCcle$Compound)
[1] "17-AAG" "AEW541" "AZD0530"
[4] "AZD6244" "Erlotinib" "Irinotecan"
[7] "L-685458" "LBW242" "Lapatinib"
[10] "Nilotinib" "Nutlin-3" "PD-0325901"
[13] "PD-0332991" "PF2341066" "PHA-665752"
[16] "PLX4720" "Paclitaxel" "Panobinostat"
[19] "RAF265" "Sorafenib" "TAE684"
[22] "TKI258" "Topotecan" "ZD-6474"
可以看到 Cancer Genome Project (CGP)里面的也是一千多个细胞系的2万个基因的表达量矩阵，关键是它药物超过了200种。
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)
这个核心函数pRRopheticPredict，有非常多的参数，不过绝大部分都是有默认值，如下所示：
pRRopheticPredict(testMatrix, drug, tissueType = "all", batchCorrect = "eb",
powerTransformPhenotype = TRUE, removeLowVaryingGenes = 0.2,
minNumSamples = 10, selection = -1, printOutput = TRUE,
removeLowVaringGenesFrom = "homogenizeData", dataset = "cgp2014")
我们首先使用 默认的 cgp2014数据集，来预测 前面的 硼替佐米（Bortezomib）药物信息，看看是否一致。
library(pRRophetic)
library(ggplot2)
data(bortezomibData)
dim(exprDataBortezomib)
exprDataBortezomib[1:4,1:4]
boxplot(exprDataBortezomib[,1:4])
# 264个样品的表达量矩阵
table(studyResponse)
predictedPtype <- pRRopheticPredict(testMatrix=exprDataBortezomib,
drug="Bortezomib",
tissueType = "all",
batchCorrect = "eb",
selection=1,
dataset = "cgp2014")
如果是实际使用，我们应该是读入自己的表达量哦，不过我们这里是演示如何使用它，所以直接用exprDataBortezomib矩阵来测试pRRophetic包的核心函数pRRopheticPredict即可。
因为就预测一个药物，所以速度很快：
11683 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
2324 low variabilty genes filtered.
Fitting Ridge Regression model... Done
Calculating predicted phenotype...Done
让我们看看，预测情况和时间情况的差异：
boxplot(predictedPtype)
df <- stack(list(NR = predictedPtype[((studyResponse == 'PGx_Responder = NR') & bortIndex)],
R = predictedPtype[((studyResponse == 'PGx_Responder = R') & bortIndex)]))
head(df)
ggplot(data = df,
aes(y = values,
x = ind))+
geom_boxplot(alpha = 0.3,
fill = c('#e94753','#47a4e9'))+
theme_bw()+
ylab('Predicted Bortezomib Sensitivity') +
xlab('Clinical Response')
如下所示：
总体来说，确实是R的响应组，里面的病人预测得到的药物越敏感！
这个时候因为是每个药物都需要走前面的 pRRophetic包的核心函数pRRopheticPredict，所以可以写循环啦，而且可以加入并行机制。
这个时候我们选择 cgp2016数据集，全部的代码如下所示：
library(parallel)
library(pRRophetic)
library(ggplot2)
data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016)
data(cgp2016ExprRma)
possibleDrugs2016 <- unique( drugData2016$Drug.name)
possibleDrugs2016
# 用system.time来返回计算所需时间
head(possibleDrugs2016)
system.time({
cl <- makeCluster(8)
results <- parLapply(cl,possibleDrugs2016,
function(x){
library(pRRophetic)
data(bortezomibData)
predictedPtype=pRRopheticPredict(
testMatrix=exprDataBortezomib,
drug=x,
tissueType = "all",
batchCorrect = "eb",
selection=1,
dataset = "cgp2016")
return(predictedPtype)
}) # lapply的并行版本
res.df <- do.call('rbind',results) # 整合结果
stopCluster(cl) # 关闭集群
})
这个函数运行取决于你的计算资源，需要半个小时左右。慢慢等吧，喝一杯咖啡吧，如果可以的话希望你在咱们《生信技能树》公众号任意教程末尾打赏一杯咖啡也行，我们一起慢慢喝，慢慢等！
如果是自己的真实表达量矩阵，需要走 pRRophetic包的核心函数pRRopheticPredict进行预测，只需要读入进去，替换前面的 testMatrix=exprDataBortezomib 即可！
不同的数据库资源作为函数的训练集，得到的结果必然是不一样的哦！而且函数也可以调整很多参数。
这个时候其实可以看看 前面的2021年7月新鲜出炉的 药物预测R包之oncoPredict 结果跟本次介绍的药物预测R包之pRRophetic的一致性！
联系客服
