药物敏感性分析之pRRophetic
药物敏感性分析是生信数据挖掘常用的技能之一，目前做药敏分析最常见的就是两个R包：pRRophetic和oncoPredict。
这两个包的作者都是同一个人，oncoPredict可以看做是pRRophetic的升级版。两个R包的使用基本上是一样的思路，只不过使用的训练数据集不同而已。
在介绍R包的使用之前，需要大家先了解一下常用的药物敏感性数据库，最好是去到这些数据库的主页看看或者读一读相关的文献，对这些数据库有一个大致的了解。
- 常用药敏数据库
- pRRophetic方法学介绍
- 安装
- 预测不同组别患者对化疗药物的敏感性
- 其他示例
- pRRopheticCV
- 使用CCLE示例数据
- 自定义训练集
- 自带数据探索
- 预测全的药物的敏感性
- 参考
常用药敏数据库
药敏数据库非常多，但最常用的无非就是GDSC/CTRP/CCLE等，在珠江肿瘤公众号中早就介绍过这些数据库了，所以我就不重复了，大家可以参考以下链接。
以下链接介绍了GDSC、CTRP、CCLE、NCI-60、DepMap、Pharmacodb等数据库，是非常棒的参考资料：
pRRophetic方法学介绍
这个R包的思路其实很简单，就是根据已知的细胞系表达矩阵和药物敏感性信息作为训练集建立模型，然后对新的表达矩阵进行预测。已知的信息就是从直接从上面介绍的数据库下载的，pRRophetic包使用的是CGP和CCLE的数据，但是CCLE的药敏数据只有24种药物和500多个细胞系，数据量比较少，所以通常大家使用的都是CGP的数据。
作者专门发了一篇文章，详细介绍该包背后的方法和原理：Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines。
作者把上面这篇文献中的方法变成了pRRophetic包，又发了一篇文章，非常妙：pRRophetic: An R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels
其中有一个简化版的方法学介绍，我截图如下，如果想要详细了解，建议阅读原文献哦：
简单来说，使用的表达矩阵是芯片数据，训练集和测试集分别进行quantile-normalization, 去除方差低的基因，用每个基因作为预测变量，药物的IC50作为结果变量拟合岭回归，然后使用拟合好的模型对新的数据进行预测。
安装
这个R包非常古老，虽然文章里说会不断更新，但是很明显没有更新过。
需要首先安装依赖包，然后通过以下链接下载pRRophetic_0.5.tar.gz这个压缩包，进行本地安装。
链接：https://osf.io/dwzce/?action=download
#安装依赖包
BiocManager::install(c("car", "ridge", "preprocessCore", "genefilter", "sva"))
#本地安装
install.packages("E:/R/R包/pRRophetic_0.5.tar.gz", repos = NULL, type = "source")
这个包太老了，有些版本比较新的R可能安装不了，我使用的R4.2.0安装没有任何问题。
但是安装之后还是不能使用，因为它太古老了，可能会遇到以下报错：
# 报错：
Error in if (class(exprMatUnique) == "numeric") { :
the condition has length > 1
# 或者
Error in if (class(testExprData) != "matrix") stop("ERROR: \"testExprData\" must be a matrix.") : the condition has length > 1
遇到了不要惊慌，毕竟果子老师已经帮我们解决这个问题，按照果子老师的介绍重新安装即可：
基因表达量预测药物反应的R包pRRophetic近期报错解决方案
预测不同组别患者对化疗药物的敏感性
在包的github中作者给了一个使用示例：https://github.com/paulgeeleher/pRRophetic/blob/master/vignetteOutline.pdf
下面我们结合这个示例简单介绍下这个包的使用。
通常我们会根据某种方法把所有样本分为不同的组（比如最常见的高风险组/低风险组，或者不同的分子亚型等），然后想看看不同的组对某种药物的敏感性。
这个包就可以帮你做这样的事情，而且只需要你提供自己的表达矩阵即可，它默认会使用cgp2014的数据作为训练集建立模型，然后对你的表达矩阵进行预测，这样你就可以得到每个样本的IC50值。
除了cgp2014，你还可以选择cgp2016作为训练数据，cgp2016有251种药物，cgp2014只有138种。
前面介绍过的GDSC（Genomics of Drug Sensitivity in Cancer），是CGP项目（Cancer Genome Project）的一部分。CGP的官网：https://www.cancerrxgene.org/。
加载R包：
library(pRRophetic)
在预测对某个药物的敏感性前，最好先评估数据的正态性，因为CGP中的许多药物的IC50并不是呈正态分布的，此时是不适合使用线性模型的。
用R包自带的硼替佐米数据做个演示，先看下硼替佐米这个药的IC50是不是符合正态分布：
data("bortezomibData")
pRRopheticQQplot("Bortezomib")
从这个QQ图来看其实不是非常符合，但还算可以，我们就认为它符合吧。
然后就可以用pRRopheticPredict预测对这个药物的敏感性了，这也是这个包最重要的函数。
我们这里用的是示例表达矩阵，你用的时候只需要换成自己的表达矩阵即可。
exprDataBortezomib是一个标准的表达矩阵，行是基因，列是样本：
dim(exprDataBortezomib) #22283个基因，264个样本
## [1] 22283 264
exprDataBortezomib[1:4,1:4]
## GSM246523 GSM246524 GSM246525 GSM246526
## <NA> 235.5230 498.2220 309.2070 307.5690
## RFC2 41.4470 69.0219 69.3994 36.9310
## HSPA6 84.8689 56.8352 49.4388 54.6669
## PAX8 530.4490 457.9310 536.1780 325.3630
预测：
predictedPtype <- pRRopheticPredict(testMatrix = exprDataBortezomib, #表达矩阵
drug = "Bortezomib", # 药物
tissueType = "blood",
batchCorrect = "eb", #训练集和测试集数据整合方法，默认eb，即使用combat
powerTransformPhenotype = T, # 是否进行幂转换
selection=1, # 遇到名字重复的基因取平均值
dataset = "cgp2014")
##
## 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
tissueType指定你想用CGP细胞系中的哪些类型肿瘤作为训练集，默认是all；
结果是一个命名向量，就是每个样本的IC50值：
head(predictedPtype)
## GSM246523 GSM246524 GSM246525 GSM246526 GSM246527 GSM246528
## -6.808324 -5.557028 -5.382334 -3.999054 -6.330220 -4.751816
这个示例数据中所有的样本可以被分为两组，一组是NR组，另一组是R组，通常你的表达矩阵也会分组的，比如根据某个方法分成高风险组和低风险组，一样的意思。
我们就可以对两组的IC50做个t检验：
t.test(predictedPtype[((studyResponse == "PGx_Responder = NR") & bortIndex)],
predictedPtype[((studyResponse == "PGx_Responder = R") & bortIndex)],
alternative="greater")
##
## Welch Two Sample t-test
##
## data: predictedPtype[((studyResponse == "PGx_Responder = NR") & bortIndex)] and predictedPtype[((studyResponse == "PGx_Responder = R") & bortIndex)]
## t = 4.1204, df = 165.24, p-value = 2.984e-05
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 0.3975589 Inf
## sample estimates:
## mean of x mean of y
## -4.372173 -5.036370
还可以画个图展示：
library(ggplot2)
library(ggpubr)
df <- stack(list(NR=predictedPtype[((studyResponse == "PGx_Responder = NR") & bortIndex)],
R=predictedPtype[((studyResponse == "PGx_Responder = R") & bortIndex)]))
ggboxplot(df, x="ind",y="values",fill="ind",alpha=0.3,palette = "lancet",
ylab="Predicted Bortezomib Sensitivity",
xlab="Clinical Response"
)+
stat_compare_means(method = "t.test")+
theme(legend.position = "none")
这张图就是文献中最常见的图了。
下面再展示下预测对厄洛替尼的敏感性，这个药物的IC50明显不符合正态分布：
pRRopheticQQplot("Erlotinib")
所以此时我们使用pRRopheticLogisticPredict函数预测样本的IC50值：
predictedPtype_erlotinib <- pRRopheticLogisticPredict(exprDataBortezomib,
"Erlotinib",
selection=1)
##
## 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
## Fitting model, may take some time....
后面的分析就都是一样的了~
其他示例
pRRopheticCV
为了说明这个包的预测结果的准确性，还可以使用pRRopheticCV函数查看预测结果和真实结果的一致性，使用5折交叉验证：
cvOut <- pRRopheticCV("Bortezomib", cvFold=5, testExprData=exprDataBortezomib)
##
## 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
##
## 1 of 5 iterations complete.
## 2 of 5 iterations complete.
## 3 of 5 iterations complete.
## 4 of 5 iterations complete.
## 5 of 5 iterations complete.
查看结果：
summary(cvOut)
##
## Summary of cross-validation results:
##
## Pearsons correlation: 0.44 , P = 6.37210297840637e-15
## R-squared value: 0.2
## Estimated 95% confidence intervals: -4.21, 3.36
## Mean prediction error: 1.61
真实结果和预测结果的相关性只有0.42，还给出了P值、R^2、预测错误率等信息，可以画个图展示下真实结果和预测结果：
plot(cvOut)
使用CCLE示例数据
CCLE中只有24个药物，500+细胞系，用的很少，数据量比CGP少太多了。该包自带了一个CCLE数据ccleData，其使用方法和上面完全一样，就不重复介绍了。
自定义训练集
指定训练用的表达矩阵和对应的样本类别，再提供一个表达矩阵，就可以预测该表达矩阵每个样本对药物的敏感性。也就是说这个方法可以让你能够使用自己的训练数据~但是我好像并没有见到这么做的，如果大家有见过的，欢迎告诉我~
下面我们继续用硼替佐米数据作为示例进行演示。
我们先从exprDataBortezomib这个完整的表达矩阵提取一部分数据作为训练用的表达矩阵，并且也提取这部分样本的类别（有5个类别：CR、PR、MR、NC、PD）。
然后再提取一部分表达矩阵作为测试用表达矩阵，来预测这部分样本对硼替佐米的敏感性。
准备训练数据和测试数据：
# 训练用表达矩阵
trainExpr <- exprDataBortezomib[, (detailedResponse %in% c(1,2,3,4,5)) & studyIndex %in% c("studyCode = 25", "studyCode = 40")]
# 训练用样本类型
trainPtype <- detailedResponse[(detailedResponse %in% c(1,2,3,4,5)) & studyIndex %in% c("studyCode = 25", "studyCode = 40")]
# 预测用表达矩阵
testExpr <- exprDataBortezomib[, (detailedResponse %in% c(1,2,3,4,5)) & studyIndex %in% c("studyCode = 39")]
dim(testExpr) # 141个样本
## [1] 22283 141
下面就可以预测样本对药物的敏感性了：
ptypeOut <- calcPhenotype(trainExpr, trainPtype, testExpr, selection=1)
##
## 22283 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
##
## 2500 low variabilty genes filtered.
## Fitting Ridge Regression model... Done
##
## Calculating predicted phenotype...Done
这个结果就是141个样本的预测敏感性：
length(ptypeOut)
## [1] 141
head(ptypeOut)
## GSM246530 GSM246536 GSM246537 GSM246539 GSM246540 GSM246544
## 2.990402 2.615408 3.314234 2.718105 2.578793 2.823383
有了这个预测的结果，我们可以与真实的结果做一个相关性分析：
# 提取真实结果
testPtype <- detailedResponse[(detailedResponse %in% c(1,2,3,4,5)) & studyIndex %in% c("studyCode = 39")]
# 相关性分析
cor.test(testPtype, ptypeOut, alternative="greater")
##
## Pearson's product-moment correlation
##
## data: testPtype and ptypeOut
## t = 2.1512, df = 139, p-value = 0.01659
## alternative hypothesis: true correlation is greater than 0
## 95 percent confidence interval:
## 0.04142448 1.00000000
## sample estimates:
## cor
## 0.1795014
还可以做t检验：
t.test(ptypeOut[testPtype %in% c(3,4,5)], ptypeOut[testPtype %in% c(1,2)],
alternative="greater")
##
## Welch Two Sample t-test
##
## data: ptypeOut[testPtype %in% c(3, 4, 5)] and ptypeOut[testPtype %in% c(1, 2)]
## t = 2.0182, df = 137.43, p-value = 0.02276
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 0.02533599 Inf
## sample estimates:
## mean of x mean of y
## 2.646449 2.505269
自带数据探索
这个包自带的所有数据可以在包的安装目录中查看，就是这几个：
rm(list = ls())
library(pRRophetic)
加载数据查看一下：
data(cgp2016ExprRma)
dim(cgp2016ExprRma)
## [1] 17419 1018
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
这个是cgp2016版本的表达矩阵，其中行是基因，列是细胞系，一共17419个基因，1018个细胞系。
data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016)
length(unique(drugData2016$Drug.name))
## [1] 251
head(unique(drugData2016$Drug.name),30)
## [1] "Erlotinib" "Rapamycin" "Sunitinib"
## [4] "PHA-665752" "MG-132" "Paclitaxel"
## [7] "Cyclopamine" "AZ628" "Sorafenib"
## [10] "VX-680" "Imatinib" "TAE684"
## [13] "Crizotinib" "Saracatinib" "S-Trityl-L-cysteine"
## [16] "Z-LLNle-CHO" "Dasatinib" "GNF-2"
## [19] "CGP-60474" "CGP-082996" "A-770041"
## [22] "WH-4-023" "WZ-1-84" "BI-2536"
## [25] "BMS-536924" "BMS-509744" "CMK"
## [28] "Pyrimethamine" "JW-7-52-1" "A-443654"
上面这个是cgp2016版本的细胞系和药物敏感性信息，包含了每种细胞系对每种药物的IC50等信息，可以看到其中一共有251种药物，cgp2014只有138种药物（可以通过?pRRopheticPredict查看帮助文档确定）。
data(drugAndPhenoCgp)
这里面是一些原始文件，貌似平常用不到，大家感兴趣可以自己探索下。
可以看到其中还有一个ccleData，其实和上面用到的硼替佐米数据是一样的，只不过一个来自于CGP，另一个来自于CCLE而已，就不展示了。
预测全部药物的敏感性
假如我们要对自己的表达矩阵预测所有药物的敏感性，只需要把所有的药物提取出来，写个循环即可，这里以cgp2016的药物为例。
以下这段代码来自生信技能树：药物预测R包之pRRophetic
耗时巨长！！
library(parallel)
library(pRRophetic)
# 加载cgp2016的药敏信息
data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016) # drugData2016
data(cgp2016ExprRma) # cgp2016ExprRma
data(bortezomibData)
#提取cgp2016的所有药物名字
possibleDrugs2016 <- unique( drugData2016$Drug.name)
#possibleDrugs2016
# 用system.time来返回计算所需时间
#head(possibleDrugs2016)
system.time({
cl <- makeCluster(8)
results <- parLapply(cl,possibleDrugs2016[1:10],#只用前10个测试，用全部时间太长
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) # 关闭集群
})
画个图看看，画图之前需要一些数据格式的转换，就是常规的长宽转换，加名字即可。
然后使用ggplot系列包画图、添加显著性，一气呵成，非常简单，所以R语言基础是非常有必要的。
library(tidyr)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(ggsci)
plot_df <- res.df %>%
as.data.frame() %>%
t() %>%
bind_cols(studyResponse) %>%
bind_cols(bortIndex) %>%
filter(!studyResponse == "PGx_Responder = IE", bortIndex == TRUE)
names(plot_df) <- c(possibleDrugs2016[1:10],"studyResponse","bortIndex")
plot_df %>%
pivot_longer(1:10,names_to = "drugs",values_to = "ic50") %>%
ggplot(., aes(studyResponse,ic50))+
geom_boxplot(aes(fill=studyResponse))+
scale_fill_jama()+
theme(axis.text.x = element_text(angle = 45,hjust = 1),
axis.title.x = element_blank(),
legend.position = "none")+
facet_wrap(vars(drugs),scales = "free_y",nrow = 2)+
stat_compare_means()
easy！这次内容挺多的，下次再介绍oncoPredict吧。
参考
生信技能树：药物预测R包之pRRophetic