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

摘要

从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与肿瘤的关系尚不明确，两种药物与肝细胞癌的关系均有待后续实验的验证。