成人午夜直播,天天弄天天操

有了前面的教程：藥物預(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ù)庫資源。也介紹了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》

但是它一直沒有把包提交到cran或者bioconductor，而是自己維護(hù)：

website (http://genemed.uchicago.edu/~pgeeleher/pRRophetic)
GitHub (https://github.com/paulgeeleher/pRRophetic).

文章短小精悍，就5個(gè)參考文獻(xiàn)，而且它并不使用我們介紹的Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 兩個(gè)數(shù)據(jù)庫資源，而是使用 Cancer Genome Project (CGP)計(jì)劃里面的表達(dá)量矩陣和藥物處理信息，該數(shù)據(jù)庫有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èn)識(shí)R包pRRophetic的數(shù)據(jù)

安裝成功之后，在你的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ù)。不過它畢竟是沒有提交到cran或者bioconductor，還是有點(diǎn)擔(dān)心?。?/p>

首先查看硼替佐米（Bortezomib）

Bortezomib (PS-341) 是一種可逆性和選擇性的蛋白酶體 (proteasome) 抑制劑，通過靶向蘇氨酸殘基有效抑制 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ù)集的來源文獻(xiàn)啦，但是藥物處理結(jié)局事件是很明顯的，主要是區(qū)分R和NR，應(yīng)該是有無響應(yīng)的簡(jiǎn)單情況。

然后查看ccle數(shù)據(jù)資源

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ù)庫里面的表達(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)

可以看到 Cancer Genome Project (CGP)里面的也是一千多個(gè)細(xì)胞系的2萬個(gè)基因的表達(dá)量矩陣，關(guān)鍵是它藥物超過了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)

pRRophetic包的核心函數(shù)pRRopheticPredict的使用

這個(gè)核心函數(shù)pRRopheticPredict，有非常多的參數(shù)，不過絕大部分都是有默認(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ù)集，來預(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á)量哦，不過我們這里是演示如何使用它，所以直接用exprDataBortezomib矩陣來測(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')

如下所示：

總體來說，確實(shí)是R的響應(yīng)組，里面的病人預(yù)測(cè)得到的藥物越敏感！

預(yù)測(cè)全部的藥物

這個(gè)時(shí)候因?yàn)槭敲總€(gè)藥物都需要走前面的 pRRophetic包的核心函數(shù)pRRopheticPredict，所以可以寫循環(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來返回計(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í)左右。慢慢等吧，喝一杯咖啡吧，如果可以的話希望你在咱們《生信技能樹》公眾號(hào)任意教程末尾打賞一杯咖啡也行，我們一起慢慢喝，慢慢等！

如果是自己的真實(shí)表達(dá)量矩陣，需要走 pRRophetic包的核心函數(shù)pRRopheticPredict進(jìn)行預(yù)測(cè)，只需要讀入進(jìn)去，替換前面的 testMatrix=exprDataBortezomib 即可！

不同的數(shù)據(jù)庫資源作為函數(shù)的訓(xùn)練集，得到的結(jié)果必然是不一樣的哦！而且函數(shù)也可以調(diào)整很多參數(shù)。

這個(gè)時(shí)候其實(shí)可以看看前面的2021年7月新鮮出爐的藥物預(yù)測(cè)R包之oncoPredict 結(jié)果跟本次介紹的藥物預(yù)測(cè)R包之pRRophetic的一致性！

本站僅提供存儲(chǔ)服務(wù)，所有內(nèi)容均由用戶發(fā)布，如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊舉報(bào)。