新闻资讯
频繁你的抒发矩阵也会分组的
💡专注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
本站仅提供存储奇迹，通盘履行均由用户发布，如发现存害或侵权履行，请点击举报。下一篇：没有了