液体冷却
💡专注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
本站仅提供存储管事,所有这个词实践均由用户发布,如发现存害或侵权实践,请点击举报。