液体冷却
💡专注R话语在🩺生物医学中的使用
设为“星标”,精彩可以过
药物敏锐性分析是生信数据挖掘常用的时间之一,当今作念药敏分析最常见的便是两个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
本站仅提供存储奇迹,通盘本色均由用户发布,如发现存害或侵权本色,请点击举报。