Bliss协同矩阵分析
步骤：
- 通过高内涵或者酶标仪的下机数据，提取为矩阵读值，并标出相应浓度
- 计算生成不同药物浓度组合协同矩阵的抑制率文件，有标准格式
- 转成长数据，并保存为“ld-”+同名的.xlsx文件
- 利用synergyfinder R包，计算协同分数并绘制出网站类似的协同分析图
- 同时画出热图，关于抑制率矩阵和协同分数矩阵
1. 通过高内涵或者酶标仪的下机数据，提取为矩阵读值，并标出相应浓度
2. 计算生成不同药物浓度组合协同矩阵的抑制率文件，有标准格式
药物1，药物2，浓度单位
接下来是矩阵信息，A4为空表头，A列为药物1的浓度梯度，第4行为药物2的浓度梯度
抑制率矩阵文件，标准格式
数据如下：
Drug1: Ribavirin
Drug2: Mycophenolate mofetil
ConcUnit: μM
0.3 0.1 0.037 0
30.00 99.85766446 96.62107804 82.69385482 58.77838975
10.00 96.87171236 67.22260041 30.10706108 9.635497246
3.33 85.55603688 38.22946965 12.53481032 0.355838851
1.11 78.10508076 29.32730986 12.27489325 -1.983414815
0.37 74.41054521 23.84120304 13.02679621 -0.07116777
0.00 71.89491924 21.61643666 5.891453679 0
将多个该文件放在一个.xlsx文件中，含多个工作表
但是尽量是细胞病毒一样的才放在一起，因为暂时没有额外定义细胞病毒参数，只能靠文件名中明确细胞和病毒类型
3. 转成长数据，并保存为“ld-”+同名的.xlsx文件
多个工作表对应多个工作表，如下为其中一个
转换后的长数据
R代码如下：文件名为：source("./useSynergy_20240903.R")
library(synergyfinder)
library(tidyverse)
library(reshape2)
library(readxl)
library(plotly)
library(grid)
library(cowplot)
#---------------------函数-start------------------------------------------------------------
#矩阵宽数据转长数据
width2long <- function(matrix,sheet_name){
require(synergyfinder)
require(plotly)
require(cowplot)
require(tidyverse)
require(reshape2)
df <- readxl::read_xlsx(matrix,col_names = FALSE,sheet = sheet_name)
cr = dim(df)
dft <- df[4:cr[1],1:cr[2]]
dft[1,1] <- "conc1"
dg2 <- dft[1,]
colnames(dft) <- dg2
dft <- dft[2:dim(dft)[1],]
dft_long <- melt(dft,id.vars = c("conc1"),
measure.vars=dg2[2:cr[2]],
variable.name=c("conc2"),
value.name = "response")
dft_long$block_id =1
dft_long$cell_line_name = "-"
dft_long$drug1=rep(df[1,2],length(dft_long$conc1))
dft_long$drug2=rep(df[2,2],length(dft_long$conc1))
dft_long$conc_unit="μM"
dft_long$response <- as.numeric(dft_long$response)
dft_long$conc1 <- as.double(dft_long$conc1)
dft_long$conc2 <- as.double(as.character(dft_long$conc2))
dft_long$block_id <- as.numeric(dft_long$block_id)
return(dft_long)
}
4. 利用synergyfinder R包，计算协同分数并绘制出网站类似的协同分析图
继续定义绘图函数one_step_draw <- function(longd,cell,virus)
该函数需要长数据，细胞，病毒参数，
R代码如下：同样在文件名为：source("./useSynergy_20240903.R")
# -----------------一键绘图（长数据）-----------------------
#传入长数据，一键出等高线图，单药抑制图，旧版热图
#longd代表整理好的长数据
#cell代表所使用的细胞系，体现在文件名中
#virus代表所使用的毒株，体现在文件名中
one_step_draw <- function(longd,cell,virus){
require(synergyfinder)
require(plotly)
require(cowplot)
require(tidyverse)
df <- longd
df$response <- round(df$response,2)
res <- ReshapeData(
data = df,
data_type = "inhibition",
impute = TRUE,
impute_method = NULL,
noise = TRUE,
seed = 1)
res <- CalculateSynergy(
data = res,
method = c("ZIP", "HSA", "Bliss", "Loewe"),
Emin = NA,
Emax = NA,
correct_baseline = "non")
res <- CalculateSensitivity(
data = res,
correct_baseline = "non")
drug_title = paste0(res$drug_pairs$drug1,"-",res$drug_pairs$drug2,"_",cell,"_",virus,
res$drug_pairs$drug1,"_",max(res$response$conc1),"-",
res$drug_pairs$drug2,"_",max(res$response$conc2))
px = Plot2DrugContour(
data = res,
plot_block = 1,
drugs = c(1, 2),
plot_value = "Bliss_synergy",
dynamic = TRUE,
plot_title = paste0(drug_title,"\n","Bliss_Synery Score",":",round(res$drug_pairs$Bliss_synergy,2)),
text_size_scale = 1.5)
# plotly::save_image(px,paste0(drug_title,"_contour",".png"),width = 800,height = 800)
# PlotDoseResponse(data = res,block_ids = c(1),drugs = c(1,2),save_file = TRUE,file_type = "png",
# file_name = drug_title)
PlotDoseResponse(
data = res,
block_ids = c(1),
drugs = c(1,2),
save_file = TRUE,
file_type = "png",
file_name = paste0(drug_title,"_curve"))
bs <- res$synergy_scores$Bliss_synergy %>% as.data.frame()
bs$cocn1 <- df$conc1
bs$conc2 <- df$conc2
colnames(bs) <- c("value","conc1","conc2")
bs <- filter(bs,bs$conc1 !=0 & bs$conc2 != 0)
bs$conc1 <- as.character.numeric_version(bs$conc1)
bs$conc2 <- as.character.numeric_version(bs$conc2)
bs$value <- round(bs$value,1)
p <- ggplot(bs,aes(conc1,conc2))+
geom_tile(aes(fill=value),colour = "black", size = 1)+
scale_fill_gradient2(low = "#5C5DAF",mid = "white",high = "#EA2E2D",limits = c(-max(abs(bs$value)), max(abs(bs$value))))+
geom_text(aes(label=value),col ="black",size=10,parse = T)+ #框内字体
theme_minimal()+
theme(
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),# 去掉x 轴
axis.text = element_text(size = 32,face = "bold",color = "black"),
legend.key.height = unit(1.37,"in"),
legend.key.width = unit(0.8,"in"),
legend.title = element_text(size = 26,face = "bold",color = "black",vjust = 2),
legend.text = element_text(size = 32),
plot.title = element_text(vjust = 20,size = 18,face = "bold",color = "black"),
axis.title = element_text(size = 28,face = "bold",color = "black",vjust = 6),
axis.title.x = element_text(vjust = -4),
axis.title.y = element_text(vjust = 6))+
labs(title = drug_title,fill="Synery \nScore",x =paste0(df$drug1[1]," (μM)"),y=paste0(df$drug2[2]," (μM)"))
ggsave(paste0(drug_title,"_heatmap",".png"), egg::set_panel_size(p, width=unit(8, "in"), height=unit(8, "in")),
width = 14, height = 14, units = 'in', dpi = 300)
dft <- df
dft$conc2 <- round(df$conc2,2)
dft$conc1 <- as.character.numeric_version(dft$conc1)
dft$conc2 <- as.character.numeric_version(dft$conc2)
dft$response <- round(dft$response,1)
p1 <- ggplot(dft,aes(conc1,conc2))+geom_tile(aes(fill=response))+
scale_fill_gradient2(low = "#04ff04",mid = "#ffffff",high = "#ff0000",limits = c(-max(df$response)-10, max(df$response)+10))+
geom_text(aes(label=response),col ="black",size=10,parse = T)+
theme_minimal()+
labs(x =paste0(df$drug1[1]," (μM)"),y=paste0(df$drug2[2]," (μM)"),fill="Inhibition (%)",
title = paste0(drug_title,"\n","\n","\n"," Dose Response Matrix"))+
theme(plot.title = element_text(vjust = 6,size = 38,face = "bold",color = "black"),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
axis.title = element_text(size = 28,face = "bold",color = "black"),
axis.title.x = element_text(vjust = -3),
axis.title.y = element_text(vjust = 6),
axis.text = element_text(size = 26,face = "bold",color = "black"),
legend.key.height = unit(1.46,"in"),
legend.key.width = unit(0.8,"in"),
legend.title = element_text(size = 26,face = "bold",color = "black",vjust = 2),
legend.text = element_text(size = 32), plot.margin = margin(60,60,60,60));p
ggsave(paste0(drug_title,"_response",".png"), egg::set_panel_size(p1, width=unit(8, "in"), height=unit(8, "in")),
width = 14, height = 14, units = 'in', dpi = 300)
}
如果只画协同矩阵分数图，可以参考以下函数：
synery_heatmap <- function(longd,cell,virus)
# 只有协同分数矩阵图
#传入长数据，出新的协同抑制热图
synery_heatmap <- function(longd,cell,virus){
require(synergyfinder)
require(plotly)
require(cowplot)
require(tidyverse)
#dft <- width2long(matrix = martix)
dft <- longd
res <- ReshapeData(
data = dft,
data_type = "inhibition",
impute = TRUE,
impute_method = NULL,
noise = TRUE,
seed = 1)
res <- CalculateSynergy(
data = res,
method = c("ZIP", "HSA", "Bliss", "Loewe"),
Emin = NA,
Emax = NA,
correct_baseline = "non")
res <- CalculateSensitivity(
data = res,
correct_baseline = "non"
)
drug_title = paste0(res$drug_pairs$drug1,"-",res$drug_pairs$drug2,"_",cell,"_",virus,"_",
res$drug_pairs$drug1,"_",max(res$response$conc1),"-",
res$drug_pairs$drug2,"_",max(res$response$conc2))
bs <- res$synergy_scores$Bliss_synergy %>% as.data.frame()
bs$cocn1 <- dft$conc1
bs$conc2 <- dft$conc2
colnames(bs) <- c("value","conc1","conc2")
bs <- filter(bs,bs$conc1 !=0 & bs$conc2 != 0)
bs$conc1 <- as.character.numeric_version(bs$conc1)
bs$conc2 <- as.character.numeric_version(bs$conc2)
bs$value <- round(bs$value,2)
p <- ggplot(bs,aes(conc1,conc2))+
geom_tile(aes(fill=value),colour = "black", size = 1)+
scale_fill_gradient2(low = "#5C5DAF",mid = "white",high = "#EA2E2D",limits = c(-max(bs$value), max(bs$value)))+
geom_text(aes(label=value),col ="black",size=12,parse = T)+
theme_minimal()+
theme(axis.title.x=element_blank(), # 去掉 title
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),# 去掉x 轴
axis.title.y=element_blank(),
axis.text = element_text(size = 32,face = "bold",color = "black"),
legend.key.height = unit(1.46,"in"),
legend.title = element_blank(),
legend.text = element_text(size = 32))+
scale_x_discrete(position = "top")+
labs(title = drug_title)
fn = paste0(drug_title,".jpg")
ggsave(fn,p,width =11,height = 8.55,units = "in")
}
