王炸！生信联合药物分析新套路，用这个数据库，立马高级！

cellMiner数据库的使用及药物敏感性预测

大家好，我是阿琛。在36策中，我们学习了“引药生变”的相关内容，其中2个问题4个类型，概括了基础科研中变量药物的全部研究策略。在内容，里面提到，无论是搞分子，还是研究药物，最终最好的结局也是做诊断和治疗的靶点。当然，在生信研究中同样也不例外，无论是单基因研究，还是建模分析，整个文章的落脚点都可以往药物与治疗水平靠一靠，毫无疑问文章的内涵和吸引力都会提升一个档次。

CellMiner数据库，主要是通过国家癌症研究所癌症研究中心(NCI)所列出的60种癌细胞为基础而建立的。该数据库最初发表于2009年，后于2012年在Cancer Research杂志上进行了更新，题目为“CellMiner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set”。大家后期在使用该数据库记得应用相关文献。

NCI-60细胞系是目前使用最广泛的用于抗癌药物测试的癌细胞样本群。大家可以通过它查询到 NCI-60细胞系中已确认的22379个基因，以及20503个已分析的化合物的数据（包括多种已获美国食品和药物监督局批准的药物，以及临床试验中的药物分子）。

下面，我们来看下相关数据的下载和药物敏感性分析过程。

1.CellMiner数据库的使用

1. 进入CellMiner数据库主页（https://discover.nci.nih.gov/cellminer/home.do）；

2. 点击Download Data Sets，进入数据下载界面；在下载界面中，可以看到两个不同的板块，分别是原始数据Raw Data Set和经过处理后的数据Processed Data Set；

3. 在此，我们直接选择经过处理后的数据Processed Data Set，勾选RNA表达数据（RNA: RNA-seq）和药物数据（Compound activity: DTP NCI-60）；

4. 点击按钮Get Processed 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$`FDA status` %in% c("FDA approved", "Clinical trial"),]rt1 <- rt1[,-c(1, 3:6)]write.table(rt1, file = "drug.txt",sep = "\t",row.names = F,quote = F)

为了保证分析结果的可靠性，选取经过临床试验（Clinical trial）和FDA批准（FDA approved）的药物结果。当然，你也可以选择将所有的药物进行保留。最终，一共得到263个药物结果，并将其保存为txt文件用于后续分析。

3. 基因表达数据的准备

rt2 <- read_excel(path = "RNA__RNA_seq_composite_expression.xls", skip = 9)colnames(rt2) <- rt2[1,]rt2 <- rt2[-1,-c(2:6)]write.table(rt2, file = "geneExp.txt",sep = "\t",row.names = F,quote = F)

同样的，读取表达数据，并对其进行整理和保存，用于后续的分析。

4. 药物敏感性分析

rm(list = ls())

4.1 引用包

library(impute)library(limma)

首先，加载分析使用的R包，包括impute包和limma包。

4.2 读取药物输入文件

rt <- read.table("drug.txt",sep="\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$datadrug<- avereps(drug)

考虑到药物敏感性数据中存在部分NA缺失值，通过impute.knn()函数来评估并补齐药物数据。其中，impute.knn()函数是一个使用最近邻平均来估算缺少的表达式数据的函数。

4.3 读取表达输入文件

exp <- read.table("geneExp.txt", sep="\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("gene.txt",sep="\t",header=F,check.names=F)genelist <- as.vector(gene[,1])genelist

将提前准备的目标基因列表进行读取；结果显示，包括FANCD2，BRCA1，ABCC1，TP53，EGFR基因。

genelist <- gsub(" ","",genelist)genelist <- intersect(genelist,row.names(exp))exp <- exp[genelist,]

结果显示：

4.5 药物敏感性计算

outTab <-data.frame()

首先，新建一个空的数据框，用于保存后续分析结果。

for(Gene in row.names(exp)){ x <- as.numeric(exp[Gene,])#对药物循环for(Drug in row.names(drug)){ y <- as.numeric(drug[Drug,]) corT <- cor.test(x,y,method="pearson") cor <- corT$estimate pvalue <- corT$p.valueif(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="drugCor.txt", sep="\t", row.names=F, quote=F)

最后，将相关性分析结果进行输出保存。

4.6 可视化

library(ggplot2)library(ggpubr)

在可视化分析之前，首先加载绘图使用的ggplot2包和ggpubr包。下面，我们提供两种可视化的方法。

方法一：散点图

plotList_1<- list()corPlotNum<- 16if(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=0if(as.numeric(outTab[i,4])<0.001){ pvalue="p<0.001" }else{ 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",se=FALSE, formula=y~x)+ labs(x="Expression",y="IC50",title = paste0(Gene,", ",Drug),subtitle = paste0("Cor=",cor,", ",pvalue))+ theme(axis.ticks = element_blank(), axis.text.y = element_blank(),axis.text.x = element_blank())+ theme_bw() plotList_1[[i]]=p1}

使用ggplot()函数结合for循环，逐个绘制散点图，进行可视化展示。

方法二：箱线图

plotList_2<- list()corPlotNum<- 16if(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 <- ifelse(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 = "wilcox.test", #设置统计方法 symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns"))) plotList[[i]]=p1}

通过for循环，结合使用ggpubr包中的ggboxplot()函数，逐个绘制箱线图，并使用stat_compare_means()函数进行两组间统计分析，从而进行可视化展示。

nrow<- ceiling(sqrt(corPlotNum))ncol<- ceiling(corPlotNum/nrow)ggarrange(plotlist=plotList_1,nrow=nrow,ncol=ncol)ggarrange(plotlist=plotList_2,nrow=nrow,ncol=ncol)

将绘制得到的结果按4x4的分布进行排列组合，并输出结果。

结果显示：