R/overlap_info.R

In KEGGlincs: Visualize all edges within a KEGG pathway and overlay LINCS data

#' Get overlap information for pairs of gene knock-outs from LINCS data
#' @export
#' @importFrom utils combn data
#' @import hgu133a.db
#' @import KOdata
#' @param KGML_file An object of formal class KEGGPathway
#' @param KEGG_mappings The data.frame object generated by the function 
#' expand_KEGG_mappings
#' @param cell_type Choose from the set of cell lines: 
#' (A375,A549,ASC,HA1E,HCC515,HEK293T,HEKTE,HEPG2,HT29,MCF7,NCIH716,NPC,PC3,
#' SHSY5Y,SKL,SW480,VCAP) 
#' @param data_type Choose from data types: (100_full, 100_bing, 50_lm)
#' @param pert_time Choose from (6,24,48,96,120,144,168)
#' @param only_mapped A logical indicator; if set to FALSE will return 'de-novo'
#' edges that 'exist' in data but are not documented in KEGG
#' @param affy_based A logical indicator; if set to TRUE will return 
#' lists/counts based on probeID instead of gene symbol.
#' @param keep_counts_only A logical indicator; if set to FALSE will return 
#' data frame with lists [of gene symbols or probe ids] as well as counts
#' @param add_fisher_information A logical indicator; by default the 
#' relationships are analyzed for strength of correlation via 
#' Fisher's Exact Test
#' @param p.adjust.method For available methods, type 'p.adjust.methods' into 
#' command promt and press enter.
#' @return A data frame where each row corresponds to information for pairs of 
#' experimental gene knock-outs from LINCS data (found in selected pathway).   
#' @examples 
#' p53_KGML <- get_KGML("hsa04115")
#' p53_KEGG_mappings <-  expand_KEGG_mappings(p53_KGML)
#' p53_edges <- expand_KEGG_edges(p53_KGML, p53_KEGG_mappings)
#' 
#' summary <- path_genes_by_cell_type(p53_KEGG_mappings)
#' p53_HA1E_data <- overlap_info(p53_KGML, p53_KEGG_mappings, 
#'                                "HA1E", data_type = "100_bing", 
#'                                only_mapped = FALSE)

overlap_info <- 
    function(KGML_file, KEGG_mappings, cell_type, 
                        data_type = "100_full",pert_time = 96,
                        only_mapped = TRUE, affy_based = FALSE,
                        keep_counts_only = TRUE, 
                        add_fisher_information = TRUE, 
                        p.adjust.method = "BH"){
    data("KO_data", envir = environment())
    KO_data <- get("KO_data")
    if(!cell_type %in% unique(KO_data$cell_id)){
        warning("Knock-out data is not available for selected cell-line \n")
        overlaps <- data.frame("null" = NA)
        return(overlaps)
        }
    KO_by_CT <- KO_data[KO_data$pert_time == pert_time & 
                        KO_data$cell_id == cell_type,]
    expanded_edges <- expand_KEGG_edges(KGML_file, KEGG_mappings)
    keeps <- c("entry1symbol", "entry2symbol")
    expanded_edges <- expanded_edges[expanded_edges$type != "maplink",keeps]
    if(nrow(expanded_edges) == 0 & only_mapped){
        warning("No documented edges are found in data; 
            only data for de-novo edges can be generated \n")
        overlaps <- data.frame("null" = NA)
        return(overlaps)
    }

    keeps <- c("entryACCESSION", "entrySYMBOL")
    path_genes <- KEGG_mappings[KEGG_mappings$entryTYPE == "gene", keeps]
    names(path_genes)[2] <- "SYMBOL"
    path_genes <- data.frame("ENTREZID" = unlist(path_genes$entryACCESSION), 
                            "SYMBOL" = unlist(path_genes$SYMBOL), 
                            stringsAsFactors = FALSE)
    path_genes <- path_genes[!duplicated(path_genes$ENTREZID),]
    data <- subset(path_genes, path_genes$SYMBOL %in% KO_by_CT$pert_desc)

    cat(paste0("Number of genes documented in selected pathway = ", 
                nrow(path_genes)), "\n")
    cat(paste0("Number of pathway genes in dataset = ", nrow(data),"\n"))
    cat(paste0("Coverage = ", round(nrow(data)/nrow(path_genes),4)*100,
                "%","\n"))
    if(nrow(data) == 0 | nrow(data) == 1){
        warning("Overlap data for selected pathway cannot be generated for 
            selected cell line")
        overlaps <- data.frame("null" = NA)
        return(overlaps)
    }

    conversion_key <- suppressWarnings(suppressMessages(
        AnnotationDbi::select(hgu133a.db::hgu133a.db, 
        AnnotationDbi::keys(hgu133a.db), c("SYMBOL","ENTREZID"), "PROBEID")))
  
    for (i in 1:nrow(data)){
        if (data_type == "100_full"){
            data$up[i] <- 
                strsplit(KO_by_CT$up100_full[KO_by_CT$pert_desc == 
                                                data$SYMBOL[i]], ";")
            data$down[i] <- 
                strsplit(KO_by_CT$dn100_full[KO_by_CT$pert_desc == 
                                                data$SYMBOL[i]], ";")
        }
        else if (data_type == "100_bing"){
            data$up[i] <- 
                strsplit(KO_by_CT$up100_bing[KO_by_CT$pert_desc == 
                                                data$SYMBOL[i]], ";")
            data$down[i] <- 
                strsplit(KO_by_CT$dn100_bing[KO_by_CT$pert_desc == 
                                                data$SYMBOL[i]], ";")
        }
        else if (data_type == "50_lm"){
            data$up[i] <- 
                strsplit(KO_by_CT$up50_lm[KO_by_CT$pert_desc == 
                                            data$SYMBOL[i]], ";")
            data$down[i] <- 
                strsplit(KO_by_CT$dn50_lm[KO_by_CT$pert_desc == 
                                            data$SYMBOL[i]], ";")
        }
        else {
            warning("valid data_type options: 100_full, 100_bing, 50_lm")
            return(NA)
        }
        data$up_symbol[i] <- 
            list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                            unlist(data$up[i]))])
        data$down_symbol[i] <-
            list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                            unlist(data$down[i]))])
    }
    overlaps<- data.frame(t(combn(data$SYMBOL,2)), stringsAsFactors = FALSE)
    names(overlaps) <- c("knockout1", "knockout2")
    overlaps$unique_ID <- paste0(overlaps$knockout1, ",", overlaps$knockout2)
    expanded_edges$unique_ID <- paste0(expanded_edges$entry1symbol, ",", 
                                    expanded_edges$entry2symbol)
    expanded_edges$unique_IDR <- paste0(expanded_edges$entry2symbol, ",", 
                                    expanded_edges$entry1symbol)
    pre_mapped1 <- subset(overlaps, overlaps$unique_ID %in% 
                        expanded_edges$unique_ID) ## Direction is correct
    pre_mapped2 <- subset(overlaps, overlaps$unique_ID %in% 
                        expanded_edges$unique_IDR) ## Direction is reversed
    pre_mapped2 <- pre_mapped2[,c(2,1)]
    names(pre_mapped2)[c(1,2)] = names(pre_mapped1)[c(1,2)]
    if(nrow(pre_mapped2) >= 1 & nrow(pre_mapped1) >= 1){
        pre_mapped2$unique_ID <- paste0(pre_mapped2$knockout1, ",", 
                                    pre_mapped2$knockout2)
        pre_mapped <- rbind(pre_mapped1, pre_mapped2)
    }
    if (nrow(pre_mapped2) == 0 & (nrow(pre_mapped1) == 0)){
        pre_mapped <- data.frame("unique_ID" = NA)
    }
    if (nrow(pre_mapped1) == 0 & nrow(pre_mapped2) >= 1){
        pre_mapped2$unique_ID <- paste0(pre_mapped2$knockout1, ",", 
                                    pre_mapped2$knockout2)
        pre_mapped <- pre_mapped2
    }
    if (nrow(pre_mapped2) == 0 & nrow(pre_mapped1) >= 1){
        pre_mapped <- pre_mapped1
    }
    if(is.na(pre_mapped[1,1]) & only_mapped){
        warning("No documented edges are found in data; only data for 
                de-novo edges can be generated \n")
        overlaps <- data.frame("null" = NA)
        return(overlaps)
    }
    pre_mapped$pre_mapped <- 1
    keeps <-c("knockout1", "knockout2", "pre_mapped")
    pre_mapped <- pre_mapped[,keeps]
    if (!only_mapped){
        un_mapped <- subset(overlaps, !overlaps$unique_ID %in% 
                        expanded_edges$unique_ID & !overlaps$unique_ID %in% 
                        expanded_edges$unique_IDR)
        un_mapped$pre_mapped <- 0
        un_mapped <- un_mapped[,keeps]
        overlaps <- rbind(pre_mapped, un_mapped)
    }
    else {
        overlaps <- pre_mapped
    }
    for (i in 1:nrow(overlaps)){
        overlaps$affy.genesUP[i] <- 
            list(intersect(unlist(data$up[which(data$SYMBOL == overlaps[i,1])]),
                        unlist(data$up[which(data$SYMBOL == overlaps[i,2])])))
    overlaps$affy.genesDOWN[i] <- 
        list(intersect(unlist(data$down[which(data$SYMBOL == overlaps[i,1])]), 
                        unlist(data$down[which(data$SYMBOL == overlaps[i,2])])))
    overlaps$affy.genesUPK1_DOWNK2[i] <- 
        list(intersect(unlist(data$up[which(data$SYMBOL == overlaps[i,1])]), 
                        unlist(data$down[which(data$SYMBOL == overlaps[i,2])])))
    overlaps$affy.genesDOWNK1_UPK2[i] <- 
        list(intersect(unlist(data$down[which(data$SYMBOL == overlaps[i,1])]), 
                        unlist(data$up[which(data$SYMBOL == overlaps[i,2])])))
    overlaps$num.affy.genesUP[i] <- length(unlist(overlaps$affy.genesUP[i]))
    overlaps$num.affy.genesDOWN[i] <- length(unlist(overlaps$affy.genesDOWN[i]))
    overlaps$num.affy.genesUPK1_DOWNK2[i] <- 
        length(unlist(overlaps$affy.genesUPK1_DOWNK2[i]))
    overlaps$num.affy.genesDOWNK1_UPK2[i] <- 
        length(unlist(overlaps$affy.genesDOWNK1_UPK2[i]))
    }
    if (!affy_based){
        for(i in 1:nrow(overlaps)){
            overlaps$genes.symbolsUP[i]<- 
                list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                        unlist(overlaps$affy.genesUP[i]))])
            overlaps$genes.symbolsDOWN[i]<- 
                list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                        unlist(overlaps$affy.genesDOWN[i]))])
        overlaps$genes.symbolsUPK1_DOWNK2[i]<- 
            list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                unlist(overlaps$affy.genesUPK1_DOWNK2[i]))])
        overlaps$genes.symbolsDOWNK1_UPK2[i]<- 
            list(conversion_key$SYMBOL[which(conversion_key$PROBEID %in% 
                                unlist(overlaps$affy.genesDOWNK1_UPK2[i]))])
        overlaps$num.genes.symbolsUP[i] <- 
            length(unlist(overlaps$genes.symbolsUP[i]))
        overlaps$num.genes.symbolsDOWN[i] <- 
            length(unlist(overlaps$genes.symbolsDOWN[i]))
        overlaps$num.genes.symbolsUPK1_DOWNK2[i] <- 
            length(unlist(overlaps$genes.symbolsUPK1_DOWNK2[i]))
        overlaps$num.genes.symbolsDOWNK1_UPK2[i] <- 
            length(unlist(overlaps$genes.symbolsDOWNK1_UPK2[i]))
        }
        keeps = c("knockout1","knockout2","num.genes.symbolsUP", 
                    "num.genes.symbolsDOWN", "num.genes.symbolsUPK1_DOWNK2", 
                    "num.genes.symbolsDOWNK1_UPK2", "genes.symbolsUP", 
                    "genes.symbolsDOWN", "genes.symbolsUPK1_DOWNK2", 
                    "genes.symbolsDOWNK1_UPK2", "pre_mapped")
        overlaps <- overlaps[,keeps]
    }
    else {
        keeps = c("knockout1","knockout2","num.affy.genesUP", 
                "num.affy.genesDOWN", "num.affy.genesUPK1_DOWNK2", 
                "num.affy.genesDOWNK1_UPK2", "affy.genesUP", "affy.genesDOWN", 
                "affy.genesUPK1_DOWNK2", "affy.genesDOWNK1_UPK2", "pre_mapped")
        overlaps <- overlaps[, keeps]
    }
    names(overlaps)[c(3:6)] <- c("UP", "DOWN", "UK1_DK2", "DK1_UK2")
    method = p.adjust.method
    if (keep_counts_only){
        keeps = c("knockout1","knockout2", "UP", "DOWN", "UK1_DK2", "DK1_UK2", 
                "pre_mapped")
        overlaps <- overlaps[,keeps]
        }
    if(add_fisher_information){
        overlaps <- get_fisher_info(overlaps, method)
    }
    return(overlaps)
}