# Complete_Sample_Prep.R

# This script is intended to pair genomics, transcriptomics, proteomics and drug response data

# mainly from the DepMap resource.

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

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

require(data.table)

# ==== Cell line info cleanup ====

depmap_samples <- fread(paste0(path, "Data/DepMap/21Q2/sample_info.csv"))

# Subset relevant (machine learning) columns 

depmap_samples <- depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")]

fwrite(depmap_samples, paste0(path, "Data/DRP_Training_Data/DepMap_21Q2_Line_Info.csv"))

# depmap_samples <- fread("Data/DRP_Training_Data/DepMap_21Q2_Line_Info.csv")

# ==== Expression data cleanup ====

ccle_exp <- fread(paste0(path, "Data/DepMap/21Q2/CCLE_expression.csv"))

max(ccle_exp[, -1])

min(ccle_exp[, -1])

dim(ccle_exp)

ccle_exp[1:5, 1:20]

# Change column names to only contain HGNC name: replace everything after first word with ""

colnames(ccle_exp) <- gsub(" .+", "", colnames(ccle_exp))

colnames(ccle_exp)[1] <- "DepMap_ID"

# Merge with sample info to have cell line name in addition to DepMap ID

ccle_exp <- merge(ccle_exp, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")], by = "DepMap_ID")

ccle_exp[, DepMap_ID := NULL]

ccle_exp[1:5, 1:20]

# Move cell line name to the first column: just giving the column name to the function moves it to first place

setcolorder(ccle_exp, neworder = sort(colnames(ccle_exp)))

setcolorder(ccle_exp, neworder = c("stripped_cell_line_name", "primary_disease"))

ccle_exp[1:5, 1:20]

# Save

fwrite(ccle_exp, paste0(path, "Data/DRP_Training_Data/DepMap_21Q2_Expression.csv"), sep = ',')

ccle_exp <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_21Q2_Expression.csv"))

ccle_exp

# DIMENSIONS OF EXPRESSION DATA: 1375 X 19178

rm(ccle_exp)

# ==== Copy number data cleanup ====

ccle_cn <- fread(paste0(path, "Data/DepMap/21Q2/CCLE_gene_copy_number.csv"))

dim(ccle_cn)

ccle_cn[1:5, 1:10]

# Change column names to only contain HGNC name: replace everything after first word with ""

colnames(ccle_cn) <- gsub(" .+", "", colnames(ccle_cn))

colnames(ccle_cn)[1] <- "DepMap_ID"

# Merge with sample info to have cell line name in addition to DepMap ID

ccle_cn <- merge(ccle_cn, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")], by = "DepMap_ID")

ccle_cn[, DepMap_ID := NULL]

setcolorder(ccle_cn, neworder = sort(colnames(ccle_cn)))

setcolorder(ccle_cn, neworder = c("stripped_cell_line_name", "primary_disease"))

ccle_cn[1:5, 1:20]

which(is.na(ccle_cn), arr.ind = T)

ccle_cn[1407, 26569]

# Replace NA with 0

setnafill(ccle_cn, fill = 0, cols = unique(which(is.na(ccle_cn), arr.ind = T)[,2]))

dim(ccle_cn)

anyNA(ccle_cn)

sum(is.na(ccle_cn))

which(is.na(ccle_cn))

# Save

fwrite(ccle_cn, paste0(path, "Data/DRP_Training_Data/DepMap_21Q2_CopyNumber.csv"), sep = ',')

# DIMENSIONS OF COPY NUMBER DATA: 1740 X 27563

rm(ccle_cn)

gc()

# ==== Proteomic data cleanup ====

ccle_prot <- fread(paste0(path, "Data/DepMap/20Q2/CCLE_protein_quant_current_normalized.csv"))

dim(ccle_prot)

ccle_prot[1:5, 1:10]

ccle_prot[1:5, 48:60]

# Subset only the Uniprot accession (since its unique unlike HGNC) and the cell line experimental data

ccle_prot <- ccle_prot[, c(6, 49:ncol(ccle_prot)), with = F]

colnames(ccle_prot) <- gsub("\\_.+", "", colnames(ccle_prot))

colnames(ccle_prot)[1] <- "Uniprot_Acc"

# Transpose the data.table to match with other data type tables

t <- transpose(ccle_prot, make.names = "Uniprot_Acc")

# Check if transpose worked as intended

as.numeric(unlist(t[1,])) == as.numeric(unlist(ccle_prot[,2]))

as.numeric(unlist(t[2,])) == as.numeric(unlist(ccle_prot[,3]))

# Add cell lines

t$stripped_cell_line_name <- colnames(ccle_prot)[-1]

# Merge with sample info to have cell line name in addition to DepMap ID

t <- merge(t, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")], by = "stripped_cell_line_name")

# Move to front

setcolorder(t, neworder = c("DepMap_ID", "stripped_cell_line_name", "primary_disease"))

t[1:5, 1:10]

# Save

fwrite(t, paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_ProteinQuant.csv"), sep = ',')

### Get proteins that are observed in all cell lines

# Create the same transposed table as above

# Remove all rows and columns that have any NA in them

prot_nona <- na.omit(ccle_prot)

which(is.na(prot_nona))

# Transpose the data.table to match with other data type tables

t <- transpose(prot_nona, make.names = "Uniprot_Acc")

# Check if transpose worked as intended

as.numeric(unlist(t[1,])) == as.numeric(unlist(prot_nona[,2]))

as.numeric(unlist(t[2,])) == as.numeric(unlist(prot_nona[,3]))

# Add cell lines

t$stripped_cell_line_name <- colnames(prot_nona)[-1]

# Merge with sample info to have cell line name in addition to DepMap ID

t <- merge(t, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")], by = "stripped_cell_line_name")

# Move to front

setcolorder(t, neworder = c("DepMap_ID", "stripped_cell_line_name", "primary_disease"))

t[1:5, 1:10]

# Now we have ~5000 proteins that are available in all samples

dim(t)

# We have 3 duplicates

sum(duplicated(t$stripped_cell_line_name))

# Save

fwrite(t, paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv"), sep = ',')

# ccle_prot <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv"))

dim(ccle_prot)

ccle_prot[1:5, 1:5]

ccle_prot[, DepMap_ID := NULL]

# fwrite(ccle_prot, paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv"), sep = ',')

# DIMENSIONS OF PROTEIN QUANTITY DATA: 378 X 5155

ccle_prot <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv"))

anyDuplicated(ccle_prot$stripped_cell_line_name)

# ==== Mutation data cleanup ====

rm(list = ls(pattern = "ccle"))

require(data.table)

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

ccle_mut <- fread(paste0(path, "Data/DepMap/21Q2/CCLE_mutations.csv"))

table(ccle_mut$isCOSMIChotspot)

table(ccle_mut$isTCGAhotspot)

table(ccle_mut$Variant_Type)

length(unique(ccle_mut$DepMap_ID))

dim(ccle_mut)

ccle_mut[1,]

colnames(ccle_mut)

# Calculate number of mutations per cell line

temp <- ccle_mut[, c("Variant_Type", "DepMap_ID")]

temp[, nMut := .N, by = "DepMap_ID"]

temp

unique(temp$Variant_Type)

# For simplicity, extract only SNP data for now: this discards ~90,000 mutations

# ccle_mut <- ccle_mut[Variant_Type == "SNP"]

dim(ccle_mut)

t <- ccle_mut[, c("DepMap_ID", "Chromosome", "Strand", "Start_position", "End_position")]

dim(unique(t))

length(unique(ccle_mut$DepMap_ID))

# Keep relevant columns/features

# Aside: Should the sequence change be provided, or just whether the SNP is deleterious or not?

ccle_mut <- ccle_mut[, c("DepMap_ID", "Hugo_Symbol", "Chromosome", "Start_position", "End_position", "Strand",

             "Variant_Classification", "Variant_Type", "isDeleterious",

             "isTCGAhotspot", "isCOSMIChotspot", "Genome_Change", "cDNA_Change")]

dim(ccle_mut)

length(unique(ccle_mut$DepMap_ID))

table(ccle_mut$isDeleterious)

table(ccle_mut$isTCGAhotspot)

table(ccle_mut$isCOSMIChotspot)

# ==== CCLE Mut Overlap with COSMIC CGC ====

# Perhaps it's best to use the mutations in genes that COSMIC considers important, like another paper in

# the field (~500 genes)

# Or, we can use a binary vector for genes and whether they have a deleterious mutation: this will result in 

# ~20,000 parameters

length(unique(ccle_mut$Hugo_Symbol))

length(unique(ccle_mut[isCOSMIChotspot == T]$Hugo_Symbol))

length(unique(ccle_mut[isTCGAhotspot == T]$Hugo_Symbol))

length(unique(ccle_mut[isDeleterious == T]$Hugo_Symbol))

tcga_hotspot_genes <- unique(ccle_mut[isTCGAhotspot == T]$Hugo_Symbol)

# Read COSMIC Cancer Gene Census data

cgc <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/COSMIC/cancer_gene_census.csv")

dim(cgc)

cgc[1:5, 1:20]

length(unique(cgc$`Gene Symbol`))

length(unique(cgc$HGVSG))

# Get Genes in this census

cgc_genes <- unique(cgc$`Gene Symbol`)

cgc[Tier == 1]

length(unique(cgc$`Genome Location`))  # 922,732

# rm(cgc)

# Subset DepMap mutations based on the CGC genes

sum(unique(ccle_mut$Hugo_Symbol) %in% unique(cgc_genes))

ccle_mut <- ccle_mut[Hugo_Symbol %in% cgc_genes]

length(unique(ccle_mut$DepMap_ID))

sum(ccle_mut$isDeleterious)

ccle_mut[Variant_Classification == "Missense_Mutation"]

length(unique(ccle_mut[isDeleterious == T]$Hugo_Symbol))

ccle_mut[isDeleterious == T]

# TODO: Use CGC to check for overlap with CCLE cell lines, then collapse to whether each of the 700 genes for

# that cell line has a mutation listed in the CGC

length(unique(cgc$`Mutation genome position`))  # ~922,000 unique mutations

unique(ccle_mut$NCBI_Build)  # CCLE is with GRCh 37

unique(cgc$GRCh)  # CGC has GRCh 38

# We must "lift over" the mutations from 37 to 38 before checking for overlap

if (!require(liftOver)) {

    BiocManager::install("liftOver")

    require(liftOver)

    require(rtracklayer)

}

# liftOver requires a chain file to convert 37 to 38: http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/

chain_path <- paste0(path, "Data/hg19ToHg38.over.chain")

grch_37_38_chain <- import.chain(chain_path)

# Must add "chr" to start of chromosome names

ccle_mut$Chromosome <- paste0("chr", ccle_mut$Chromosome)

# Must convert positions to GRanges

ccle_mut_gr <- makeGRangesFromDataFrame(df = ccle_mut, keep.extra.columns = T,

                         seqnames.field = "Chromosome", start.field = "Start_position",

                         end.field = "End_position", strand.field = "Strand")

length(unique(ccle_mut_gr$DepMap_ID))

# Lift over

lifted_ccle_mut <- liftOver(x = ccle_mut_gr, chain = grch_37_38_chain)

# Convert GRangesList to GRanges

lifted_ccle_mut <- unlist(lifted_ccle_mut)

# Convert back to data.table

lifted_ccle_mut <- as.data.table(lifted_ccle_mut)

# Note: Genome_Change is now out of date!

# Remove chr from seqnames

lifted_ccle_mut$seqnames <- gsub("chr", "", lifted_ccle_mut$seqnames)

# Can find the overlap of Mutation genome position in CGC with a newly created column based on CCLE positions

lifted_ccle_mut[, Mutation_Position := paste0(seqnames, ':', start, '-', end)]

ccle_mut$seqnames <- gsub("chr", "", ccle_mut$Chromosome)

ccle_mut[, Mutation_Position := paste0(seqnames, ':', as.character(Start_position), '-', as.character(End_position))]

length(unique(lifted_ccle_mut$DepMap_ID))

sum(ccle_mut$Mutation_Position %in% unique(cgc$`Genome Location`))

# Now find the overlap with CGC (which already has GRCh38)

subset <- lifted_ccle_mut[Mutation_Position %in% unique(cgc$`Genome Location`)]

table(subset$Variant_Type)

length(unique(subset$DepMap_ID))

# IMPORTANT! There is a loss of 8 cell lines (which do not have a mutation that is in

# CGC) using the Tier 1 data only

# Alternative (March 2021) ====

# Take those mutations that are COSMIC or TCGA hotspots, ignoring CGC

subset <- ccle_mut[isTCGAhotspot | isCOSMIChotspot]

uniqueN(subset$Hugo_Symbol)

table(ccle_mut[isTCGAhotspot == T]$Variant_Classification)

uniqueN(ccle_mut[isTCGAhotspot == T]$Hugo_Symbol)

uniqueN(ccle_mut[isCOSMIChotspot == T]$Hugo_Symbol)

uniqueN(subset$isTCGAhotspot)

### Create a vector of mutations for each cell line with the CGC genes

length(unique(subset$Hugo_Symbol))

sub_dcast <- dcast.data.table(data = subset[, c("DepMap_ID", "Hugo_Symbol")],

                 formula = DepMap_ID ~ Hugo_Symbol, fun.aggregate = length, value.var = "DepMap_ID")

dim(sub_dcast)

sub_dcast[1:5, 1:50]

sum(sub_dcast$A1BG)

sum(sub_dcast$A1CF)

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

sub_dcast <- merge(sub_dcast, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")],

                  by = "DepMap_ID")

setcolorder(sub_dcast, c("DepMap_ID", "stripped_cell_line_name", "primary_disease"))

sub_dcast[1:5, 1:50]

# Save

fwrite(sub_dcast, paste0(path, "Data/DRP_Training_Data/DepMap_21Q2_Mutations_by_Cell.csv"), sep = ',')

dim(cgc_muts)

cgc_muts[1:5, 1:5]

typeof(cgc_muts[1,2])

temp <- fread("Data/DRP_Training_Data/DepMap_21Q2_Mutations_by_Cell.csv")

dim(temp)

temp[1:5, 1:50]

# # Attach the cell line name and primary disease

# # cgc_muts <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_CGC_Mutations_by_Cell.csv"))

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

# cgc_muts <- merge(cgc_muts, depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease")],

#                   by = "DepMap_ID")

# setcolorder(cgc_muts, neworder = c("stripped_cell_line_name", colnames(cgc_muts)[-ncol(cgc_muts)]))

# 

# # Save

# fwrite(cgc_muts, paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_CGC_Mutations_by_Cell.csv"), sep = ',')

# cgc_muts <- fread(paste0(path, "Data/DRP_Training_Data/DepMap_20Q2_CGC_Mutations_by_Cell.csv"))

# cgc_muts[1:5, 1:5]

# cgc_muts[, DepMap_ID := NULL]

# DIMENSIONS OF CGC MUTATIONAL DATA: 1733 X 697

# ==== miRNA Data Cleanup ====

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

require(data.table)

depmap_samples <- fread("Data/DRP_Training_Data/DepMap_21Q2_Line_Info.csv")

ccle_mirna <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_miRNA_20181103.gct")

dim(ccle_mirna)

anyNA(ccle_mirna)

ccle_mirna[1:5, 1:5]

min(ccle_mirna[, -c(1:2)], na.rm = T)

max(ccle_mirna[, -c(1:2)], na.rm = T)

ccle_mirna <- transpose(ccle_mirna, keep.names = "Name")

dim(ccle_mirna)

ccle_mirna[1:5, 1:5]

ccle_mirna$Name

sum(duplicated(unlist(ccle_mirna[2, ])))

ccle_mirna <- ccle_mirna[-1,]

ccle_mirna[1:5, 1:5]

colnames(ccle_mirna) <- unlist(ccle_mirna[1,])

ccle_mirna <- ccle_mirna[-1,]

# Clean cell line name

ccle_mirna$Description <- gsub(pattern = "\\_.+", replacement = "", ccle_mirna$Description)

colnames(ccle_mirna)[1] <- "stripped_cell_line_name"

ccle_mirna <- merge(ccle_mirna, depmap_samples[, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                    by = "stripped_cell_line_name")

dim(ccle_mirna)

ccle_mirna[1:5, 1:5]

setcolorder(ccle_mirna, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

# fwrite(ccle_mirna, paste0(path, "Data/DRP_Training_Data/DepMap_2019_miRNA.csv"), sep = ',')

rm(ccle_mirna)

# ccle_mirna <- fread("Data/DRP_Training_Data/DepMap_2019_miRNA.csv")

# dim(ccle_mirna)

# ==== Metabolomics Data Cleanup ====

depmap_samples <- fread("Data/DRP_Training_Data/DepMap_21Q2_Line_Info.csv")

ccle_metab <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_metabolomics_20190502.csv")

dim(ccle_metab)

ccle_metab[1:5, 1:5]

min(ccle_metab[, -c(1:2)], na.rm = T)

max(ccle_metab[, -c(1:2)], na.rm = T)

anyNA(ccle_metab)

sum(is.na(ccle_metab))

which(is.na(ccle_metab), arr.ind = T)

ccle_metab[which(is.na(ccle_metab), arr.ind = T)]

ccle_metab[554, 2]  # DepMap_ID is NA

min(ccle_metab[, -c(1:2)])

ccle_metab <- merge(ccle_metab[, -1], depmap_samples[, c("DepMap_ID", "stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                    by = "DepMap_ID")

setcolorder(ccle_metab, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

ccle_metab$DepMap_ID <- NULL

dim(ccle_metab)

fwrite(ccle_metab, paste0(path, "Data/DRP_Training_Data/DepMap_2019_Metabolomics.csv"), sep = ',')

rm(ccle_metab)

temp <- fread("Data/DRP_Training_Data/DepMap_2019_Metabolomics.csv")

dim(temp)

# ==== RPPA Data Cleanup ====

ccle_rppa <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_RPPA_20181003.csv")

dim(ccle_rppa)

ccle_rppa[1:5, 1:5]

anyNA(ccle_rppa)

min(ccle_rppa[, -c(1:2)], na.rm = T)  # has negative values

max(ccle_rppa[, -c(1:2)], na.rm = T)

ccle_rppa$V1 <- gsub(pattern = "\\_.+", replacement = "", ccle_rppa$V1)

colnames(ccle_rppa)[1] <- "stripped_cell_line_name"

ccle_rppa <- merge(ccle_rppa, depmap_samples[, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                    by = "stripped_cell_line_name")

dim(ccle_rppa)

ccle_rppa[1:5, 1:10]

setcolorder(ccle_rppa, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

fwrite(ccle_rppa, paste0(path, "Data/DRP_Training_Data/DepMap_2019_RPPA.csv"), sep = ',')

rm(ccle_rppa)

temp <- fread("Data/DRP_Training_Data/DepMap_2019_RPPA.csv")

dim(temp)

# ==== Chromatin Profiling Data Cleanup ====

ccle_chrom <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_GlobalChromatinProfiling_20181130.csv")

dim(ccle_chrom)

ccle_chrom[1:5, 1:10]

min(ccle_chrom[, -c(1:2)], na.rm = T)  # has negative values

max(ccle_chrom[, -c(1:2)], na.rm = T)

anyNA(ccle_chrom)

sum(is.na(ccle_chrom))  # 842 NA values

unique(which(is.na(ccle_chrom), arr.ind = T)[,2])

length(unique(which(is.na(ccle_chrom), arr.ind = T)[,2]))  # 26 columns have NAs

# Convert NA to 0

setnafill(ccle_chrom, fill = 0, cols = unique(which(is.na(ccle_chrom), arr.ind = T)[,2]))

anyNA(ccle_chrom)

ccle_chrom$CellLineName <- gsub(pattern = "\\_.+", replacement = "", ccle_chrom$CellLineName)

colnames(ccle_chrom)[1] <- "stripped_cell_line_name"

dim(ccle_chrom)

ccle_chrom <- merge(ccle_chrom, depmap_samples[, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                   by = "stripped_cell_line_name")

ccle_chrom$BroadID <- NULL

dim(ccle_chrom)

ccle_chrom[1:5, 1:10]

setcolorder(ccle_chrom, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

# fwrite(ccle_chrom, paste0(path, "Data/DRP_Training_Data/DepMap_2019_ChromatinProfiling.csv"), sep = ',')

rm(ccle_chrom)

temp <- fread(paste0("Data/DRP_Training_Data/DepMap_2019_ChromatinProfiling.csv"))

dim(temp)

# ==== Fusion Data Cleanup ====

ccle_fusion <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_fusions.csv")

dim(ccle_fusion)

ccle_fusion[1:5, 1:17]

length(unique(ccle_fusion$FusionName))

length(unique(ccle_fusion$DepMap_ID))

unique(ccle_fusion$SpliceType)

quantile(ccle_fusion$FFPM)

ccle_fusion$CellLineName <- gsub(pattern = "\\_.+", replacement = "", ccle_fusion$CellLineName)

colnames(ccle_fusion)[1] <- "stripped_cell_line_name"

dim(ccle_fusion)

ccle_fusion <- merge(ccle_fusion, depmap_samples[, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                    by = "stripped_cell_line_name")

ccle_fusion$BroadID <- NULL

dim(ccle_fusion)

ccle_fusion[1:5, 1:10]

setcolorder(ccle_fusion, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

fwrite(ccle_fusion, paste0(path, "Data/DRP_Training_Data/DepMap_2019_GeneFusion.csv"), sep = ',')

rm(ccle_fusion)

# ==== Exon Usage Ratio Data Cleanup ====

require(data.table)

setDTthreads(8)

ccle_exon <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_RNAseq_ExonUsageRatio_20180929.gct")

dim(ccle_exon)

ccle_exon[1:10, 1:17]

transpose(ccle_exon, keep.names = "exon")

length(unique(ccle_exon$FusionName))

length(unique(ccle_exon$DepMap_ID))

unique(ccle_exon$SpliceType)

quantile(ccle_exon$FFPM)

ccle_exon$CellLineName <- gsub(pattern = "\\_.+", replacement = "", ccle_exon$CellLineName)

colnames(ccle_exon)[1] <- "stripped_cell_line_name"

dim(ccle_exon)

ccle_exon <- merge(ccle_exon, depmap_samples[, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype")],

                     by = "stripped_cell_line_name")

ccle_exon$BroadID <- NULL

dim(ccle_exon)

ccle_exon[1:5, 1:10]

setcolorder(ccle_exon, c("stripped_cell_line_name", "primary_disease", "lineage", "lineage_subtype"))

# fwrite(ccle_exon, paste0(path, "Data/DRP_Training_Data/DepMap_2019_ExonUsageRatio.csv"), sep = ',')

rm(ccle_exon)

# ==== RRBS Profiling Data Cleanup ====

require(data.table)

setDTthreads(8)

# === TSS

ccle_tss <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE_RRBS_TSS1kb_20181022.txt")

dim(ccle_tss)

ccle_tss[1:5, 1:5]

length(unique(ccle_tss$cluster_id))

ccle_tss <- transpose(ccle_tss[, -2], keep.names = "cluster_id")

colnames(ccle_tss) <- unlist(ccle_tss[1,])

ccle_tss <- ccle_tss[-1,]

# === Promoter

ccle_tss <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DepMap/Extra/CCLE")

dim(ccle_tss)

ccle_tss[1:5, 1:5]

length(unique(ccle_tss$cluster_id))

ccle_tss <- transpose(ccle_tss[, -2], keep.names = "cluster_id")

colnames(ccle_tss) <- unlist(ccle_tss[1,])

ccle_tss <- ccle_tss[-1,]

# === Enhancers

# ==== GDSC Cell Line Characterization Data Cleanup ====

require(data.table)

gdsc_line_info <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/GDSC/GDSC_Line_Info.csv")

gdsc_line_info <- gdsc_line_info[-c(1,2,3, 1005), c("V1", "V2", "V9")]

colnames(gdsc_line_info) <- c("stripped_cell_line_name", "COSMIC_ID", "primary_disease")

gdsc_line_info$stripped_cell_line_name <- gsub("[^[:alnum:] ]", "", gdsc_line_info$stripped_cell_line_name)

gdsc_line_info$stripped_cell_line_name <- gsub(" ", "", gdsc_line_info$stripped_cell_line_name)

gdsc_line_info$stripped_cell_line_name <- toupper(gdsc_line_info$stripped_cell_line_name)

# Match with DepMap cell lines

depmap_line_info <- fread("Data/DRP_Training_Data/DepMap_21Q2_Line_Info.csv")

sum(gdsc_line_info$stripped_cell_line_name %in% depmap_line_info$stripped_cell_line_name) / nrow(gdsc_line_info)  # 98%

gdsc_line_info[!(gdsc_line_info$stripped_cell_line_name %in% depmap_line_info$stripped_cell_line_name)]

gdsc_line_info[stripped_cell_line_name == "EOL1CELL"] <- "EOL1"

gdsc_line_info[stripped_cell_line_name == "SR"] <- "SR786"

gdsc_line_info[stripped_cell_line_name == "U266"] <- "U266B1"

gdsc_line_info[stripped_cell_line_name == "G292"] <- "G292CLONEA141B1"

gdsc_line_info[stripped_cell_line_name == "NCISNU1"] <- "SNU1"

gdsc_line_info[stripped_cell_line_name == "NCISNU5"] <- "SNU5"

gdsc_line_info[stripped_cell_line_name == "NCISNU16"] <- "SNU16"

gdsc_line_info[stripped_cell_line_name == "7860"] <- "786O"

gdsc_line_info[stripped_cell_line_name == "U031"] <- "UO31"

gdsc_line_info[stripped_cell_line_name == "H3255"] <- "NCIH3255"

gdsc_line_info[stripped_cell_line_name == "NCIH510A"] <- "NCIH510"

gdsc_line_info[stripped_cell_line_name == "U251"] <- "U251MG"

gdsc_line_info[stripped_cell_line_name == "WM793B"] <- "WM793"

gdsc_line_info[stripped_cell_line_name == "OVCAR3"] <- "NIHOVCAR3"

gdsc_line_info[stripped_cell_line_name == "SCC90"] <- "UPCISCC090"

gdsc_line_info[COSMIC_ID == "1299064"] <- "TDOTT"

gdsc_line_info[COSMIC_ID == "930299"] <- "TT"

gdsc_line_info[COSMIC_ID == "909976"] <- "KMH2"

gdsc_line_info[COSMIC_ID == "1298167"] <- "KMHDASH2"

sum(gdsc_line_info$stripped_cell_line_name %in% depmap_line_info$stripped_cell_line_name) / nrow(gdsc_line_info)  # 100%

# ==== GDSC Expression Data Cleanup ====

gdsc_exp <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/GDSC/GDSC_Expression.txt")

gdsc_exp[1:5, 1:5]

# Delete invalid gene symbol rows

gdsc_exp <- gdsc_exp[GENE_SYMBOLS != ""]

# Rename Expression columns

colnames(gdsc_exp) <- gsub("DATA.", "", colnames(gdsc_exp))

gdsc_exp <- gdsc_exp[, !"GENE_title", with = F]

# Must transpose to have cell lines on each row

gdsc_exp <- transpose(gdsc_exp, keep.names = "GENE_SYMBOLS")

gdsc_exp[1:5, 1:5]

colnames(gdsc_exp) <- unlist(gdsc_exp[1,])

gdsc_exp <- gdsc_exp[-1,]

colnames(gdsc_exp)[1] <- "COSMIC_ID"

# gdsc_exp[COSMIC_ID == "1299064"][1:5, 1:5]

# gdsc_exp[COSMIC_ID == "930299"][1:5, 1:5]

# Merge with cell line names

gdsc_exp <- merge(gdsc_exp, gdsc_line_info, by = "COSMIC_ID")

gdsc_exp$COSMIC_ID <- NULL

setcolorder(gdsc_exp, order(colnames(gdsc_exp)))

setcolorder(gdsc_exp, c("stripped_cell_line_name", "primary_disease"))

gdsc_exp[1:5, 1:5]

# Match genes with those from DepMap Cell lines

depmap_exp <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DRP_Training_Data/DepMap_21Q2_Expression.csv")

# depmap_exp <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/Data/DRP_Training_Data/DepMap_21Q2_Training_Expression.csv")

depmap_exp[1:5, 1:5]

dim(gdsc_exp)

dim(depmap_exp)

t1 <- colnames(gdsc_exp)

t2 <- colnames(depmap_exp)

sum(t1 %in% t2) / ncol(gdsc_exp)  # 91%

colnames(gdsc_exp)[!colnames(gdsc_exp) %in% colnames(depmap_exp)]

# Update gene names with biomaRt

# require(biomaRt)

# ensembl = useMart(

#   biomart = "ENSEMBL_MART_ENSEMBL",

#   # host = "feb2014.archive.ensembl.org",

#   path = "/biomart/martservice",

#   dataset = "hsapiens_gene_ensembl"

# )

# atts <- listAttributes(ensembl, what = "name")

# atts[grep(pattern = "gene", x = atts, ignore.case = T)]

# depmap_dict <- biomaRt::getBM(attributes = c("ensembl_gene_id",

#                                           "external_gene_name", "wikigene_name"),

#                            filters = "external_gene_name", values = colnames(depmap_exp),

#                            mart = ensembl)

# depmap_dict <- as.data.table(depmap_dict)

# depmap_dict <- unique(depmap_dict[, 2:3])

# 

# gdsc_dict <- biomaRt::getBM(attributes = c("ensembl_gene_id",

#                                              "external_gene_name", "wikigene_name"),

#                               filters = "external_gene_name", values = colnames(gdsc_exp),

#                               mart = ensembl)

# gdsc_dict <- as.data.table(gdsc_dict)

# Get updated gene names from: https://www.genenames.org/tools/multi-symbol-checker/

all_gene_names <- data.table(unique(c(colnames(depmap_exp), colnames(gdsc_exp)))[-(1:2)])

fwrite(all_gene_names, "GDSC_DepMap_All_Gene_Names.txt")

hgnc_check <- fread("/Users/ftaj/OneDrive - University of Toronto/Drug_Response/hgnc-symbol-check.csv", skip = 1)

depmap_hgnc_check <- hgnc_check[Input %in% colnames(depmap_exp)[-(1:2)]]

names(depmap_exp)[match(depmap_hgnc_check$Input, names(depmap_exp))] <- depmap_hgnc_check$`Approved symbol`

gdsc_hgnc_check <- hgnc_check[Input %in% colnames(gdsc_exp)[-(1:2)]]

names(gdsc_exp)[match(gdsc_hgnc_check$Input, names(gdsc_exp))] <- gdsc_hgnc_check$`Approved symbol`

t1 <- colnames(gdsc_exp)

t2 <- colnames(depmap_exp)

sum(t1 %in% t2) / ncol(gdsc_exp)  # 97.3%

sum(t2 %in% t1) / ncol(depmap_exp)  # 88.6%

sum(t2 %in% t1)  # 16988 

sum(t1 %in% t2)  # 16948 

colnames(gdsc_exp)[!colnames(gdsc_exp) %in% colnames(depmap_exp)]

fwrite(gdsc_exp, "Data/DRP_Training_Data/GDSC_Expression.csv")

# Subset both datasets's genes based on data available to both

shared <- intersect(colnames(gdsc_exp), colnames(depmap_exp))

length(shared)

gdsc_exp_sub <- gdsc_exp[, shared, with = F]

gdsc_exp_sub[1:5, 1:5]

depmap_exp_sub <- depmap_exp[, shared, with = F]

depmap_exp_sub[1:5, 1:5]

# Save

fwrite(gdsc_exp_sub, "Data/DRP_Training_Data/GDSC_Training_Expression.csv")

fwrite(depmap_exp_sub, "Data/DRP_Training_Data/DepMap_21Q2_Training_Expression.csv")

# ==== Drug Sensitivity Data Cleanup ====

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

require(data.table)

require(webchem)

# BiocManager::install("ChemmineR")

require(ChemmineR)

options(chemspider_key = "N98K4aOip0VpcSc8F9GilqIIktLt0hux")

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

ctrp <- fread("Data/DRP_Training_Data/CTRP_AUC_SMILES.txt")

gdsc1 <- fread("Data/DRP_Training_Data/GDSC1_AUC_SMILES.txt")

gdsc2 <- fread("Data/DRP_Training_Data/GDSC2_AUC_SMILES.txt")

# Clean up duplicate with missing pubchem

cpd_info_1 <- fread(paste0(path, "Data/GDSC/GDSC1_Drug_Info.csv"))

cpd_info_1[drug_name == unique(cpd_info_1[, c("drug_name", "pubchem")])[anyDuplicated(unique(cpd_info_1[, c("drug_name", "pubchem")])$drug_name),]$drug_name]

cpd_info_1 <- cpd_info_1[drug_id != 476]

cpd_info_1 <- cpd_info_1[drug_id != 1490]

cpd_info_1 <- cpd_info_1[drug_id != 1496]

cpd_info_1 <- cpd_info_1[drug_id != 1386]

cpd_info_1 <- cpd_info_1[drug_id != 1402]

cpd_info_1 <- cpd_info_1[drug_id != 1393]

nrow(cpd_info_1[pubchem == "-"])

sum(cpd_info_1$drug_name %in% unique(ctrp$cpd_name))

# Subset for valid pubchem IDs

cpd_info_1 <- cpd_info_1[pubchem != "-"]

cpd_info_1 <- cpd_info_1[pubchem != "none"]

cpd_info_1 <- cpd_info_1[pubchem != "several"]

cpd_info_1$pubchem <- as.numeric(cpd_info_1$pubchem)

cpd_1_smiles <- webchem::pc_prop(cid = cpd_info_1$pubchem, properties = "CanonicalSMILES")

cpd_info_1 <- merge(cpd_info_1, cpd_1_smiles, by.x = "pubchem", by.y = "CID")

# Save

fwrite(cpd_info_1, "Data/GDSC/GDSC1_VALID_Drug_Info.csv")

cpd_info_2 <- fread(paste0(path, "Data/GDSC/GDSC2_Drug_Info.csv"))

cpd_info_2[drug_name == unique(cpd_info_2[, c("drug_name", "pubchem")])[anyDuplicated(unique(cpd_info_2[, c("drug_name", "pubchem")])$drug_name),]$drug_name]

cpd_info_2 <- cpd_info_2[drug_id != 1811]

cpd_info_2 <- cpd_info_2[drug_id != 1806]

cpd_info_2 <- cpd_info_2[drug_id != 1819]

cpd_info_2 <- cpd_info_2[drug_id != 1816]

cpd_info_2[pubchem == "25227436, 42602260"]$pubchem <- "25227436"

cpd_info_2[pubchem == "11719003, 58641927"]$pubchem <- "11719003"

cpd_info_2[pubchem == "66577015, 16654980"]$pubchem <- "66577015"

cpd_info_2[pubchem == "11719003, 58641927"]$pubchem <- "11719003"

nrow(cpd_info_2[pubchem == "-"])

sum(cpd_info_2$pubchem %in% cpd_info_1$pubchem) / nrow(cpd_info_2)

cpd_info_2 <- cpd_info_2[pubchem != "-"]

cpd_info_2 <- cpd_info_2[pubchem != "none"]

cpd_info_2 <- cpd_info_2[pubchem != "several"]

cpd_info_2$pubchem <- as.numeric(cpd_info_2$pubchem)

cpd_2_smiles <- webchem::pc_prop(cid = cpd_info_2$pubchem, properties = "CanonicalSMILES")

cpd_info_2 <- merge(cpd_info_2, cpd_2_smiles, by.x = "pubchem", by.y = "CID")

# Save

fwrite(cpd_info_2, "Data/GDSC/GDSC2_VALID_Drug_Info.csv")

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

# ==== GDSC ====

require(stringr)

gdsc1 <- fread(paste0(path, "Data/GDSC/GDSC1_Fitted_Dose_Response.csv"))

sum(unique(gdsc1$CELL_LINE_NAME) %in% depmap_samples$stripped_cell_line_name) / length(unique(gdsc1$CELL_LINE_NAME))  # 0.22

sum(toupper(unique(gdsc1$CELL_LINE_NAME)) %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(gdsc1$CELL_LINE_NAME))  # 0.24

sum(str_remove_all(toupper(unique(gdsc1$CELL_LINE_NAME)), "-") %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(gdsc1$CELL_LINE_NAME))  # 0.9696049

dim(gdsc1)  # 310K Combinations

colnames(gdsc1)

sum(gdsc1$AUC == 0)

min(gdsc1$AUC)

max(gdsc1$AUC)

# Count unique combinations in GDSC1

length(unique(unique(gdsc1[, c("DRUG_NAME", "CELL_LINE_NAME")])$CELL_LINE_NAME))  # 987

length(unique(unique(gdsc1[, c("DRUG_NAME", "CELL_LINE_NAME")])$DRUG_NAME))  # 345

nrow(unique(unique(gdsc1[, c("DRUG_NAME", "CELL_LINE_NAME")]))) # 292,849

gdsc1_final <- merge(unique(gdsc1[, c("DRUG_NAME", "CELL_LINE_NAME", "AUC")]), unique(cpd_info_1[, c("drug_name", "CanonicalSMILES")]), by.x = "DRUG_NAME", by.y = "drug_name")

colnames(gdsc1_final) <- c("cpd_name", "ccl_name", "area_under_curve", "cpd_smiles")

# Save

fwrite(gdsc1_final, "Data/DRP_Training_Data/GDSC1_AUC_SMILES.txt")

unique(gdsc1_pubchem$DRUG_NAME)

# gdsc1_cs_ids <- webchem::get_csid(query = unique(gdsc1$DRUG_NAME), from = "name", match = "all", verbose = T)

gdsc1_cs_ids <- webchem::cir_query(identifier = unique(gdsc1$DRUG_NAME), representation = "smiles", verbose = T, )

# Count unique combinations in GDSC2

gdsc2 <- fread(paste0(path, "Data/GDSC/GDSC2_Fitted_Dose_Response.csv"))

sum(unique(gdsc2$CELL_LINE_NAME) %in% depmap_samples$stripped_cell_line_name) / length(unique(gdsc2$CELL_LINE_NAME))  # 0.2311496

sum(toupper(unique(gdsc2$CELL_LINE_NAME)) %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(gdsc2$CELL_LINE_NAME))  # 0.2546354

sum(str_remove_all(toupper(unique(gdsc2$CELL_LINE_NAME)), "-") %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(gdsc2$CELL_LINE_NAME))  # 0.9678616

gdsc2_cpd_smiles <- webchem::cir_query(identifier = unique(gdsc2$DRUG_NAME), representation = "smiles", verbose = T)

dim(gdsc2)  # 135K Combinations

colnames(gdsc2)

length(unique(unique(gdsc2[, c("DRUG_NAME", "CELL_LINE_NAME")])$CELL_LINE_NAME))  # 809

length(unique(unique(gdsc2[, c("DRUG_NAME", "CELL_LINE_NAME")])$DRUG_NAME))  # 192

nrow(unique(unique(gdsc2[, c("DRUG_NAME", "CELL_LINE_NAME")]))) # 131,108

gdsc2_final <- merge(unique(gdsc2[, c("DRUG_NAME", "CELL_LINE_NAME", "AUC")]), unique(cpd_info_2[, c("drug_name", "CanonicalSMILES")]), by.x = "DRUG_NAME", by.y = "drug_name")

colnames(gdsc2_final) <- c("cpd_name", "ccl_name", "area_under_curve", "cpd_smiles")

# Save

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

# Count overlap of drugs and cell lines

sum(unique(gdsc1$DRUG_NAME) %in% unique(gdsc2$DRUG_NAME)) # Drug Overlap: 88

sum(unique(gdsc1$CELL_LINE_NAME) %in% unique(gdsc2$CELL_LINE_NAME))  # Cell Line Overlap: 808

# ==== CTRP ====

require(data.table)

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

# NOTE: Newer and better AUC calculation in PharmacoGx.R file!

ctrp_curves <- fread(paste0(path, "Data/CTRP/v20.data.curves_post_qc.txt"))

exper_data <- fread(paste0(path, "Data/CTRP/v20.meta.per_experiment.txt"))

cell_data <- fread(paste0(path, "Data/CTRP/v20.meta.per_cell_line.txt"))

table(cell_data$ccl_availability)

# Merge sensitivity, experimental and cell line data

temp <- merge(unique(ctrp_curves[, c("experiment_id", "master_cpd_id")]),

              unique(exper_data[, c("experiment_id", "master_ccl_id")]),

              by = "experiment_id")

ctrp <- merge(temp, cell_data[, c("master_ccl_id", "ccl_name")], by = "master_ccl_id")

sum(unique(ctrp$ccl_name) %in% depmap_samples$stripped_cell_line_name) / length(unique(ctrp$ccl_name))  # 0.9492672

sum(toupper(unique(ctrp$ccl_name)) %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(ctrp$ccl_name))  # 0.9492672

sum(str_remove_all(toupper(unique(ctrp$ccl_name)), "-") %in% toupper(depmap_samples$stripped_cell_line_name)) / length(unique(ctrp$ccl_name))  # 0.9503946

# Add compound information

cpd_data <- fread(paste0(path, "Data/CTRP/v20.meta.per_compound.txt"))

ctrp <- merge(ctrp, cpd_data[, c("master_cpd_id", "cpd_name", "cpd_smiles")], by = "master_cpd_id")

# Add AUC curve information

ctrp_auc <- fread(paste0(path, "Data/CTRP/v20.data.curves_post_qc.txt"))

ctrp <- merge(ctrp, ctrp_auc[, c("experiment_id", "master_cpd_id", "area_under_curve")], by = c("experiment_id", "master_cpd_id"))

# Save

fwrite(ctrp, paste0(path, "Data/DRP_Training_Data/CTRP_AUC_SMILES.txt"))

# Add primary disease information. NOTE: This removes some DR data as 45 cell lines in CTRPv2 cannot be paired with DepMap!!!

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

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

sum(unique(ctrp$ccl_name) %in% unique(line_info$stripped_cell_line_name))  # 150

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

ctrp$other_ccl_name <- str_replace(toupper(ctrp$ccl_name), "-", "")

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

ctrp$other_ccl_name <- NULL

setcolorder(ctrp, neworder = c("cpd_name", "ccl_name", "primary_disease", "area_under_curve", "cpd_smiles"))

fwrite(ctrp, "Data/DRP_Training_Data/CTRP_AUC_SMILES.txt")

# Experiment ID 

unique(ctrp[, c("master_ccl_id", "experiment_id")])

length(unique(ctrp$master_ccl_id))

length(unique(ctrp$experiment_id))

length(unique(ctrp$ccl_name))

length(unique(ctrp$master_cpd_id))

# Check overlap with GDSC 1 and 2

sum(unique(ctrp$ccl_name) %in% gdsc1$CELL_LINE_NAME)

sum(unique(ctrp$ccl_name) %in% gdsc2$CELL_LINE_NAME)

dim(ctrp)  # 395K Combinations

# ==== Chemical Data Cleanup ====

require(data.table)

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

chembl <- fread(paste0(path, "Data/chembl_27_chemreps.txt"))

dim(chembl)

# ==== EDA ======

require(data.table)

require(stringr)

require(ggplot2)

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

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

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

exp <- fread("Data/DRP_Training_Data/DepMap_21Q2_Expression.csv")

mut <- fread("Data/DRP_Training_Data/DepMap_21Q2_Mutations_by_Cell.csv")

cnv <- fread("Data/DRP_Training_Data/DepMap_21Q2_CopyNumber.csv")

prot <- fread("Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv")

mirna <- fread("Data/DRP_Training_Data/DepMap_2019_miRNA.csv")

metab <- fread("Data/DRP_Training_Data/DepMap_2019_Metabolomics.csv")

hist <- fread("Data/DRP_Training_Data/DepMap_2019_ChromatinProfiling.csv")

rppa <- fread("Data/DRP_Training_Data/DepMap_2019_RPPA.csv")

# pdb_table <- fread("Data/cell_annotation_table_1.1.1.csv")

# pdb_sub <- pdb_table[, c("CTRPv2.cellid", "CCLE.cellid")]

# pdb_sub <- pdb_sub[!is.na(CTRPv2.cellid) & !is.na(CCLE.cellid)]

exp[1:5., 1:5]

length(unique(ctrp$ccl_name))

length(unique(ctrp$cpd_name))

sum(unique(ctrp$ccl_name) %in% line_info$stripped_cell_line_name) / length(unique(ctrp$ccl_name))

ccl_names = toupper(ctrp$ccl_name)

ccl_names = unique(str_replace(ccl_names, "-", ""))

length(ccl_names)

sum(ccl_names %in% line_info$stripped_cell_line_name) / length(ccl_names)

line_info[!(stripped_cell_line_name %in% ccl_names)]

ctrp[ccl_name %like% "NIHOVCAR3"]

ctrp[ccl_name %like% "HEL"]

# sum(exp$stripped_cell_line_name %in% pdb_sub$CCLE.cellid) 

# Remove hyphens and convert all to upper case

# pdb_ccl_names = pdb_sub$CCLE.cellid

# pdb_ccl_names = str_replace(toupper(pdb_ccl_names), "-", "")

ctrp$ccl_name = str_replace(toupper(ctrp$ccl_name), "-", "")

exp_ccl_names = exp$stripped_cell_line_name

exp_ccl_names = str_replace(toupper(exp_ccl_names), "-", "")

mut_ccl_names = mut$stripped_cell_line_name

mut_ccl_names = str_replace(toupper(mut_ccl_names), "-", "")

cnv_ccl_names = cnv$stripped_cell_line_name

cnv_ccl_names = str_replace(toupper(cnv_ccl_names), "-", "")

sum(exp_ccl_names %in% ccl_names) / length(unique(ccl_names))

# sum(exp_ccl_names %in% pdb_ccl_names) / length(unique(pdb_ccl_names))

sum(mut_ccl_names %in% ccl_names) / length(unique(ccl_names)) * length(unique(ccl_names))

# sum(mut_ccl_names %in% pdb_ccl_names) / length(unique(pdb_ccl_names)) * length(unique(pdb_ccl_names))

ctrp[ccl_name %in% mut_ccl_names[mut_ccl_names %in% ccl_names]]   ### 308K!!!!! Not 144K

# ctrp[ccl_name %in% mut_ccl_names[mut_ccl_names %in% pdb_ccl_names]]   ### 302K!!!!! Not 144K

sum(cnv_ccl_names %in% ccl_names) / length(unique(ccl_names))

sum(exp_ccl_names %in% cnv_ccl_names) / length(unique(exp_ccl_names))

sum(cnv_ccl_names %in% exp_ccl_names) / length(unique(cnv_ccl_names))

dir.create(path = "Plots")

dir.create(path = "Plots/DepMap")

ggplot(data = line_info) +

  geom_bar(mapping = aes(x = primary_disease), stat = "count") +

  xlab("Primary Disease") +

  ylab("# of cell lines") + 

  ggtitle(label = "Proportion of Cancer Types in DepMap Data (overall)", subtitle = "20Q2 Version") +

  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(filename = "Plots/DepMap/DepMap_Cell_Lines_Proportion.pdf", device = "pdf")

prot[, 1:5]

unique(prot$stripped_cell_line_name)

mut$stripped_cell_line_name = str_replace(toupper(mut$stripped_cell_line_name), "-", "")

cnv$stripped_cell_line_name = str_replace(toupper(cnv$stripped_cell_line_name), "-", "")

exp$stripped_cell_line_name = str_replace(toupper(exp$stripped_cell_line_name), "-", "")

prot$stripped_cell_line_name = str_replace(toupper(prot$stripped_cell_line_name), "-", "")

mirna$stripped_cell_line_name = str_replace(toupper(mirna$stripped_cell_line_name), "-", "")

metab$stripped_cell_line_name = str_replace(toupper(metab$stripped_cell_line_name), "-", "")

hist$stripped_cell_line_name = str_replace(toupper(hist$stripped_cell_line_name), "-", "")

rppa$stripped_cell_line_name = str_replace(toupper(rppa$stripped_cell_line_name), "-", "")

ctrp$ccl_name = str_replace(toupper(ctrp$ccl_name), "-", "")

mut_line_info <- line_info[stripped_cell_line_name %in% unique(mut$stripped_cell_line_name)]  

cnv_line_info <- line_info[stripped_cell_line_name %in% unique(cnv$stripped_cell_line_name)]  

exp_line_info <- line_info[stripped_cell_line_name %in% unique(exp$stripped_cell_line_name)]  

prot_line_info <- line_info[stripped_cell_line_name %in% unique(prot$stripped_cell_line_name)]

mirna_line_info <- line_info[stripped_cell_line_name %in% unique(mirna$stripped_cell_line_name)]

metab_line_info <- line_info[stripped_cell_line_name %in% unique(metab$stripped_cell_line_name)]

hist_line_info <- line_info[stripped_cell_line_name %in% unique(hist$stripped_cell_line_name)]

rppa_line_info <- line_info[stripped_cell_line_name %in% unique(rppa$stripped_cell_line_name)]

ctrp_line_info <- line_info[stripped_cell_line_name %in% unique(ctrp$ccl_name)]

mut_line_info <- mut_line_info[, c("stripped_cell_line_name", "primary_disease")]

mut_line_info$data_type <- "Mutational"

cnv_line_info <- cnv_line_info[, c("stripped_cell_line_name", "primary_disease")]

cnv_line_info$data_type <- "Copy Number"

exp_line_info <- exp_line_info[, c("stripped_cell_line_name", "primary_disease")]

exp_line_info$data_type <- "Gene Expression"

prot_line_info <- prot_line_info[, c("stripped_cell_line_name", "primary_disease")]

prot_line_info$data_type <- "Protein Quantification"

mirna_line_info <- mirna_line_info[, c("stripped_cell_line_name", "primary_disease")]

mirna_line_info$data_type <- "miRNA Expression"

metab_line_info <- metab_line_info[, c("stripped_cell_line_name", "primary_disease")]

metab_line_info$data_type <- "Metabolomics"

hist_line_info <- hist_line_info[, c("stripped_cell_line_name", "primary_disease")]

hist_line_info$data_type <- "Histone Profiling"

rppa_line_info <- rppa_line_info[, c("stripped_cell_line_name", "primary_disease")]

rppa_line_info$data_type <- "RPPA"

ctrp_line_info <- ctrp_line_info[, c("stripped_cell_line_name", "primary_disease")]

ctrp_line_info$data_type <- "Dose-Response"

datatype_line_info <- rbindlist(list(mut_line_info, cnv_line_info, exp_line_info, prot_line_info,

                                     mirna_line_info, metab_line_info, hist_line_info, rppa_line_info,

                                     ctrp_line_info))

datatype_line_info[, count_per_disease := .N, by = "primary_disease"]

# .(count = .N, var = sum(VAR)), by = MNTH

ggplot(data = datatype_line_info) +

  geom_bar(mapping = aes(x = reorder(primary_disease, -count_per_disease),

                         fill = data_type), stat = "count", position = "dodge") +

  xlab("Primary Disease") +

  ylab("Number of cell lines") +

  labs(fill = "Data Type") +

  # ggtitle(label = "Proportion of Cancer Types in DepMap Data", subtitle = "By data type, 21Q2 Version - Overlap with CTRPv2: 79%") +

  theme(axis.text.x = element_text(angle = 45, hjust = 1),

        text = element_text(face = "bold"),

        legend.position = c(.9,.55))

ggsave(filename = "Plots/DepMap/DepMap_CTRP_Cell_Lines_Proportion.pdf", device = "pdf")

BiocManager::install("VennDiagram")

require(VennDiagram)

library(RColorBrewer)

myCol <- brewer.pal(5, "Pastel2")

# NOTE: The CTRPv2 here is from before ctrp was merged with cell line info to add primary disease!

venn.diagram(x = list(mut_line_info$stripped_cell_line_name,

                      cnv_line_info$stripped_cell_line_name,

                      exp_line_info$stripped_cell_line_name,

                      prot_line_info$stripped_cell_line_name,

                      unique(ctrp$ccl_name)),

             category.names = c("Mutational", "Copy Number", "Gene Expression", "Protein Quantification", "CTRPv2 Dose-Response"),

             filename = "Plots/DepMap/DepMap_CTRP_Cell_Lines_Venn.png",

             imagetype = "png",

             output = TRUE,

             height = 3000 ,

             width = 3000 ,

             resolution = 600,

             # Circles

             lwd = 2,

             # lty = 'blank',

             fill = myCol,

             # Numbers

             cex = .6,

             fontface = "bold",

             fontfamily = "sans",

             # Set names

             cat.cex = 0.6,

             cat.fontface = "bold",

             cat.default.pos = "outer",

             cat.pos = c(0, 0, -130, 150, 0),

             cat.dist = c(0.2, 0.2, 0.2, 0.2, 0.2),

             cat.fontfamily = "sans",

             # rotation = 1

)

# ==== Bottleneck Cancer Cell Lines ====

require(data.table)

require(stringr)

require(ggplot2)

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

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

mut <- fread("Data/DRP_Training_Data/DepMap_21Q2_Mutations_by_Cell.csv")

cnv <- fread("Data/DRP_Training_Data/DepMap_21Q2_CopyNumber.csv")

exp <- fread("Data/DRP_Training_Data/DepMap_21Q2_Expression.csv")

prot <- fread("Data/DRP_Training_Data/DepMap_20Q2_No_NA_ProteinQuant.csv")

mirna <- fread("Data/DRP_Training_Data/DepMap_2019_miRNA.csv")

metab <- fread("Data/DRP_Training_Data/DepMap_2019_Metabolomics.csv")

hist <- fread("Data/DRP_Training_Data/DepMap_2019_ChromatinProfiling.csv")

rppa <- fread("Data/DRP_Training_Data/DepMap_2019_RPPA.csv")

bottleneck_keys <- Reduce(intersect, list(ctrp$ccl_name, mut$stripped_cell_line_name, cnv$stripped_cell_line_name,

                                          exp$stripped_cell_line_name, prot$stripped_cell_line_name, mirna$stripped_cell_line_name,

                                          metab$stripped_cell_line_name, hist$stripped_cell_line_name, rppa$stripped_cell_line_name)

)

bottleneck_keys <- data.table(keys = bottleneck_keys)

fwrite(bottleneck_keys, "Data/DRP_Training_Data/bottleneck_keys.csv")

bottleneck_keys <- Reduce(intersect, list(ctrp$ccl_name, mut$stripped_cell_line_name, cnv$stripped_cell_line_name,

                                          exp$stripped_cell_line_name, prot$stripped_cell_line_name, mirna$stripped_cell_line_name,

                                          metab$stripped_cell_line_name, hist$stripped_cell_line_name, rppa$stripped_cell_line_name)

)

bottleneck_datatype_line_info <- datatype_line_info[stripped_cell_line_name %in% bottleneck_keys]

bottleneck_datatype_line_info[, count_per_disease := .N, by = "primary_disease"]

ggplot(data = bottleneck_datatype_line_info) +

  geom_bar(mapping = aes(x = reorder(primary_disease, -count_per_disease),

                         fill = data_type), stat = "count", position = "dodge") +

  xlab("Cancer Type") +

  ylab("Number of Cell Lines") +

  labs(fill = "Data Type") +

  # ggtitle(label = "Proportion of Cancer Types in DepMap Data (Bottleneck Subset)", subtitle = "By data type, 21Q2 Version - Overlap with CTRPv2: 79%") +

  theme(axis.text.x = element_text(angle = 45, hjust = 1),

        legend.position = c(.9,.55))

# ggsave(filename = "Plots/DepMap/DepMap_CTRP_Cell_Lines_Bottleneck_Proportion.pdf", device = "pdf")

# Compare bottleneck to original data

uniqueN(datatype_line_info$primary_disease)

uniqueN(datatype_line_info$stripped_cell_line_name)

uniqueN(bottleneck_datatype_line_info$primary_disease)

uniqueN(bottleneck_datatype_line_info$stripped_cell_line_name)

unique(bottleneck_datatype_line_info$stripped_cell_line_name)

ctrp[ccl_name %in% unique(bottleneck_datatype_line_info$stripped_cell_line_name)]

ctrp[ccl_name %in% unique(datatype_line_info$stripped_cell_line_name)] 

nrow(ctrp[ccl_name %in% unique(bottleneck_datatype_line_info$stripped_cell_line_name)]) /

nrow(ctrp[ccl_name %in% unique(datatype_line_info$stripped_cell_line_name)])   # 0.4038907

# Show cell line counts in bottleneck and original data side by side

bottleneck_datatype_line_info$`Data Subset` <- "Bottleneck"

datatype_line_info$`Data Subset` <- "Original"

both_data_type_line_info <- rbindlist(list(bottleneck_datatype_line_info, datatype_line_info))

# Change bar order

require(dplyr)

bar_level_df <- data.frame(x1 = c("Original", "Bottleneck"))

colnames(bar_level_df) <- "Data Subset"

both_data_type_line_info <- left_join(bar_level_df,  

                      both_data_type_line_info,

                      by = "Data Subset")

both_data_type_line_info <- as.data.table(both_data_type_line_info)

ggplot(data = both_data_type_line_info) +

  geom_bar(mapping = aes(x = reorder(primary_disease, -count_per_disease),

                         y = count_per_disease,

                         fill = `Data Subset`), stat = "identity", position = "dodge") +

  xlab("Cancer Type") +

  ylab("Number of Cell Line Profiles") +

  labs(fill = "Data Subset Type") +

  # scale_x_discrete(limits = c("Original", "Bottleneck")) +

  # ggtitle(label = "Proportion of Cancer Types in DepMap Data (Bottleneck Subset)", subtitle = "By data type, 21Q2 Version - Overlap with CTRPv2: 79%") +

  theme(axis.text.x = element_text(size = 14, angle = 45, hjust = 1, face = "bold"),

        axis.text.y = element_text(size = 14, face = "bold"),

        axis.title = element_text(size = 16, face = "bold"),

        legend.title = element_text(size = 16),

        legend.text = element_text(size=14),

        legend.position = c(.9,.55))

ggsave(filename = "Plots/DepMap/DepMap_CTRP_Cell_Lines_Bottleneck_Proportion_Comparison.pdf",

       device = "pdf")

# both_data_type_line_info <- merge(bottleneck_datatype_line_info, datatype_line_info, by)

bottleneck_datatype_line_info

# proteomics_datatype_line_info <- datatype_line_info[data_type == "Protein Quantification"]

# ggplot(data = proteomics_datatype_line_info) +

#   geom_bar(mapping = aes(x = reorder(primary_disease, -count_per_disease),

#                          fill = data_type), stat = "count", position = "dodge") +

#   xlab("Primary Disease") +