# PharmacoGx.R

# BiocManager::install("PharmacoGx")

require(PharmacoGx)

require(data.table)

path = "/Users/ftaj/OneDrive - University of Toronto/Drug_Response/"

dir.create(paste0(path, "Data/DRP_Training_Data"))

depmap_samples <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_Line_Info.csv"))

depmap_samples$stripped_cell_line_name

PharmacoGx::availablePSets()

# CTRPv2 ====

??summarizeSensitivityProfiles

require(PharmacoGx)

ctrp2 <- downloadPSet("CTRPv2_2015")

ctrp2_aac <- summarizeSensitivityProfiles(object = ctrp2, sensitivity.measure = "aac_recomputed", summary.stat = "median", fill.missing = T,)

fwrite(as.data.table(ctrp2@drug, keep.rownames = T), "Data/DRP_Training_Data/CTRP_DRUG_INFO.csv")

drug_info <- ctrp2@drug[, c("gene_symbol_of_protein_target", "target_or_activity_of_compound", "cpd_status")]

drug_info[1:5,]

drug_info <- as.data.table(drug_info, keep.rownames = T)

drug_info[rn %like% "Paclitaxel"]

drug_info[gene_symbol_of_protein_target != "" & cpd_status == "clinical"]$rn

# ctrp2@drug[1:20, c("drugid", "cpd_name")]

# ctrp2_ic50 <- PharmacoGx::summarizeSensitivityProfiles(object = ctrp2, sensitivity.measure = "ic50_recomputed", summary.stat = "median", fill.missing = T)

# Convert to data.table

ctrp2_aac <- as.data.table(ctrp2_aac, keep.rownames = T)

# ctrp2_ic50 <- as.data.table(ctrp2_ic50, keep.rownames = T)

dim(ctrp2_aac)

# dim(ctrp2_ic50)

# Convert to long format

ctrp2_aac <- melt.data.table(data = ctrp2_aac)

# ctrp2_ic50 <- melt.data.table(data = ctrp2_ic50)

# Remove NAs

ctrp2_aac <- ctrp2_aac[!is.na(value)]

# ctrp2_ic50 <- ctrp2_ic50[!is.na(value)]

colnames(ctrp2_aac) <- c("cpd_name", "ccl_name", "area_above_curve")

# colnames(ctrp2_ic50) <- c("cpd_name", "ccl_name", "ic50")

# Add smiles data

dim(ctrp2_aac)  ### 363,634 dose-response curves

dim(ctrp2_ic50)

dim(ctrp2)

ctrp2_smiles <- data.table(drugid = ctrp2@drug$drugid, cpd_smiles = ctrp2@drug$cpd_smiles)

sum(unique(ctrp2@drug$drugid) %in% unique(ctrp2_aac$cpd_name)) / length(unique(ctrp2_aac$cpd_name))

length(unique(ctrp2_aac$cpd_name))

length(unique(ctrp2_ic50$cpd_name))

length(unique(ctrp2_smiles$drugid))

sum(unique(ctrp2_smiles$drugid) %in% unique(ctrp2_aac$cpd_name)) / length(unique(ctrp2_aac$cpd_name))

sum(unique(ctrp2_smiles$drugid) == unique(ctrp2_ic50$cpd_name))

sum(tolower(unique(ctrp2_smiles$drugid)) == tolower(unique(ctrp2_aac$cpd_name)))

sum(tolower(unique(ctrp2_smiles$drugid)) == tolower(unique(ctrp2_ic50$cpd_name)))

changed = !(tolower(unique(ctrp2_smiles$drugid)) == tolower(unique(ctrp2_aac$cpd_name)))

changed = !tolower(unique(ctrp2_smiles$drugid)) == tolower(unique(ctrp2_ic50$cpd_name))

temp = data.table(before = unique(ctrp2_smiles$drugid)[changed], after = unique(ctrp2_aac$cpd_name)[changed])

temp = data.table(before = unique(ctrp2_smiles$drugid)[changed], after = unique(ctrp2_ic50$cpd_name)[changed])

View(temp)

ctrp2_smiles

ctrp2_aac

unique(ctrp2_smiles$drugid)[!(unique(ctrp2_smiles$drugid) == unique(ctrp2_aac$cpd_name))]

unique(ctrp2_aac$cpd_name)[!(unique(ctrp2_smiles$cpd_name) == unique(ctrp2_aac$cpd_name))]

drug_table <- data.table(canonical_name = rownames(drugInfo(ctrp2)),

                         common_name = drugInfo(ctrp2)[, c("cpd_name")],

                         cpd_smiles = drugInfo(ctrp2)[, c("cpd_smiles")])

final_ctrpv2_aac <- merge(ctrp2_aac, drug_table[, c("canonical_name", "cpd_smiles")], by.x = "cpd_name", by.y = "canonical_name")

final_ctrpv2_ic50 <- merge(ctrp2_ic50, drug_table[, c("canonical_name", "cpd_smiles")], by.x = "drugid", by.y = "canonical_name")

# Save

fwrite(final_ctrpv2_aac, paste0(path, "Data/DRP_Training_Data/CTRP_AAC_SMILES.txt"))

fwrite(final_ctrpv2_ic50, paste0(path, "Data/DRP_Training_Data/CTRP_IC50_SMILES.txt"))

final_ctrpv2_aac <- fread(paste0(path, "Data/DRP_Training_Data/CTRP_AAC_SMILES.txt"))

final_ctrpv2_ic50 <- fread(paste0(path, "Data/DRP_Training_Data/CTRP_IC50_SMILES.txt"))

# Add disease information

require(stringr)

line_info <- fread("Data/DRP_Training_Data/DepMap_20Q2_Line_Info.csv")

final_ctrpv2_aac <- fread("Data/DRP_Training_Data/CTRP_AAC_SMILES.txt")

# Remove all hyphens, convert to upper case

line_info$other_ccl_name <- str_replace_all(toupper(line_info$stripped_cell_line_name), "-", "")

final_ctrpv2_aac$other_ccl_name <- str_replace_all(toupper(final_ctrpv2_aac$ccl_name), "-", "")

# Remove all spaces

line_info$other_ccl_name <- str_replace_all(toupper(line_info$other_ccl_name), " ", "")

final_ctrpv2_aac$other_ccl_name <- str_replace_all(toupper(final_ctrpv2_aac$other_ccl_name), " ", "")

sum(unique(final_ctrpv2_aac$other_ccl_name) %in% unique(line_info$other_ccl_name)) / length(unique(final_ctrpv2_aac$other_ccl_name))  ### 0.89177!!!

# A quarter of the data cannot be paired...

# Find what cannot be paired

final_ctrpv2_aac[other_ccl_name %in% unique(final_ctrpv2_aac$other_ccl_name)[!(unique(final_ctrpv2_aac$other_ccl_name) %in% unique(line_info$other_ccl_name))]]

final_ctrpv2_aac <- merge(final_ctrpv2_aac, line_info[, c("other_ccl_name", "primary_disease")], by = "other_ccl_name")

final_ctrpv2_aac$other_ccl_name <- NULL

setcolorder(final_ctrpv2_aac, neworder = c("cpd_name", "ccl_name", "primary_disease", "area_above_curve", "cpd_smiles"))

fwrite(final_ctrpv2_aac, "Data/DRP_Training_Data/CTRP_AAC_SMILES.txt")

# ==== Update August 2021 =====

# Fix Tipifarnbi SMILES

require(data.table)

ctrp_aac <- fread("Data/DRP_Training_Data/CTRP_AAC_SMILES.txt")

ctrp_aac[cpd_name == "Tipifarnib"]$cpd_smiles <- "CN1C=NC=C1C(C2=CC=C(C=C2)Cl)(C3=CC4=C(C=C3)N(C(=O)C=C4C5=CC(=CC=C5)Cl)C)N"

fwrite(ctrp_aac, "Data/DRP_Training_Data/CTRP_AAC_SMILES.txt")

# ctrp_aac <- fread("Data/DRP_Training_Data/CTRP_AAC_SMILES.txt")

# ==== Update September 2021 ====

# Add canonical SMILES for some drugs

require(data.table)

gdsc1_aac <- fread("Data/DRP_Training_Data/GDSC1_AAC_SMILES.txt")

gdsc1_aac[cpd_name == "Tamoxifen"]$cpd_smiles <- "CCC(=C(C1=CC=CC=C1)C2=CC=C(C=C2)OCCN(C)C)C3=CC=CC=C3"

gdsc1_aac[cpd_name == "JW-7-52-1"]$cpd_smiles <- "CCC(=O)N1CCN(CC1)C2=C(C=C(C=C2)N3C(=O)C=CC4=CN=C5C=CC(=CC5=C43)C6=CC7=CC=CC=C7N=C6)C(F)(F)F"

gdsc1_aac[cpd_name == "KIN001-135"]$cpd_smiles <- "COC1=C(C=C2C(=C1)N=CN2C3=CC(=C(S3)C#N)OCC4=CC=CC=C4S(=O)(=O)C)OC"

fwrite(gdsc1_aac, "Data/DRP_Training_Data/GDSC1_AAC_SMILES.txt")

gdsc2_aac <- fread("Data/DRP_Training_Data/GDSC2_AAC_SMILES.txt")

gdsc2_aac[cpd_name == "Tamoxifen"]$cpd_smiles <- "CCC(=C(C1=CC=CC=C1)C2=CC=C(C=C2)OCCN(C)C)C3=CC=CC=C3"

fwrite(gdsc2_aac, "Data/DRP_Training_Data/GDSC2_AAC_SMILES.txt")

# CCLE ====

# ccle <- downloadPSet("CCLE")

# ccle_auc <- summarizeSensitivityProfiles(pSet = ccle, sensitivity.measure = "auc_recomputed", summary.stat = "median", fill.missing = T)

# # Convert to data.table

# ccle_auc <- as.data.table(ccle_auc, keep.rownames = T)

# dim(ccle_auc)

# # Convert to long format

# ccle_auc <- melt.data.table(data = ccle_auc)

# # Remove NAs

# ccle_auc <- ccle_auc[!is.na(value)]

# colnames(ccle_auc) <- c("cpd_name", "ccl_name", "area_under_curve")

# drug_table <- data.table(canonical_name = rownames(drugInfo(ccle)), common_name = drugInfo(ccle)[, c("drug.name")], pubchem = drugInfo(ccle)[, "PubCHEM"])

# ccle_smiles <- webchem::pc_prop(cid = drug_table$pubchem, properties = "CanonicalSMILES")

# # Remove drugs without pubchem entries

# sum(is.na(drug_table$pubchem))  # 15

# drug_table <- drug_table[!is.na(pubchem)]

# drug_table <- merge(drug_table, ccle_smiles, by.x = "pubchem", by.y = "CID")

# 

# 

# ccle_final <- merge(ccle_auc, drug_table[, c("canonical_name", "CanonicalSMILES")], by.x = "cpd_name", by.y = "canonical_name")

# # Save

# fwrite(ccle_final, paste0(path, "Data/DRP_Training_Data/CCLE_AUC_SMILES.txt"))

# GDSC1 ====

gdsc2 <- downloadPSet("GDSC")

gdsc2_auc <- summarizeSensitivityProfiles(pSet = gdsc2, sensitivity.measure = "auc_recomputed", summary.stat = "median", fill.missing = T)

# Convert to data.table

gdsc2_auc <- as.data.table(gdsc2_auc, keep.rownames = T)

dim(gdsc2_auc)

# Convert to long format

gdsc2_auc <- melt.data.table(data = gdsc2_auc)

# Remove NAs

gdsc2_auc <- gdsc2_auc[!is.na(value)]

colnames(gdsc2_auc) <- c("cpd_name", "ccl_name", "area_under_curve")

drug_table <- data.table(canonical_name = rownames(drugInfo(gdsc2)), common_name = drugInfo(gdsc2)[, c("drug.name")], pubchem = drugInfo(gdsc2)[, "PubCHEM"])

gdsc2_smiles <- webchem::pc_prop(cid = drug_table$pubchem, properties = "CanonicalSMILES")

# Remove drugs without pubchem entries

sum(is.na(drug_table$pubchem))  # 15

drug_table <- drug_table[!is.na(pubchem)]

drug_table <- merge(drug_table, gdsc2_smiles, by.x = "pubchem", by.y = "CID")

gdsc2_final <- merge(gdsc2_auc, drug_table[, c("canonical_name", "CanonicalSMILES")], by.x = "cpd_name", by.y = "canonical_name")

# Save

fwrite(gdsc2_final, paste0(path, "Data/DRP_Training_Data/GDSC2_AUC_SMILES.txt"))

gdsc2_final <- fread(paste0(path, "Data/DRP_Training_Data/GDSC2_AUC_SMILES.txt"))

colnames(gdsc2_final)[4] <- "cpd_smiles"