cellMiner数据库的使用及药物敏感性预测—科研工具箱

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$data
drug <- 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.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=\"drugCor.txt\", sep=\"\\t\", row.names=F, quote=F)

最后，将相关性分析结果进行输出保存，后续进行可视化