有了前面的教程：藥物預測之認識表達量矩陣和藥物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的一致性！
