基于高通量数据的药物重定位

利用芯片数据分析疾病的潜在药物(旧药新用)

摘要

从GEO数据库下载两套与肝细胞疾病相关的数据(这里选择了两份肝细胞癌的数据),分别对数据进行差异表达分析,绘制HeatMap图。将差异表达的探针信息分为上调和下调组,分别提交给Connectivity Map数据库,对结果进行筛选,挑选出最佳的药物模型。

1.数据处理

下载肝细胞疾病相关的数据,分别对应NCBI GEO数据库中的GSE98383和GSE84402,其中,GSE98383的数据存在多组分类,这里我们选择前35个样本作为实验对象。
为描述方便,以下数据处理过程仅描述GSE84402。

1.1 载入数据

setwd("E:/生信/大四/技能训练/4yxq/GSE84402_data")
geneCELs=list.celfiles('E:/生信/大四/技能训练/4yxq/GSE84402_RAW/',full.name=T)
data.raw<-ReadAffy(filenames=geneCELs)

1.2 数据的标准化

使用rma算法对数据进行标准化处理

data.raw.rma <- rma(data.raw)

1.3 去重复探针

calls <- mas5calls(data.raw) # 获得探针的PMA信息
calls <- exprs(calls) # 提取包含探针的PMA信息的表达矩阵
head(calls)
present <- rowSums(calls == "P") # 统计每个探针在所有样本中分类为"P"的个数
present <- which(present != 0) # 提取至少在一个样本中表达的探针
rmaFiltered <- data.raw.rma[present,] # 留下至少在一个样本中表达的探针

1.4 去无对应基因的探针

library(hgu133plus2.db)
columns(hgu133plus2.db) #查看芯片注释包提供的注释信息
probe2symbol <- toTable(hgu133plus2SYMBOL) #提取注释信息+
exprSet <- as.data.frame(exprs(rmaFiltered)) #提取表达矩阵,转换为data.frame对象
exprSet$probe_id <- rownames(exprSet) #将行名转化为probe_id
probeset <- rownames(exprSet) #提取行名为探针名
Symbol <- as.character(as.list(hgu133plus2SYMBOL[probeset])) #提取探针对应的GENESYMBOL
exprSet_symbol <- cbind(Symbol,exprSet)
exprSet_symbol <- aggregate(x = exprSet_symbol[,-1],by = list(exprSet_symbol$Symbol), FUN = median)
exprSet_symbol <- na.omit(exprSet_symbol)
rownames(exprSet_symbol) <- exprSet_symbol$probe_id
exprSet_symbol$probe_id <- NULL
dim(exprSet_symbol)
#[1] 10769 29
exprSet_symbol_1 <- exprSet_symbol[,-1]

2.差异表达分析

2.1 t检验并做BH校正

HCC<-exprSet_symbol_1[,c(seq(1,by=2,length=14))]
NON<-exprSet_symbol_1[,c(seq(2,by=2,length=14))]
p_value<-c()
for(i in 1:nrow(exprSet_symbol_1)){
p_value[i]<-t.test(HCC[i,],NON[i,])$p.value
}
p_adjust <- p.adjust(p_value,method ="BH" ,n=length(p_value))
exprSet_symbol_1$p_value<-p_value
exprSet_symbol_p_value_0001 <- exprSet_symbol_1[p_value<0.001,]
dim(exprSet_symbol_p_value_0001)
#[1] 1062 30
exprSet_symbol_1$p_adjust<-p_adjust
exprSet_symbol_p_adjust_0001 <- exprSet_symbol_1[p_adjust<0.001,]
dim(exprSet_symbol_p_adjust_0001)
#[1] 388 30

2.2利用log2 fold change确定上下调探针

HCC_p_adjust_0001<-exprSet_symbol_p_adjust_0001[,c(seq(1,by=2,length=14))]
NON_p_adjust_0001<-exprSet_symbol_p_adjust_0001[,c(seq(2,by=2,length=14))]
HCC_avg<-apply(HCC_p_adjust_0001,1,mean)
NON_avg<-apply(NON_p_adjust_0001,1,mean)
exprSet_symbol_p_adjust_0001$HCC_avg<-HCC_avg
exprSet_symbol_p_adjust_0001$NON_avg<-NON_avg
exprSet_symbol_adjust_fold_up <- exprSet_symbol_p_adjust_0001[log2(HCC_avg/NON_avg)>0,]
exprSet_symbol_adjust_fold_down <- exprSet_symbol_p_adjust_0001[log2(HCC_avg/NON_avg)<0,]

2.3去除Connectivity Map数据库中没有的探针

ad_no<-read.table("adjust_notin.txt",header = T)
exprSet_symbol_adjust_fold_up$gene<-rownames(exprSet_symbol_adjust_fold_up)
adjust_fold_up<-merge(exprSet_symbol_adjust_fold_up, ad_no, by = 'gene',all.x = TRUE, sort = TRUE)
write.csv(adjust_fold_up, file="E:/生信/大四/技能训练/4yxq/GSE84402_data/adjust_up.csv")
exprSet_symbol_adjust_fold_down$gene<-rownames(exprSet_symbol_adjust_fold_down)
adjust_fold_down<-merge(exprSet_symbol_adjust_fold_down, ad_no, by = 'gene',all.x = TRUE, sort = TRUE)
write.csv(adjust_fold_down, file="E:/生信/大四/技能训练/4yxq/GSE84402_data/adjust_down.csv")

2.4 相关基因表达情况

TRIM21 SEL1L IBSP FAZH
四个基因结果均不存在显著差异表达

3.绘制HeatMap图

treats<-c(rep(c(“HCC”,”normal”),time=14))
exprSet_symbol_p_adjust_0001<-read.csv(“result_p_adjust_0001.csv”)
exprSet_symbol_p_adjust_0001<-as.matrix(exprSet_symbol_p_adjust_0001)
exprSet_symbol_p_adjust_a<-exprSet_symbol_p_adjust_0001[,1:28]
library(ggplot2)
library(pheatmap) #使用pheatmap绘图需要引入(Drawing with pheatmap needs to be introduced)
heatmap<-function(data,mainname,path){
pdf(path,onefile = FALSE,width=10, height=10) #保存地址(save PATH),onefile删除了第一页的空白
annotation_col = data.frame(Type = factor(treats))
rownames(annotation_col) = paste(colnames(data)) #构建annotation_col
#绘制heatmap(Drawing the heatmap)
pheatmap(data,
#cluster_col = FALSE,
scale=”row”, #按行聚类(Clustering by row)
treeheight_row=100,treeheight_col=20, #去除聚类树(Remove the clustering tree)
#加颜色(color)
color=colorRampPalette(c(‘green’,’black’,’red’))(50),
show_rownames=F, #去行名(Remove rowname)
annotation_col = annotation_col, #将轨道赋给heatmap(Assign orbits to heatmap)
main=mainname,
border_color=NA
)
dev.off()

}
heatmap(data=exprSet_symbol_p_adjust_a,mainname=”Heatmap Of t.test p_value<0.001”,path=”E:/生信/大四/技能训练/4yxq/GSE84402_data/heatmap.pdf”)

4.HeatMap结果

4.1 GSE84402结果

4.2 GSE98383结果

差异表达结果大致合格

5.相关药物预测

将差异表达基因传递给Connectivity Map数据库
按照p<0.05&mean从小到大对药物进行排序,在这里展示前五名药物结果

5.1 GSE84402结果

rank cmap name mean n enrichment p specificity percent non-null
18 DL-thiorphan -0.872 2 -0.977 0.00111 0 100
35 MS-275 -0.875 2 -0.956 0.00408 0.1235 100
82 quinostatin -0.823 2 -0.912 0.01539 0.1679 100
130 1,4-chrysenequinone -0.762 2 -0.874 0.03175 0.0781 100
47 triflusal -0.768 3 -0.86 0.00555 0.0154 100

5.2 GSE98383结果

rank cmap name mean n enrichment p specificity percent non-null
187 sanguinarine -0.76 2 -0.833 0.05531 0.0813 100
217 5248896 -0.73 2 -0.815 0.06786 0.0658 100
228 blebbistatin -0.747 2 -0.805 0.07487 0.1149 100
229 MS-275 -0.721 2 -0.804 0.07549 0.284 100
250 quinostatin -0.765 2 -0.789 0.08678 0.3435 100
quinostatin

5.3 药物分析

MS-275、quinostatin三种药物在两个试验集的结果中都排在Top5
MS-275(Entinostat)

与乳腺癌相关,目前同样被用于研究非小细胞肺癌(小细胞癌以外的癌症)
quinostatin

调节细胞生长的mTOR信号网络的小分子抑制剂

两种药物中,MS-275有被提到与非小细胞肺癌相关,quinostatin与肿瘤的关系尚不明确,两种药物与肝细胞癌的关系均有待后续实验的验证。