cellMiner药物敏感性分析

cellMiner药物敏感性分析

1.数据下载

https://discover.nci.nih.gov/cellminer/home.do

下载并解压，放在工作目录下

2.药物数据整理

rm(list=ls())
library(stringr)
library(rio)
drug0 = rio::import("DTP_NCI60_ZSCORE.xlsx", skip = 8)
drug0 = drug0[,-c(67,68)]
k1 = drug0$`Drug name`!="-";table(k1)

## k1
## FALSE  TRUE 
##  3766 21087

table(drug0$`FDA status`)

## 
##              - Clinical trial   FDA approved 
##          23993            546            314

k2 = drug0$`FDA status`!="-";table(k2)

## k2
## FALSE  TRUE 
## 23993   860

drug0 = drug0[k1&k2,]
drug = apply(drug0[,-(1:6)],2,as.numeric)
rownames(drug)= drug0$`Drug name`

drug[1:4,1:4]

##                  BR:MCF7 BR:MDA-MB-231 BR:HS 578T BR:BT-549
## METHOTREXATE        0.70         -1.22      -1.89     -0.88
## 6-THIOGUANINE       0.48         -0.31      -1.09     -0.44
## 6-MERCAPTOPURINE    0.70         -0.44      -0.55     -1.44
## Colchicine          0.32            NA         NA        NA

# 缺失值处理

library(impute)
library(limma)

drug = impute.knn(drug)$data #报错

## Error in impute.knn(drug): a column has more than 80 % missing values!

a = apply(drug, 2, function(x){sum(is.na(x))/length(x)})
tail(sort(a),10) # 确实有超过80%NA 的列，要删掉

##    RE:A498  RE:CAKI-1 BR:HS 578T   LE:K-562  LC:HOP-92  BR:BT-549   BR:T-47D 
## 0.04651163 0.04883721 0.05348837 0.05465116 0.05697674 0.06046512 0.06860465 
##      LE:SR    LC:EKVX   ME:MDA-N 
## 0.07906977 0.19186047 0.84534884

drug = drug[,-which.max(a)]
drug = impute.knn(drug)$data 
drug = avereps(drug)

3.表达矩阵整理

exp0 = rio::import("RNA__RNA_seq_composite_expression.xls",skip = 10)
exp = as.matrix(exp0[,-(1:6)])
rownames(exp) = exp0$`Gene name d`
exp = avereps(exp)
exp = exp[,-which.max(a)] #上面删掉了一列，这里也必须删掉
identical(colnames(exp),colnames(drug))

## [1] TRUE

5.药敏分析

其实本质上就是基因表达量与药物IC50值的相关性分析

tinyarray 里的cor.one函数支持一列与所有列的相关性分析

g = "EGFR"
dat = t(rbind(exp[g,],drug))
colnames(dat)[1] = g
dat[1:4,1:4]

##                EGFR METHOTREXATE 6-THIOGUANINE 6-MERCAPTOPURINE
## BR:MCF7       0.242         0.70          0.48             0.70
## BR:MDA-MB-231 3.913        -1.22         -0.31            -0.44
## BR:HS 578T    2.428        -1.89         -1.09            -0.55
## BR:BT-549     3.215        -0.88         -0.44            -1.44

library(tinyarray)
re = cor.one(dat,g)
k = re$p.value<0.05 & abs(re$cor)>0.5;table(k)

## k
## FALSE  TRUE 
##   752    15

remini = re[k,]
remini

##                    p.value        cor obsnumber
## Noscapine     1.740620e-05 -0.5278020        59
## Tamoxifen     1.920620e-06 -0.5748512        59
## ciclosporin   2.433162e-05 -0.5199893        59
## auranofin     1.013795e-07 -0.6280054        59
## Staurosporine 6.910490e-06  0.5483959        59
## Dasatinib     5.932590e-07  0.5972467        59
## XAV-939       3.766265e-05  0.5095013        59
## Sapitinib     1.906967e-05  0.5256921        59
## Pipamperone   1.515924e-05 -0.5309705        59
## BMS-690514    3.847482e-06  0.5607794        59
## Bafetinib     3.823063e-05 -0.5091358        59
## spebrutinib   8.142489e-07  0.5913741        59
## S-63845       5.427220e-06 -0.5535653        59
## AZD-5991      5.042577e-05 -0.5022987        59
## BLU-667       5.958846e-07  0.5971656        59

可以根据相关系数与p值筛选药物-基因对啦

6.可视化

相关性点图

library(ggpubr)
pdat = data.frame(dat[,c("Tamoxifen","EGFR")],
                  check.names = F)
colnames(pdat)

## [1] "Tamoxifen" "EGFR"

sp1 <- ggscatter(pdat, x = "Tamoxifen", y = "EGFR",
                 add = "reg.line",  # Add regressin line 
                 add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
                 conf.int = TRUE # Add confidence interval
) + stat_cor(method = "pearson")
sp1

箱线图

pdat$group = ifelse(pdat$EGFR>median(pdat$EGFR),"high","low")

ggboxplot(pdat,"group","Tamoxifen",fill = "group")+
  stat_compare_means(comparisons = list(c("high","low")),
                     aes(label = after_stat(p.signif)))+
  scale_fill_manual(values = c("#2874C5", "#f87669"))