| 发布日期：2024-08-09 09:08 点击次数：174 |
💡专注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的数据：
图片
appCTRP提供了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个药物：
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包还有其他几个函数，咱们概况先容一下。
福彩快乐8第2024175期(上周三)开奖回顾：07 09 12 15 17 19 32 33 40 47 48 49 55 58 62 65 66 69 70 73，其中奖号五区比为4：3：4：4：5。
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)
成果亦然自动保存在刻下责任目次中。
文献中的这张图应该便是凭证上头的成果画出来的，环球我方谈论下喽~
图片
放弃。
本文对后两个函数仅仅按照匡助文档启动了一遍软件开发价格，如果要详备了解各式细节，暴戾环球去阅读文献哈~
本站仅提供存储就业，悉数实质均由用户发布，如发现存害或侵权实质，请点击举报。