药物预测之差异基因为什么不行
有了前面的教程：药物预测之认识表达量矩阵和药物IC50 的背景知识铺垫，认识了Cancer Therapeutics Response Portal (CTRP) 和 Genomics of Drug Sensitivity in Cancer (GDSC) 两个数据库资源。也介绍了2021年7月新鲜出炉的 药物预测R包之oncoPredict的用法以及一个超级老旧的R包的用法，见：药物预测R包之pRRophetic
现在可以尝试一下理解药物预测的原理啦，首先呢在前面的教程 药物预测之相关性为什么不行，我们发现简简单单的表达量相关性居然都可以勉勉强强得到还算是合理的结果啊！现在让我们一起看看差异基因能不能进行药物预测!
本次仍然是使用Genomics of Drug Sensitivity in Cancer (GDSC) 这个数据库资源里面的Bortezomib这个药物的表达量矩阵和药物信息，参考前面的教程 药物预测之相关性为什么不行，代码是：
rm(list = ls())
library(pRRophetic)
if(!file.exists('Bortezomib_in_cgp2016.Rdata')){
# pRRophetic的安装方式 见 ；https://mp.weixin.qq.com/s/avETGVE8_5I5CrvljG2Tqg
library(pRRophetic)
# 都是 pRRophetic 包的内置数据集，来源于 公共数据库哦 ！
data(cgp2016ExprRma)
dim(cgp2016ExprRma)
data(PANCANCER_IC_Tue_Aug_9_15_28_57_2016)
length(unique(drugData2016$Drug.name))
unique(drugData2016$Drug.name)
cgp2016ExprRma[1:4,1:4]
drugData2016[1:4,1:4]
Bort_drugs=drugData2016[drugData2016$Drug.name=='Bortezomib',]
dim(Bort_drugs)
kp = Bort_drugs$Cell.line.name %in% colnames(cgp2016ExprRma)
table(kp)
Bort_drugs=Bort_drugs[kp,]
pos=match( Bort_drugs$Cell.line.name ,colnames(cgp2016ExprRma) )
Bort_ExprRma=cgp2016ExprRma[,pos]
identical(colnames(Bort_ExprRma),Bort_drugs$Cell.line.name)
save(Bort_ExprRma,Bort_drugs,file = 'Bortezomib_in_cgp2016.Rdata')
}
load(file = 'Bortezomib_in_cgp2016.Rdata')
identical(colnames(Bort_ExprRma),Bort_drugs$Cell.line.name)
获取如下所示的表达量矩阵和药物信息：
> Bort_ExprRma[1:4,1:4]
MC-CAR ES3 ES5 ES7
TSPAN6 3.238273 8.690198 8.233101 8.333466
TNMD 2.982254 3.091473 2.824687 3.966757
DPM1 10.235491 9.992487 10.015884 9.793991
SCYL3 4.856061 4.572198 4.749715 3.976923
> Bort_drugs[1:4,1:4]
Drug.name Drug.Id Cell.line.name Cosmic.sample.Id
14878 Bortezomib 104 MC-CAR 683665
14879 Bortezomib 104 ES3 684055
14880 Bortezomib 104 ES5 684057
14881 Bortezomib 104 ES7 684059
做多次差异分析
这个时候 Bortezomib这个药物的 IC50 可以把全部的细胞系进行高低分组 ：
> fivenum(Bort_drugs$IC50)
[1] -9.9718999 -6.7907169 -5.9342808 -4.8832931 -0.5923361
分组有多种策略，比如中位值划分，以及top10%和bottom10%的划分，我们这里把多种组合都做一下差异。参考：不同癌症内部按照estimate的两个打分值高低分组看蛋白编码基因表达量差异，代码如下所示：
load(file = 'Bortezomib_in_cgp2016.Rdata')
identical(colnames(Bort_ExprRma),Bort_drugs$Cell.line.name)
Bort_ExprRma[1:4,1:4]
Bort_drugs[1:4,1:4]
fivenum(Bort_drugs$IC50)
od=order(Bort_drugs$IC50,decreasing = T)
library(limma)
deg_by_limma <- function(exprSet, design, contrast.matrix){
# step 1线性模型拟合
fit <- lmFit(exprSet, design)
# step 2贝叶斯检验
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
# step 3
allDiffOut <- topTable(fit2, coef = 1, n=Inf)
nrDEG <- na.omit(allDiffOut)
# write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
return(nrDEG)}
deg_list <- lapply(1:5, function(i){
# i=1
res_ExprRma=Bort_ExprRma[,head(od,i*floor(length(od)/10))]
sens_ExprRma=Bort_ExprRma[,tail(od,i*floor(length(od)/10))]
expr=cbind(res_ExprRma,sens_ExprRma)
group_list=c(rep('res',ncol(res_ExprRma)),
rep('sens',ncol(sens_ExprRma)))
options(digits = 4) # 设置全局的数字有效位数为4
# 指定任意两组进行差异分析
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(expr)
# contrast.matrix这个矩阵声明哪个分组去和哪个分组进行差异分析比较
contrast.matrix <- makeContrasts('sens-res',
levels = design)
allDiff <- deg_by_limma(expr, design, contrast.matrix)
return(allDiff)
})
这样我们就做了五次差异分析， 各种差异分析的上下调基因统计一下，代码如下：
logFC_cut = 1
padj_cut = 0.05
colnames(deg_list[[1]])
stat = lapply( deg_list , function(DEG){
DEG_sig = subset(DEG, abs(logFC)> logFC_cut &
adj.P.Val< padj_cut)
tb=table(DEG_sig$logFC>0)
up_total = as.numeric(tb["TRUE"])
down_total = as.numeric(tb["FALSE"])
return(c(up_total, down_total))
})
stat = as.data.frame(do.call(rbind,stat))
colnames(stat) = c("UP","DOWN")
rownames(stat) = names(deg_list)
stat$Type = rownames(stat)
stat
# stat_reshaped=reshape2::melt(stat)
# head(stat_reshaped)
stat_up = subset(stat, select = -DOWN)
stat_up$group = "UP"
stat_down = subset(stat, select = -UP)
stat_down$label = stat_down$DOWN
stat_down$DOWN = -1*stat_down$DOWN
stat_down$group = "DOWN"
p_multi_DEG=ggplot() +
geom_bar(data = stat_up, aes(x=Type, y=UP, fill=group),stat = "identity",position = 'dodge') +
geom_text(data = stat_up, aes(x=Type, y=UP, label=UP, vjust=-0.25)) +
geom_bar(data = stat_down, aes(x=Type, y=DOWN, fill=group),stat = "identity",position = 'dodge')+
geom_text(data = stat_down, aes(x=Type, y=DOWN, label=label, vjust= 1.25)) +
scale_fill_manual(values=c("#0072B5","#BC3C28")) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.8, hjust=0.8, size = 12)) +
xlab(label = "") + ylab(label = "DEG numbers") +
ggtitle("logFC cut-off : 1\nPadj cut-off : 0.05")
p_multi_DEG
结果如下所示：
可以看到，按照中位值把全部的样品划分为高低两个分组然后差异基因有统计学显著的数量就很少，目前看来，top30%的和bottom 30%的样品分组，得到的差异基因数量是最多的。
不过也并不是说，差异基因数量越多越好啦。这里仅仅是展现一个分析方法而已。
使用ssgsea方法来根据前面的上下调基因进行预测
首先获取需要预测的表达量矩阵；
library(GSVA)
library(limma)
library(GSEABase)
library(data.table)
data(bortezomibData)
dim(exprDataBortezomib)
exprDataBortezomib[1:4,1:4]
boxplot(exprDataBortezomib[,1:4])
# 264个样品的表达量矩阵
table(studyResponse)
然后把前面的上下调基因做出基因集格式：
colnames(deg_list[[1]])
upG = lapply( deg_list , function(DEG){
#DEG=deg_list[[1]]
DEG_sig = subset(DEG, abs(logFC)> logFC_cut &
adj.P.Val< padj_cut)
return(rownames(DEG_sig[DEG_sig$logFC>0,]))
})
names(upG)=paste0('upG_top_per_',1:5)
downG = lapply( deg_list , function(DEG){
#DEG=deg_list[[1]]
DEG_sig = subset(DEG, abs(logFC)> logFC_cut &
adj.P.Val< padj_cut)
return(rownames(DEG_sig[DEG_sig$logFC < 0,]))
})
names(downG)=paste0('downG_top_per_',1:5)
cg=c(upG,downG)
cg
as.matrix(exprDataBortezomib)[1:4,1:4]
geneSet <- GeneSetCollection(mapply(function(geneIds, keggId) {
GeneSet(geneIds, geneIdType=EntrezIdentifier(),
collectionType=KEGGCollection(keggId),
setName=keggId)
}, cg, names(cg)))
geneSet
然后运行ssgsea即可：
ssgseaScore=gsva(as.matrix(exprDataBortezomib), geneSet,
method='ssgsea', kcdf='Gaussian',
abs.ranking=TRUE)
ssgseaScore=as.data.frame(t(ssgseaScore))
head(ssgseaScore)
library(pheatmap)
pheatmap(ssgseaScore)
可以看到上调基因集在几乎全部的样品里面都是正向富集，而下调基因集在几乎全部的样品里面都是负向富集 ：
理论上 ssgsea 算法应该是无需考虑不同表达量矩阵批次效应的， 前面的数据库训练得到的药物敏感组相当于药物耐受组的上调基因，居然在我们的待预测表达量矩阵里面全部是正向富集，意味着我们的待预测样品基本上全部的药物敏感！
这个与实际情况不符合，待预测表达量矩阵本来就有药物信息，如下所示：
> table(studyResponse)
studyResponse
PGx_Responder = IE PGx_Responder = NR PGx_Responder = R
25 126 113
我们任意选择一个差异基因列表的ssgsea的打分划分高低分组，看看是否跟待预测表达量矩阵本来就有的药物信息匹配：
studyResponse FALSE TRUE
PGx_Responder = IE 13 12
PGx_Responder = NR 67 59
PGx_Responder = R 52 61
可以看到，完全没有任何预测价值，ssgsea的打分高低与否跟PGx_Responder不匹配。
对全部的10个基因集进行预测，代码如下所示：
t(
apply(ssgseaScore, 2, function(x){
# x=ssgseaScore[,1]
table(studyResponse,x>median(x) )
})
)
结果如下所示：
[,1] [,2] [,3] [,4] [,5] [,6]
upG_top_per_1 13 67 52 12 59 61
upG_top_per_2 16 67 49 9 59 64
upG_top_per_3 15 65 52 10 61 61
upG_top_per_4 11 63 58 14 63 55
upG_top_per_5 11 63 58 14 63 55
downG_top_per_1 12 69 51 13 57 62
downG_top_per_2 17 64 51 8 62 62
downG_top_per_3 16 65 51 9 61 62
downG_top_per_4 16 65 51 9 61 62
downG_top_per_5 15 66 51 10 60 62
可以看到，全部的10个基因集都完全没有任何预测价值，ssgsea的打分高低与否跟PGx_Responder不匹配。
结果很有意思
无论我们使用什么样的阈值分组来对Bortezomib这个药物划分敏感和耐受，进行表达量差异分析，得到了统计学显著的上下调基因列表，去看我们待预测的样品里面是否倾向于包含这样的基因列表，都没有预测价值！那我们根据药物敏感程度分组后的差异基因一般来说会作为该药物的特征基因集，如果它没有太大的意义，前面几万篇文章都错了吗？
同样的我们在一个随机群体里面 预测一个人的财富情况，选取穷人和富人进行差异分析，发现富人绝大部分都是北京人，从政，年龄大的，而穷人都是东北的，无业的，20岁。有了这样的背景知识，接下来我们开始随机预测一个人，发现它是30岁附近，无业，也是东北人，但是它并不穷啊，因为它姓王，王八的王。
那么，我们到底该如何去根据表达量信息预测其药物敏感程度呢？这个问题，跟如何根据身高体重籍贯工作等属性预测一个人的财富情况类似，都很难啊！