有了前面的教程：藥物預(yù)測(cè)之認(rèn)識(shí)表達(dá)量矩陣和藥物IC50 的背景知識(shí)鋪墊，認(rèn)識(shí)了Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 兩個(gè)數(shù)據(jù)庫(kù)資源。也介紹了2021年7月新鮮出爐的 藥物預(yù)測(cè)R包之oncoPredict
還可以嘗試一下同一個(gè)團(tuán)隊(duì)早在2014年就出品的R包之 pRRophetic ，也可以對(duì)你的表達(dá)量矩陣進(jìn)行藥物反應(yīng)預(yù)測(cè)啦！很有意思的是這個(gè)包雖然是2014就發(fā)表了，文章是：《pRRophetic: An R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels》
但是它一直沒(méi)有把包提交到cran或者bioconductor，而是自己維護(hù)：
文章短小精悍，就5個(gè)參考文獻(xiàn)，而且它并不使用我們介紹的Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 兩個(gè)數(shù)據(jù)庫(kù)資源，而是使用 Cancer Genome Project (CGP)計(jì)劃里面的表達(dá)量矩陣和藥物處理信息，該數(shù)據(jù)庫(kù)有138 anticancer drugs against 727 cell lines
使用 pRRophetic 之前先安裝，自己下載pRRophetic_0.5.tar.gz這個(gè)壓縮包（大于500M哦），然后當(dāng)前工作目錄下面使用代碼如下：
# 自己下載pRRophetic_0.5.tar.gz這個(gè)壓縮包（大于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
可以看到雖然這個(gè)包是2014的，但是里面有一些數(shù)據(jù)是2016的，應(yīng)該是作者還在維護(hù)。不過(guò)它畢竟是沒(méi)有提交到cran或者bioconductor，還是有點(diǎn)擔(dān)心啊！
Bortezomib (PS-341) 是一種可逆性和選擇性的蛋白酶體 (proteasome) 抑制劑，通過(guò)靶向蘇氨酸殘基有效抑制 20S 蛋白酶體 (Ki=0.6 nM)。Bortezomib 破壞細(xì)胞周期、誘導(dǎo)細(xì)胞凋亡以及抑制核因子 NF-κB。
簡(jiǎn)而言之：Bortezomib 是第一種蛋白酶體抑制劑，具有抗癌活性。
library(pRRophetic)
data(bortezomibData)
dim(exprDataBortezomib)
exprDataBortezomib[1:4,1:4]
boxplot(exprDataBortezomib[1:4,])
# 264個(gè)樣品的表達(dá)量矩陣
table(studyResponse)
# 264個(gè)病人的結(jié)局 （藥物處理后）
PGx_Responder = IE PGx_Responder = NR PGx_Responder = R
25 126 113
具體是什么癌癥什么病人就需要看該文章以及該數(shù)據(jù)集的來(lái)源文獻(xiàn)啦，但是藥物處理結(jié)局事件是很明顯的，主要是區(qū)分R和NR，應(yīng)該是有無(wú)響應(yīng)的簡(jiǎn)單情況。
library(pRRophetic)
# exprMatCcle (ccleData)
data(ccleData)
dim(exprMatCcle)
exprMatCcle[1:4,1:4]
boxplot(exprMatCcle[,1:4])
# 1037個(gè)細(xì)胞系的18988個(gè)基因組成的 表達(dá)量矩陣
head(sensDataCcle)
tmp=as.data.frame(table(sensDataCcle[,2:3]) )
length(unique(sensDataCcle$Compound))
length(unique(sensDataCcle$Primary.Cell.Line.Name))
可以看到ccle數(shù)據(jù)庫(kù)里面的表達(dá)量矩陣信息很豐富，有1037個(gè)細(xì)胞系的18988個(gè)基因表達(dá)。
但是藥物處理信息只有24個(gè)藥物，五百多個(gè)細(xì)胞系，如下所示：
> 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)里面的也是一千多個(gè)細(xì)胞系的2萬(wàn)個(gè)基因的表達(dá)量矩陣，關(guān)鍵是它藥物超過(guò)了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)
這個(gè)核心函數(shù)pRRopheticPredict，有非常多的參數(shù)，不過(guò)絕大部分都是有默認(rèn)值，如下所示：
pRRopheticPredict(testMatrix, drug, tissueType = "all", batchCorrect = "eb",
powerTransformPhenotype = TRUE, removeLowVaryingGenes = 0.2,
minNumSamples = 10, selection = -1, printOutput = TRUE,
removeLowVaringGenesFrom = "homogenizeData", dataset = "cgp2014")
我們首先使用 默認(rèn)的 cgp2014數(shù)據(jù)集，來(lái)預(yù)測(cè) 前面的 硼替佐米（Bortezomib）藥物信息，看看是否一致。
library(pRRophetic)
library(ggplot2)
data(bortezomibData)
dim(exprDataBortezomib)
exprDataBortezomib[1:4,1:4]
boxplot(exprDataBortezomib[,1:4])
# 264個(gè)樣品的表達(dá)量矩陣
table(studyResponse)
predictedPtype <- pRRopheticPredict(testMatrix=exprDataBortezomib,
drug="Bortezomib",
tissueType = "all",
batchCorrect = "eb",
selection=1,
dataset = "cgp2014")
如果是實(shí)際使用，我們應(yīng)該是讀入自己的表達(dá)量哦，不過(guò)我們這里是演示如何使用它，所以直接用exprDataBortezomib矩陣來(lái)測(cè)試pRRophetic包的核心函數(shù)pRRopheticPredict即可。
因?yàn)榫皖A(yù)測(cè)一個(gè)藥物，所以速度很快：
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
讓我們看看，預(yù)測(cè)情況和時(shí)間情況的差異：
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')
如下所示：
總體來(lái)說(shuō)，確實(shí)是R的響應(yīng)組，里面的病人預(yù)測(cè)得到的藥物越敏感！
這個(gè)時(shí)候因?yàn)槭敲總€(gè)藥物都需要走前面的 pRRophetic包的核心函數(shù)pRRopheticPredict，所以可以寫(xiě)循環(huán)啦，而且可以加入并行機(jī)制。
這個(gè)時(shí)候我們選擇 cgp2016數(shù)據(jù)集，全部的代碼如下所示：
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來(lái)返回計(jì)算所需時(shí)間
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) # 整合結(jié)果
stopCluster(cl) # 關(guān)閉集群
})
這個(gè)函數(shù)運(yùn)行取決于你的計(jì)算資源，需要半個(gè)小時(shí)左右。慢慢等吧，喝一杯咖啡吧，如果可以的話希望你在咱們《生信技能樹(shù)》公眾號(hào)任意教程末尾打賞一杯咖啡也行，我們一起慢慢喝，慢慢等！
如果是自己的真實(shí)表達(dá)量矩陣，需要走 pRRophetic包的核心函數(shù)pRRopheticPredict進(jìn)行預(yù)測(cè)，只需要讀入進(jìn)去，替換前面的 testMatrix=exprDataBortezomib 即可！
不同的數(shù)據(jù)庫(kù)資源作為函數(shù)的訓(xùn)練集，得到的結(jié)果必然是不一樣的哦！而且函數(shù)也可以調(diào)整很多參數(shù)。
這個(gè)時(shí)候其實(shí)可以看看 前面的2021年7月新鮮出爐的 藥物預(yù)測(cè)R包之oncoPredict 結(jié)果跟本次介紹的藥物預(yù)測(cè)R包之pRRophetic的一致性！
聯(lián)系客服
