oncoPredict 是一款用来预测药物反应的R包，背景知识有很多介绍的了，这里介绍下真实的使用场景 以及 后续联合基因表达，分子分型 或者 预后模型等的联合。
一 载入数据，R包
1，安装R包
oncoPredict是CRAN中的包，直接install.packages安装，但是大概率会遇到缺少数据库相关R包的情况，根据提示安装即可。
install.packages("oncoPredict")#BiocManager::install("GenomicFeatures")#BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")library(oncoPredict)
2，输入数据
oncoPredict需要3个输入数据：（1）训练集表达数据-GDSC2 expression data （2）训练集药物反应数据-GDSC2 pType data （3）待预测的表达数据-our expression data 。
oncoPredict包提供了两种当前经典的药物预测训练集-Cancer Therapeutics Response Portal (CTRP-V2) 和 Genomics of Drug Sensitivity in Cancer (GDSC-V2) 。可以直接使用R包oncoPredict整理好的这两个数据库的rdata文件，下载链接https://osf.io/c6tfx/files/osfstorage。
CTRP可能是测序方式，提供了 FPKM和TPM数据。GDSC应该是芯片，提供了 RMA Normalized and Log Transformed。
3，读取三个数据集
训练集使用GDSC-V2 ，预测集使用之前使用的SKCM的表达矩阵
## 药物训练集trainingExprData=readRDS(file='./DataFiles/Training Data/GDSC2_Expr (RMA Normalized and Log Transformed).rds')dim(trainingExprData) #trainingExprData[1:4,1:4]#IMPORTANT note: here I do e^IC50 since the IC50s are actual ln values/log transformed already, and the calcPhenotype function Paul #has will do a power transformation (I assumed it would be better to not have both transformations)trainingPtype<-exp(trainingPtype)trainingPtype[1:4,1:4]## 转录组数据load("RNAseq.SKCM.RData")fpkm[1:4,1:4]fpkm_in <- fpkm %>%column_to_rownames("gene") %>%as.matrix()testExprData = fpkm_intestExprData[1:4,1:4]
trainingPtype是所有细胞系的不同药物的反应矩阵，trainingExprData是对应细胞系基因的表达矩阵，二者的样本名称是对应的。
testExprData是待预测的自己的基因表达矩阵数据，行名为基因，列为样本。
二 预测药物反应
使用calcPhenoty函数进行预测。
1，配置参数
calcPhenotype的参数较多，没啥特殊需求建议使用默认的 ，参数解释相见官方脚步的参数说明：https://github.com/cran/oncoPredict/blob/master/R/CALCPHENOTYPE.R
batchCorrect<-"eb"powerTransformPhenotype<-TRUEremoveLowVaryingGenes<-0.2removeLowVaringGenesFrom<-"homogenizeData"minNumSamples=10selection<- 1printOutput=TRUEpcr=FALSEreport_pc=FALSEcc=FALSErsq=FALSEpercent=80
2，calcPhenotype计算
calcPhenotype的结果会默认输出在当前目录下的calcPhenotype_Output文件夹，结果输出为DrugPredictions.csv文件。
多次分析的话 注意覆盖原结果的问题。
calcPhenotype(trainingExprData=trainingExprData,trainingPtype=trainingPtype,testExprData=testExprData,batchCorrect=batchCorrect,powerTransformPhenotype=powerTransformPhenotype,removeLowVaryingGenes=removeLowVaryingGenes,minNumSamples=minNumSamples,selection=selection,printOutput=printOutput,pcr=pcr,removeLowVaringGenesFrom=removeLowVaringGenesFrom,report_pc=report_pc,cc=cc,percent=percent,rsq=rsq)#读取结果resultPtype <- read.csv('./calcPhenotype_Output/DrugPredictions.csv', header = T ,stringsAsFactors = F ,check.names = F)resultPtype[1:4, 1:4]# Camptothecin_1003 Vinblastine_1004 Cisplatin_1005#1 TCGA-EE-A2GJ-06A 0.355799775 0.03555260 51.93174#2 TCGA-EE-A2GI-06A 0.009369448 0.02313505 11.29654#3 TCGA-WE-A8ZM-06A 0.083040399 0.04688820 76.53759#4 TCGA-DA-A1IA-06A 0.152475833 0.38765839 17.42450
如上就得到了每个样本的所有药物反应的结果。
三 基因+，分型+，预后+
既然得到的是每个样本的药物反应结果，当然就可以进行各种联合分析，这里分别以 目标基因 和 分子分型 进行示例。
1，重点基因表达量-相关性点图
重点基因可以来自于RNAseq|WGCNA-组学数据黏合剂，代码实战-一（尽）文（力）解决文献中常见的可视化图 找到的hub基因，RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线构建的预后模型的基因，可以是单细胞scRNA分析|Marker gene 可视化 以及 细胞亚群注释--你是如何人工注释的？找到的celltype 差异基因，可以是RNAseq|免疫浸润也杀疯了，cibersoert？xCELL？ESTIMATE？你常用哪一个免疫浸润程度等等。
names(resultPtype)[1] <- "sample"#以两个基因为示例fpkm_gene_drug <- fpkm_in %>%t() %>%as.data.frame() %>%dplyr::select(c("SPP1","TREM2")) %>%rownames_to_column("sample") %>%inner_join(resultPtype)fpkm_gene_drug[1:4,1:4]library(ggstatsplot)ggscatterstats(data = fpkm_gene_drug,x = SPP1,y = `BMS-754807_2171`,xlab = "SPP1",ylab = "BMS-754807_2171",marginal = TRUE,marginal.type = "densigram",margins = "both",xfill = "blue", # 分别设置颜色yfill = "#009E73",title = "Relationship between SPP1 and BMS-754807_2171",messages = FALSE) #BMS-754807_2171
更多参数详见R|散点图+边际图（柱形图，小提琴图），颜值UP
2，聚类/预后模型-分组箱线图
分组性状可以来自于RNAseq|组学分型-ConsensusClusterPlus（一致性聚类）， NMF（非负矩阵分解）的一致性聚类结果，可以是RNAseq|Lasso构建预后模型，绘制风险评分的KM 和 ROC曲线的风险分组，可以是关心的分类性状。
以无监督聚类结果为例
km <- read.csv("resultstrain.k=2.consensusClass.csv",header = T,stringsAsFactors = F ,check.names = F)names(km) <- c("sample","cluster")fpkm_km_drug <- km %>%inner_join(resultPtype)fpkm_km_drug[1:4,1:4]ggboxplot(fpkm_km_drug, x="cluster", y="Cisplatin_1005", width = 0.6,color = "black",#轮廓颜色fill="cluster",#填充palette = "npg",xlab = F, #不显示x轴的标签bxp.errorbar=T,#显示误差条bxp.errorbar.width=0.5, #误差条大小size=1, #箱型图边线的粗细legend = "right") #图例放右边 +
更多参数详见scRNA分析|自定义你的箱线图-统计检验，添加p值，分组比较p值
根据课题需要自行结合更多的分析方式。
◆ ◆ ◆ ◆ ◆
精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）
RNAseq纯生信挖掘思路分享？不，主要是送你代码！（建议收藏）