oncoPredict包还有其他几个函数
💡专注R讲话在🩺生物医学中的使用
设为“星标”，精彩可以过
之前一经详备先容过pRRophetic包揣度药物明锐性了，可是这个包太陈腐了，我筹办好多东谈主会困在安设这一步，毕竟关于生手来说最难的便是R包安设了。
今天先容下oncoPredict，这个包是pRRophetic的升级版，使用门径和旨趣一模通常，仅仅换了以下进修数据汉典，也便是默许使用的数据库不通常了，其他王人是通常的。除此除外还增多了几个新的函数。
主邀功能是揣度药物反应和药物-基因关联，github的描写：An R package for drug response prediction and drug-gene association prediction. The prepared GDSC and CTRP matrices for the calcPhenotype() are located in the oncoPredict OSF.
oncoPredict包的作家说他们会握续更新这个包，你们信吗？
图片
CRAN和github王人袒露前次更新是2年前了~
安设
calcPhenotype
IDWAS
CNV
snv
GLDS
安设可以平直从cran安设了：
install.packages("oncoPredict")可是安设后不成平直使用，需要下载这个包配套的进修数据，不外表面上你我方准备进修数据亦然可以的~等我有技艺尝试一下。
配套进修数据下载地址：https://osf.io/c6tfx/，一共是600多M大小。
library(oncoPredict)## ## ## Attaching package: 'oncoPredict'## The following objects are masked from 'package:pRRophetic':## ## calcPhenotype， homogenizeData， summarizeGenesByMean
下载下来的进修数据便是以下几个，主淌若GDSC1和GDSC2，以及CTRP的数据：
图片
CTRP提供了RPKM和TPM两种方法的抒发矩阵，GDSC则是芯片方法的。
简略看下数据是什么款式的：
GDSC2_Expr <- readRDS(file="../000files/DataFiles/Training Data/GDSC2_Expr (RMA Normalized and Log Transformed).rds")GDSC2_Res <- readRDS(file = "../000files/DataFiles/Training Data/GDSC2_Res.rds")dim(GDSC2_Expr) #17419，805## [1] 17419 805dim(GDSC2_Res) #805，198## [1] 805 198GDSC2_Expr[1:4，1:4]## COSMIC_906826 COSMIC_687983 COSMIC_910927 COSMIC_1240138## 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.063701GDSC2_Res[1:4，1:4]## Camptothecin_1003 Vinblastine_1004 Cisplatin_1005## COSMIC_906826 -1.152528 -1.566172 7.017559## COSMIC_687983 -1.263248 -4.292974 3.286848## COSMIC_910927 -3.521093 -5.008028 2.492692## COSMIC_1240138 1.976381 NA NA## Cytarabine_1006## COSMIC_906826 2.917980## COSMIC_687983 2.790819## COSMIC_910927 -1.082517## COSMIC_1240138 NA
GDSC2_Expr是活动的抒发矩阵方法，行是基因，列是细胞系，一共17419个基因，805个细胞系；
GDSC2_Res是每个细胞系对每个药物的IC50值，行是细胞系，列是药物，一共805种细胞系，198个药物。
再看下CTRP的数据，以TPM为例：
CTRP2_Expr <- readRDS(file="../000files/DataFiles/Training Data/CTRP2_Expr (TPM， not log transformed).rds")CTRP2_Res <- readRDS(file = "../000files/DataFiles/Training Data/CTRP2_Res.rds")dim(CTRP2_Expr)## [1] 51847 829dim(CTRP2_Res)## [1] 829 545CTRP2_Expr[1:4，1:4]## CVCL_1045 CVCL_1046 CVCL_7937 CVCL_7935## DDX11L1 0.0000000 0.10511906 0.0000000 0.1444173## WASH7P 34.3887920 28.39068559 14.8423404 14.5556952## MIR1302-11 0.1167793 0.02903022 0.4054605 0.5185439## FAM138A 0.0000000 0.02432715 0.5362218 0.3674906CTRP2_Res[1:4，1:4]## CIL55 BRD4132 BRD6340 BRD9876## CVCL_1045 14.504 14.819 14.101 14.657## CVCL_1046 14.982 12.110 14.529 14.702## CVCL_7937 14.388 12.091 14.742 13.230## CVCL_7935 14.557 NA 14.546 13.258
和GDSC的数据方法一模通常的，亦然抒发矩阵和药敏信息，一共51847个基因，应该是包括mRNA和非编码RNA的，829个细胞系，545种药物。
calcPhenotype这个包主要便是3个函数，最蹙迫的一个便是calcPhenotype了。pRRophetic也有calcPhenotype函数，其实王人是通常的用法，各个参数咱们在上一篇中也先容过了，这里就未几说了。
咱们使用TCGA-BLCA的lncRNA的抒发矩阵进行演示，因为CTPR的进修数据是转录组，而且包含了非编码RNA，是以是可以用来揣度lncRNA的药物明锐性的。
底下这个示例数据我放在了粉丝QQ群，需要的加群下载即可。
最初加载lncRNA的抒发矩阵：
load(file = "../000files/testExpr_BLCA.rdata")
这个示例抒发矩阵一共12162个lncRNA，414个样本，何况被分为高风险组和低风险组，抒发矩阵已历程了log2颐养：
table(sample_group)## sample_group## high low ## 207 207dim(testExpr)## [1] 12162 414testExpr[1:4，1:4]## TCGA-ZF-A9R5-01A-12R-A42T-07 TCGA-E7-A7PW-01A-11R-A352-07## AL627309.1 0.01705863 0.00000000## AL627309.2 0.09572744 0.00000000## AL627309.5 0.08029061 0.01252150## AC114498.1 0.00000000 0.06007669## TCGA-4Z-AA7N-01A-11R-A39I-07 TCGA-DK-A2I1-01A-11R-A180-07## AL627309.1 0.00000000 0.000000000## AL627309.2 0.00000000 0.000000000## AL627309.5 0.04722951 0.004103184## AC114498.1 0.00000000 0.000000000
底下就可以一次策动整个545种药物的明锐性了，速率很慢，而且成果只可保存到现时责任目次下的calcPhenotype_Output文献夹中。
calcPhenotype(trainingExprData = CTRP2_Expr， trainingPtype = CTRP2_Res， testExprData = as.matrix(testExpr)，#需要matrix batchCorrect = 'eb'， #IC50是对数颐养的，是以抒发矩阵也用对数颐养过的 powerTransformPhenotype = F， minNumSamples = 20， printOutput = T， removeLowVaryingGenes = 0.2， removeLowVaringGenesFrom = "homogenizeData" )
成果就一个文献，便是每个样本对每一个药物的IC50值，读取进来望望：
res <- read.csv("./calcPhenotype_Output/DrugPredictions.csv")dim(res)## [1] 414 546res[1:4，1:4]## X CIL55 BRD4132 BRD6340## 1 TCGA-ZF-A9R5-01A-12R-A42T-07 14.10501 13.45597 14.30078## 2 TCGA-E7-A7PW-01A-11R-A352-07 13.81229 13.19243 14.10742## 3 TCGA-4Z-AA7N-01A-11R-A39I-07 14.02137 13.02397 14.27905## 4 TCGA-DK-A2I1-01A-11R-A180-07 14.48769 13.54434 14.41776这么就得回了414个样本对每种药物的IC50值。
有了这个成果，咱们就可以取出感好奇的药物，可视化IC50值在不同组间的各别，和之前先容的一模通常，咱们这里松弛取前10个药物：
伊吾县科艾净水器有限公司 0.55) 0px 2px 10px;line-height: 1.6 !important;">library(tidyr)library(dplyr)library(ggplot2)library(ggpubr)library(ggsci)res %>% select(1:11) %>% bind_cols(sample_group = sample_group) %>% pivot_longer(2:11， 常州依梯凯动力机械有限公司names_to = "drugs"， 福建吉美电子有限公司values_to = "ic50") %>% ggplot(.， aes(sample_group，ic50))+ geom_boxplot(aes(fill=sample_group))+ 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包还有其他几个函数，咱们简略先容一下。
IDWASIDWAS 门径是一种推广的药物反馈值填补门径，可以绵薄地在东谈主群中进行生物记号物的卤莽。IDWAS 在观念和推论上与GWAS相似，即通过使用R中的线性模子慑服填补的药物反馈值与体细胞突变或拷贝数变异之间的关联，以筹办药物-基因互相作用并识别药物反馈的生物记号物。垄断临床数据集的大样本量，这种门径可以发现细胞所有这个词据齐集莫得的新关系。此外，由于使用了患者的基因组数据，所找到的任何生物记号物王人与该患者群体谋划。
简略来说便是把柄药敏信息和表型信息识别药物的靶点。表型信息可以是突变数据或者拷贝数变异数据。
CNV咱们这里先试试拷贝数变异，平直使用TCGA-BLCA的拷贝数变异数据，平直使用easyTCGA，1行代码处分拷贝数变异数据的下载：
rm(list = ls())library(oncoPredict)library(easyTCGA)getcnv("TCGA-BLCA")下载完成后平直加载数据即可：
load(file = "G:/easyTCGA_test/output_cnv/TCGA-BLCA_CNV.rdata")blca_cnv <- datahead(blca_cnv)## # A tibble: 6 × 7## GDC_Aliquot Chromosome Start End Num_Probes Segment_Mean Sample## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> ## 1 64caa1a2-b01f-404c-ba… 1 3.30e6 1.57e8 71640 0.0136 TCGA-…## 2 64caa1a2-b01f-404c-ba… 1 1.57e8 1.57e8 2 -1.66 TCGA-…## 3 64caa1a2-b01f-404c-ba… 1 1.57e8 1.86e8 18893 0.0124 TCGA-…## 4 64caa1a2-b01f-404c-ba… 1 1.86e8 1.86e8 2 -1.79 TCGA-…## 5 64caa1a2-b01f-404c-ba… 1 1.86e8 2.48e8 39130 0.0141 TCGA-…## 6 64caa1a2-b01f-404c-ba… 2 4.81e5 2.42e8 132068 0.0091 TCGA-…
然后使用map_cnv颐养一下方法。
The mapping is accomplished by intersecting the gene with the overlapping CNV level. If the gene isn't fully #captured by the CNV， an NA will be assigned.
map_cnv(blca_cnv) 403 genes were dropped because they have exons located on both strands of the same reference sequence or on more than one reference sequence， so cannot be represented by a single genomic range. Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList object， or use suppressMessages() to suppress this message.
生成的map.RData文献会自动保存在现时责任目次下，加载进来：
load(file = "./map.RData")data<-as.data.frame(t(theCnvQuantVecList_mat))
接下来是准备药敏数据，这个数据咱们可以平直使用上头的calcPhenotype得回的成果：
drug_prediction<-as.data.frame(read.csv('./calcPhenotype_Output/DrugPredictions.csv'， header=TRUE， row.names=1))dim(drug_prediction) ## [1] 414 545drug_prediction[1:4，1:4]## CIL55 BRD4132 BRD6340 BRD9876## TCGA-ZF-A9R5-01A-12R-A42T-07 14.10501 13.45597 14.30078 14.05845## TCGA-E7-A7PW-01A-11R-A352-07 13.81229 13.19243 14.10742 13.79403## TCGA-4Z-AA7N-01A-11R-A39I-07 14.02137 13.02397 14.27905 12.32300## TCGA-DK-A2I1-01A-11R-A180-07 14.48769 13.54434 14.41776 12.53436然后就可以使用idwas策动了：
idwas(drug_prediction=drug_prediction，罐头食品 data=data， n=10， cnv=T)
成果文献会自动保存在现时责任目次中。
cnvPvalue <- read.csv("CnvTestOutput_pVals.csv"，row.names = 1)dim(cnvPvalue)## [1] 545 402cnvPvalue[15:20，1:4]## SNORD123 MIR875 MKRN2OS LOC100130075## teniposide 0.157018358 0.8511244 0.08231962 0.05243646## sildenafil 0.057896196 0.1650301 0.11829531 0.16373477## simvastatin 0.003648419 0.1894510 0.16468764 0.24833938## parbendazole 0.309342093 0.6758557 0.09965932 0.74859683## procarbazine 0.784280696 0.3716664 0.11423190 0.42647709## curcumin 0.297120449 0.6824545 0.73916619 0.23619591行是药物，列是找出来的靶点，中间是P值，应该是P值小于0.05诠释故好奇吧，详备情况群众我方阅读文献了解。
还有一个beta-value矩阵，通常的方法，就不展示了。
snv如果你使用的表型数据是突变数据，那就更简略了，这里如故以TCGA-BLCA的突变数据为例。
平直使用easyTCGA1行代码下载突变数据：
library(easyTCGA)getsnvmaf("TCGA-BLCA")药敏数据如故和上头通常的，把突变数据平直加载进来就可以进行策动了。
因为这里使用的是突变数据不是拷贝数变异，是以需要革新参数cnv=F，其他竣工通常：
load(file = "G:/easyTCGA_test/output_snv/TCGA-BLCA_maf.rdata")idwas(drug_prediction=drug_prediction， data=data， n=10， cnv=F)
成果文献亦然自动保存在现时责任目次中。
GLDS作家在原文中对GLDS的解说:
咱们之前就“GLDS 景观”过火对生物记号物卤莽的影响进行了报谈。GLDS是指在一个东谈主群（细胞系或患者）中，岂论采取何种颐养，个体频繁可以进展出更明锐或更耐药的趋势。正如原始论文所示，这一景观与多重耐药（MDR）谋划，但并不竣工斟酌。对这一变量进行修订可以更具体地慑服与特定药物谋划的生物记号物。
也便是说这个函数可以把柄你提供的表型数据和药敏数据，推测出与药物最谋划的靶点，嗅觉和上头的函数作用差未几。
表型数据也可以是拷贝数变异或者突变数据。
底下是一个示例。
最初加载药敏数据，我这里用的是示例数据和示例代码。
rm(list = ls())trainingPtype = readRDS(file = "../000files/DataFiles/Training Data/GDSC2_Res.rds")class(trainingPtype)## [1] "matrix" "array"dim(trainingPtype)## [1] 805 198trainingPtype[1:4，1:4]## Camptothecin_1003 Vinblastine_1004 Cisplatin_1005## COSMIC_906826 -1.152528 -1.566172 7.017559## COSMIC_687983 -1.263248 -4.292974 3.286848## COSMIC_910927 -3.521093 -5.008028 2.492692## COSMIC_1240138 1.976381 NA NA## Cytarabine_1006## COSMIC_906826 2.917980## COSMIC_687983 2.790819## COSMIC_910927 -1.082517## COSMIC_1240138 NA
行是细胞系，列是药物，亦然一个抒发矩阵的方法。
GLDS不成有缺失值，是以先进行缺失值插补，使用自带的completeMatrix函数，默许进行50次迭代，巨慢！！是以我建造成3了。
# 缺失值插补completeMatrix(trainingPtype， nPerms = 3)
成果文献complete_matrix_output_GDSCv2.txt会自动保存在现时责任中.
# 读取插补后的药敏数据cm<-read.table('complete_matrix_output.txt'， header=TRUE， row.names=1) dim(cm)## [1] 805 198cm[1:4，1:4]## Camptothecin_1003 Vinblastine_1004 Cisplatin_1005## COSMIC_906826 -1.152528 -1.566172 7.017559## COSMIC_687983 -1.263248 -4.292974 3.286848## COSMIC_910927 -3.521093 -5.008028 2.492692## COSMIC_1240138 1.976381 -3.957938 3.354192## Cytarabine_1006## COSMIC_906826 2.917980## COSMIC_687983 2.790819## COSMIC_910927 -1.082517## COSMIC_1240138 1.697800为了进行策动，咱们把细胞系和药物的名字王人改一下，主淌若去掉细胞系的前缀和药物的后缀。
读取细胞系的信息，内部有细胞系名字的详备信息：
cellLineDetails<-readxl::read_xlsx('../000files/DataFiles/GLDS/GDSCv2/Cell_Lines_Details.xlsx')dim(cellLineDetails)## [1] 1002 13cellLineDetails[1:4，1:4]## # A tibble: 4 × 4## `Sample Name` `COSMIC identifier` `Whole Exome Sequencing (WES)`## <chr> <dbl> <chr> ## 1 A253 906794 Y ## 2 BB30-HNC 753531 Y ## 3 BB49-HNC 753532 Y ## 4 BHY 753535 Y ## # ℹ 1 more variable: `Copy Number Alterations (CNA)` <chr>替换药敏数据的行名，也便是细胞系名字：
newRows <- substring(rownames(cm)，8) #Remove 'COSMIC'...keep the numbers after COSMIC.indices<-match(as.numeric(newRows)， as.vector(unlist(cellLineDetails[，2]))) #Refer to the cell line details file to make this replacement.newNames<-as.vector(unlist(cellLineDetails[，1]))[indices] #Reports the corresponding cell line namesrownames(cm)<-newNamesdim(cm)## [1] 805 198cm[1:4，1:4]## Camptothecin_1003 Vinblastine_1004 Cisplatin_1005 Cytarabine_1006## CAL-120 -1.152528 -1.566172 7.017559 2.917980## DMS-114 -1.263248 -4.292974 3.286848 2.790819## CAL-51 -3.521093 -5.008028 2.492692 -1.082517## H2869 1.976381 -3.957938 3.354192 1.697800
把柄作家提供的信息，把药物名字也改一下：
fix <- readxl::read_xlsx('../000files/DataFiles/GLDS/GDSCv2/gdscv2_drugs.xlsx')## New names:## · `` -> `...2`fix<-as.vector(unlist(fix[，1]))head(fix)## [1] "Camptothecin" "Vinblastine" "Cisplatin" "Cytarabine" "Docetaxel" ## [6] "Gefitinib"colnames(cm)<-as.vector(fix)drugMat<-as.matrix(cm) #Finally， set this object as the drugMat parameter. dim(drugMat) #805 samples vs. 198 drugs## [1] 805 198drugMat[1:4，1:4]## Camptothecin Vinblastine Cisplatin Cytarabine## CAL-120 -1.152528 -1.566172 7.017559 2.917980## DMS-114 -1.263248 -4.292974 3.286848 2.790819## CAL-51 -3.521093 -5.008028 2.492692 -1.082517## H2869 1.976381 -3.957938 3.354192 1.697800到这里这个药敏数据终于准备好了！
底下要准备marker矩阵，也便是突变或者拷贝数变异数据。这里亦然用的示例数据，是一个泛癌的数据，同期包含CNV和突变信息：
mutationMat<-read.csv('../000files/DataFiles/GLDS/GDSCv2/GDSC2_Pan_Both.csv')mutationMat<-mutationMat[，c(1，6，7)] #Index to these 3 columns of interest.colnames(mutationMat) #"cell_line_name" "genetic_feature" "is_mutated" ## [1] "cell_line_name" "genetic_feature" "is_mutated"head(mutationMat)## cell_line_name genetic_feature is_mutated## 1 CAL-29 CDC27_mut 0## 2 CAL-29 CDC73_mut 0## 3 CAL-29 CDH1_mut 0## 4 CAL-29 CDK12_mut 0## 5 CAL-29 CDKN1A_mut 0## 6 CAL-29 CDKN1B_mut 0上头这个文献中有一些重迭的cell_line_name和genetic_feature对，先去掉重迭的：
vec<-c()for (i in 1:nrow(mutationMat)){ vec[i]<-paste(mutationMat[i，1]，mutationMat[i，2]， sep=' ')}nonDupIndices<-match(unique(vec)， vec)mutationMat2<-mutationMat[nonDupIndices，]dim(mutationMat2)## [1] 584051 3head(mutationMat2)## cell_line_name genetic_feature is_mutated## 1 CAL-29 CDC27_mut 0## 2 CAL-29 CDC73_mut 0## 3 CAL-29 CDH1_mut 0## 4 CAL-29 CDK12_mut 0## 5 CAL-29 CDKN1A_mut 0## 6 CAL-29 CDKN1B_mut 0把空值去掉，再变为宽数据，使得行名是细胞系名字：
library(tidyverse)good<-(mutationMat2[，2]) != ""mutationMat3<-mutationMat2[good，]mutationMat4<-mutationMat3 %>% pivot_wider(names_from=genetic_feature， values_from=is_mutated)rownames(mutationMat4)<-as.vector(unlist(mutationMat4[，1])) #Make cell lines the rownames...right now they are column 1.cols<-rownames(mutationMat4)mutationMat4<-as.matrix(t(mutationMat4[，-1]))dim(mutationMat4)## [1] 1315 1389mutationMat4[1:4，1:4]## [，1] [，2] [，3] [，4]## CDC27_mut "0" "0" "0" NA ## CDC73_mut "0" "0" "0" NA ## CDH1_mut "0" "0" "0" NA ## CDK12_mut "0" "0" "0" NA
把字符型变为数值型，把NA也酿成0：
#Make sure the matrix is numeric.mutationMat<-mutationMat4mutationMat4<-apply(mutationMat4， 2， as.numeric)rownames(mutationMat4)<-rownames(mutationMat)markerMat<-mutationMat4markerMat[1:4，1:4]## [，1] [，2] [，3] [，4]## CDC27_mut 0 0 0 NA## CDC73_mut 0 0 0 NA## CDH1_mut 0 0 0 NA## CDK12_mut 0 0 0 NA# replace all non-finite values with 0markerMat[!is.finite(markerMat)] <- 0colnames(markerMat)<-colsdim(markerMat) #1315 1389## [1] 1315 1389markerMat[1:4，1:4]## CAL-29 CAL-33 697 CCNE1## CDC27_mut 0 0 0 0## CDC73_mut 0 0 0 0## CDH1_mut 0 0 0 0## CDK12_mut 0 0 0 0#保存一下#write.table(markerMat， file='markerMat.txt')markerMat<-as.matrix(read.table('markerMat.txt'， header=TRUE， row.names=1))到这里marker矩阵也准备好了。
临了还需要一个drug relatedness file文献，前边更名字便是为了和这里的名字保握一致。
drugRelatedness: This file is GDSC's updated drug relatedness file (obtained from bulk data download/all compounds screened/compounds-annotation).
drugRelatedness <- read.csv("../000files/DataFiles/GLDS/GDSCv2/screened_compunds_rel_8.2.csv")drugRelatedness<-drugRelatedness[，c(3，6)]colnames(drugRelatedness) #"DRUG_NAME" "TARGET_PATHWAY"## [1] "DRUG_NAME" "TARGET_PATHWAY"head(drugRelatedness)## DRUG_NAME TARGET_PATHWAY## 1 Erlotinib EGFR ## 2 Rapamycin PI3K/MTOR ## 3 Sunitinib RTK ## 4 PHA.665752 RTK ## 5 MG-132 Protein stability and degradation## 6 Paclitaxel Mitosis然后就可以驱动glds了：
glds(drugMat， drugRelatedness， markerMat， minMuts=5， additionalCovariateMatrix=NULL， threshold=0.7)
成果亦然自动保存在现时责任目次中。
文献中的这张图应该便是把柄上头的成果画出来的，群众我方议论下喽~
图片
适度。
本文对后两个函数仅仅按照匡助文档驱动了一遍罐头食品，如果要详备了解各式细节，提出群众去阅读文献哈~
本站仅提供存储事业，整个施行均由用户发布，如发现存害或侵权施行，请点击举报。