| 发布日期：2024-08-09 07:36 点击次数：118 |
💡专注R言语在🩺生物医学中的使用
设为“星标”，精彩可以过
[扫码下载app，中过数字彩1千万以上的专家都在这儿！]
药物明锐性分析是生信数据挖掘常用的手段之一，当今作念药敏分析最常见的即是两个R包：pRRophetic和oncoPredict。
这两个包的作家皆是吞并个东谈主，oncoPredict可以看作念是pRRophetic的升级版。两个R包的使用基本上是相通的念念路，只不外使用的西宾数据集不同良友。
在先容R包的使用之前，需要全球先了解一下常用的药物明锐性数据库，最佳是去到这些数据库的主页望望或者读一读关系的文献，对这些数据库有一个简短的了解。
常用药敏数据库
pRRophetic神志学先容
装置
瞻望不同组别患者对化疗药物的明锐性
其他示例
pRRopheticCV
使用CCLE示例数据
自界说西宾集
自带数据探索
瞻望全的药物的明锐性
参考
常用药敏数据库药敏数据库尽头多，但最常用的无非即是GDSC/CTRP/CCLE等，在珠江肿瘤公众号中早就先容过这些数据库了，是以我就不访佛了，全球可以参考以下纠合。
以下纠合先容了GDSC、CTRP、CCLE、NCI-60、DepMap、Pharmacodb等数据库，长短常棒的参考汉典：
肿瘤药敏多组学数据库(GDSC)概览肿瘤药敏多组学数据库(GDSC)的数据先容和获得GDSC与其他药敏多组学数据库GDSC与CELL数据库的药物基因组学一致性靶点抒发水平可算作靶向药物明锐性的谋划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 264exprDataBortezomib[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] 141head(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 1018cgp2016ExprRma[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] 251head(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) # drugData2016data(cgp2016ExprRma) # cgp2016ExprRmadata(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
本站仅提供存储做事，悉数本色均由用户发布，如发现存害或侵权本色，请点击举报。