王炸！⽣信联合药物分析新套路，⽤这个数据库，⽴马⾼级！
cellMiner数据库的使⽤及药物敏感性预测
⼤家好，我是阿琛。在36策中，我们学习了“引药⽣变”的相关内容，其中2个问题4个类型，概括了基础科研中变量药物的全部研究策略。在内容，⾥⾯提到，⽆论是搞分⼦，还是研究药物，最终最好的结局也是做诊断和治疗的靶点。
当然，在⽣信研究中同样也不例外，⽆论是单基因研究，还是建模分析，整个⽂章的落脚点都可以往药物与治疗⽔平靠⼀靠，毫⽆疑问⽂章的内涵和吸引⼒都会提升⼀个档次。
CellMiner数据库，主要是通过国家癌症研究所癌症研究中⼼(NCI)所列出的60种癌细胞为基础⽽建⽴的。该数据库最初发表于2009年，后于2012年在Cancer Rearch杂志上进⾏了更新，题⽬为“CellMiner: a web-bad suite of genomic and pharmacologic tools to explore tran and drug patterns in the NCI-60 cell line t”。⼤家后期在使⽤该数据库记得应⽤相关⽂献。
NCI-60细胞系是⽬前使⽤最⼴泛的⽤于抗癌药物测试的癌细胞样本群。⼤家可以通过它查询到 NCI-60细胞系中已确认的22379个基因，以及20503个已分析的化合物的数据（包括多种已获美国⾷品和药物监督局批准的药物，以及临床试验中的药物分⼦）。
下⾯，我们来看下相关数据的下载和药物敏感性分析过程。
1.CellMiner数据库的使⽤
2. 点击Download Data Sets，进⼊数据下载界⾯；在下载界⾯中，可以看到两个不同的板块，分别是原始数据Raw Data Set和经过处理后的数据Procesd Data Set；
3. 在此，我们直接选择经过处理后的数据Procesd Data Set，勾选RNA表达数据（RNA: RNA-q）和药物数据（Compound activity: DTP NCI-60）；
4. 点击按钮Get Procesd Set，进⼊下载界⾯，点击即可保存；
5.下载完成后，将其放到⼯作⽬录下，解压，并分别提取其中的Excel⽂件放于当前⼯作⽬录下；魔法星球
接下来，我们需要对下载得到两个数据⽂件进⾏整理，以⽤于后续的药物敏感性评估。
2.药物数据的准备
2.1 读取药物相关数据
•自制发膜
•
library(readxl) rt1 < -read_excel( path= "DTP_NCI60_ZSCORE.xls", skip= 7)
⾸先，使⽤readxl包中的read_excel函数，读取药物相关的数据。由于前7⾏为注释信息，因此使⽤参数skip进⾏跳过前7⾏。
•
•
colnames( rt1) < -rt1[1,] rt1< -rt1[-1,-c(67,68)]
同时，将第⼀⾏作为列名，并去除末尾两列其余信息。
2.2 筛选药物标准
•
table(rt1$ `FDA status`)
使⽤table函数药物标准，结果显⽰，其中75种经过了临床试验，188种经过FDA批准。
•
•
•
rt1 <- rt1[rt1$` FDAstatus` % in% c( "FDA approved", "Clinical trial"),] rt1 <- rt1[,- c( 1, 3: 6)] write.table(rt1, file = "",p = "\t",row.names = F,quote = F)
为了保证分析结果的可靠性，选取经过临床试验（Clinical trial）和FDA批准（FDA approved）的药物结果。当然，你也可以选择将所有的药物进⾏保留。最终，⼀共得到263个药物结果，并将其保存为txt⽂件⽤于后续分析。
3. 基因表达数据的准备
•
•
•
保持身体健康英语
•
rt2<- read_excel(path = "RNA__RNA_q_composite_expression.xls", skip = 9) colnames(rt2) <- rt2[ 1,] rt2 <- rt2[- 1,-c( 2: 6)] write.table(rt2, file = "",p = "\t",row.names = F,quote = F)
同样的，读取表达数据，并对其进⾏整理和保存，⽤于后续的分析。
4. 药物敏感性分析
•
rm( list= ls)
4.1 引⽤包
•
•
library(impute) library(limma)
⾸先，加载分析使⽤的R包，包括impute包和limma包。
4.2 读取药物输⼊⽂件
•
•
•
•
•
•
•
王维送别诗rt <- read.table( "",p= "\t",header=T,check.names=F) rt <- as.matrix(rt) rownames(rt) <- rt[, 1] drug <- rt[,
2:ncol(rt)] dimnames <- list(rownames(drug),colnames(drug)) data <- matrix( as.numeric(
as.matrix(drug)),nrow=nrow(drug),dimnames=dimnames)
⾸先，读取前⾯保存的药物敏感性结果，设置相应的⾏名，并将其转换为矩阵形式。
•
•
•
mat<- impute.knn(data) drug<- mat$data drug<- avereps(drug)
考虑到药物敏感性数据中存在部分NA缺失值，通过impute.knn函数来评估并补齐药物数据。其中，impute.knn函数是⼀个使⽤最近邻平均来估算缺少的表达式数据的函数。
4.3 读取表达输⼊⽂件
•
•
更好的拼音
•
exp<- read.table( "", p= "\t", header=T, row.names = 1, check.names=F) dim( exp) exp[ 1: 4, 1: 4]
同时，读取整理完成的NCI-60细胞系中基因表达情况。结果显⽰：其中包含了60种不同肿瘤细胞系，23805个基因的表达情况。
4.4 提取特定基因表达
•
•
•
gene <- read.table( "",p= "\t",header=F,check.names=F) genelist <- as.vector(gene[, 1]) genelist
将提前准备的⽬标基因列表进⾏读取；结果显⽰，包括FANCD2，BRCA1，ABCC1，TP53，EGFR基因。
•
•
•
裙子设计图genelist <- gsub( " ", "",genelist) genelist <- interct(genelist,row.names( exp)) exp<- exp[genelist,]
结果显⽰：
4.5 药物敏感性计算
•
就像美丽蝴蝶飞
outTab< -data.frame
⾸先，新建⼀个空的数据框，⽤于保存后续分析结果。
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
for(Gene inrow.names(exp)){ x <- as.numeric(exp[Gene,]) #对药物循环 for(Drug inrow.names(drug)){ y <-
as.numeric(drug[Drug,]) corT <- st(x,y,method= "pearson") cor <- corT$estimate pvalue <- corT$p. value if( pvalue < 0.01) { outVector <- cbind(Gene,Drug,cor,pvalue) outTab <- rbind(outTab,outVector) } } }
随后，使⽤for循环，分别计算每个基因表达与不同药物之间的Pearson相关系数。根据P值<0.01为界，对分析结果进⾏筛选，并将结果保存给变量outTab。结果显⽰：最终得到了63个相关性分析结果。
•
•
outTab <- outTab[order( as.numeric( as.vector(outTab$pvalue))),] write.table(outTab, file= "", p= "\t", row.names=F, quote=F)
最后，将相关性分析结果进⾏输出保存。
4.6 可视化
•
•
library(ggplot2) library(ggpubr)
在可视化分析之前，⾸先加载绘图使⽤的ggplot2包和ggpubr包。下⾯，我们提供两种可视化的⽅法。
⽅法⼀：散点图
•
•
•
•
•
plotList_1<- list corPlotNum<- 16 if(nrow(outTab)<corPlotNum){ corPlotNum= nrow(outTab) }
定义⼀个空的列表plotList_1⽤于保存输出结果。提取分析结果中最显著的前16个结果；当然，如果结果outTab的⾏数少于corPlotNum的话，则将其⾏数赋值给corPlotNum。
沙僧图片•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
for(i in1:corPlotNum){ Gene <- outTab[i, 1] Drug <- outTab[i, 2] x <- as.numeric(exp[Gene,]) y <- as.numeric(drug[Drug,]) cor <- sprintf( "%.03f", as.numeric(outTab[i, 3])) pvalue= 0 if( as.numeric(outTab[i, 4])< 0.001){ pvalue= "p<0.001" } el{ pvalue=paste0( "p=",sprintf( "%.03f", as.numeric(outTab[i, 4]))) } df1 <- as. data.frame(cbind(x,y)) p1=ggplot( data= df1, aes(x = x, y = y))+ geom_point(size= 1)+ stat_smooth(method= "lm",=FALSE, formula=y~x)+ labs(x= "Expression",y= "IC50",title = paste0(Gene, ", ",Drug),subtitle = paste0( "Cor=",cor, ", ",pvalue))+ theme(axis.ticks = element_blank, = element_ = element_blank)+ theme_bw plotList_1[[i]]=p1 }
使⽤ggplot函数结合for循环，逐个绘制散点图，进⾏可视化展⽰。
⽅法⼆：箱线图
•
•
•
•
•
plotList_2<- list corPlotNum<- 16 if(nrow(outTab)<corPlotNum){ corPlotNum= nrow(outTab) }
同样，定义⼀个空的列表plotList_2⽤于保存输出结果。
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
for(i in 1:corPlotNum){ Gene <- outTab[i, 1] Drug <- outTab[i, 2] x <- as.numeric(exp[Gene,]) y <-
as.numeric(drug[Drug,]) df1 <- as.data.frame(cbind(x,y)) colnames(df1)[ 2] <- "IC50" df1$group <- ifel(df1$x > median(df1$x), "high", "low") compaired <- list(c( "low", "high")) p1 <- ggboxplot(df1, x = "group", y = "IC50", fill = "group", palette = c( "#00AFBB", "#E7B800"), add = "jitter", size = 0.5, xlab = paste0( "The_expression_of_", Gene), ylab = paste0( "IC50_of_", Drug)) + stat_compare_means(comparisons = compaired, method = "st", #设置统
