R/pwOmics_vizualisation_functions.R

In pwOmics: Pathway-based data integration of omics data

#' Plot time profile clusters of dynamic analysis result.
#' 
#' @param fuzzed_matsplines result of dynamic analysis: inferred net generated 
#' by consDynamicNet function.
#' @param ... further plotting/legend parameters.
#' @return pdf file in current working directory.
#' 
#' @keywords manip
#' @export
#' @examples
#' \dontrun{
#' data(OmicsExampleData)
#' data_omics = readOmics(tp_prots = c(0.25, 1, 4, 8, 13, 18, 24), 
#' tp_genes = c(1, 4, 8, 13, 18, 24), OmicsExampleData,
#' PWdatabase = c("biocarta", "kegg", "nci", "reactome"), 
#' TFtargetdatabase = c("userspec"))
#' data_omics = readPhosphodata(data_omics, 
#' phosphoreg = system.file("extdata", "phospho_reg_table.txt", 
#' package = "pwOmics")) 
#' data_omics = readTFdata(data_omics, 
#' TF_target_path = system.file("extdata", "TF_targets.txt", 
#' package = "pwOmics"))
#' data_omics_plus = readPWdata(data_omics,  
#' loadgenelists = system.file("extdata/Genelists", package = "pwOmics")) 
#' 
#' data_omics_plus = identifyPR(data_omics_plus)
#' setwd(system.file("extdata/Genelists", package = "pwOmics"))
#' data_omics = identifyPWs(data_omics_plus)
#' data_omics = identifyTFs(data_omics)
#' data_omics = identifyRsofTFs(data_omics, 
#' noTFs_inPW = 1, order_neighbors = 10)
#' data_omics = identifyPWTFTGs(data_omics)
#' statConsNet = staticConsensusNet(data_omics)
#' consDynNet = consDynamicNet(data_omics, statConsNet)
#' timeprof = clusterTimeProfiles(consDynNet)
#' plotTimeProfileClusters(timeprof)
#' }
plotTimeProfileClusters <- function(fuzzed_matsplines) {
    
    requireNamespace("Mfuzz", quietly = TRUE)
    pdf(paste(getwd(), "/", "TimeProfiles.pdf", sep = ""))
    par(mar=c(5.1, 4.1, 4.1, 8.1), xpd = TRUE )
    plot(x = NULL, y = NULL, bty = 'L', 
         xlim = c(min(as.numeric(colnames(fuzzed_matsplines$centers))), 
                  max(as.numeric(colnames(fuzzed_matsplines$centers)))+1),
         ylim = c(min(fuzzed_matsplines$centers),max(fuzzed_matsplines$centers)), 
         xlab = "time", ylab = "normalized expression", 
         main = paste("Fuzzy c-means clustering\n with ", 
                      max(fuzzed_matsplines$cluster), " centers", sep = ""))
    colors = rainbow(max(fuzzed_matsplines$cluster))
    ind_sort = order(fuzzed_matsplines$cluster)
    
    for(s in 1: max(fuzzed_matsplines$cluster))
    {lines(as.numeric(colnames(fuzzed_matsplines$centers)), 
           fuzzed_matsplines$centers[s,], col = colors[s])}
    if(length(names(fuzzed_matsplines$cluster)[ind_sort]) < 20)
    {legend("topright", inset = c(-0.3,0), 
            legend = names(fuzzed_matsplines$cluster)[ind_sort], 
            fill = colors[fuzzed_matsplines$cluster][ind_sort], cex = 0.7)
    }else{legend("topright", inset = c(-0.3,0), 
                 legend = paste("cluster ", 1:dim(fuzzed_matsplines$centers)[1],
                                sep = ""), 
                 fill = colors[1:dim(fuzzed_matsplines$centers)[1]], cex = 0.7)        
    }
    dev.off()
}

#' Plot inferred net based on analysis analysis.
#' 
#' Dynamic analysis result is plotted to pdf file stored in current working 
#' directory.
#'  
#' @param dynConsensusNet result of dynamic analysis: inferred net generated by 
#' consDynamicNet function.
#' @param sig.level significance level for plotting edges in network 
#' (between 0 and 1).
#' @param clarify indicating if unconnected nodes should be removed; 
#' default = "TRUE".
#' @param layout igraph layout parameter; default is 
#' layout.fruchterman.reingold.
#' @param ... further plotting/legend parameters.
#' @return pdf file in current working directory.
#' 
#' @keywords manip
#' @export
#' @examples
#' \dontrun{
#' data(OmicsExampleData)
#' data_omics = readOmics(tp_prots = c(0.25, 1, 4, 8, 13, 18, 24), 
#' tp_genes = c(1, 4, 8, 13, 18, 24), OmicsExampleData,
#' PWdatabase = c("biocarta", "kegg", "nci", "reactome"), 
#' TFtargetdatabase = c("userspec"))
#' data_omics = readPhosphodata(data_omics, 
#' phosphoreg = system.file("extdata", "phospho_reg_table.txt", 
#' package = "pwOmics")) 
#' data_omics = readTFdata(data_omics, 
#' TF_target_path = system.file("extdata", "TF_targets.txt", 
#' package = "pwOmics"))
#' data_omics_plus = readPWdata(data_omics,  
#' loadgenelists = system.file("extdata/Genelists", package = "pwOmics")) 
#' 
#' data_omics_plus = identifyPR(data_omics_plus)
#' setwd(system.file("extdata/Genelists", package = "pwOmics"))
#' data_omics = identifyPWs(data_omics_plus)
#' data_omics = identifyTFs(data_omics)
#' data_omics = identifyRsofTFs(data_omics, 
#' noTFs_inPW = 1, order_neighbors = 10)
#' data_omics = identifyPWTFTGs(data_omics)
#' statConsNet = staticConsensusNet(data_omics)
#' dynConsNet = consDynamicNet(data_omics, statConsNet)
#' plotConsDynNet(dynConsNet, sig.level = 0.8)
#' }
plotConsDynNet <- function(dynConsensusNet, sig.level, clarify = "TRUE",
                           layout = layout.fruchterman.reingold, ...) {
    
    dynConsNet = dynConsensusNet[[1]]
    pdf(paste(getwd(), "/", "dynConsensusNet.pdf", sep = ""))
    dynnet_z = t(dynConsNet$z)
    dynnet <- matrix(0, nrow = nrow(dynnet_z), ncol = ncol(dynnet_z))
    colnames(dynnet) <- colnames(dynnet_z)
    rownames(dynnet) <- rownames(dynnet_z)
    cutoff <- qnorm((1 + sig.level)/2)
    if(dim(which(abs(dynnet_z) > cutoff, arr.ind = TRUE))[1] == 0)
    { stop("No nodes in the network have a significance higher than the
           selected value. Please choose a less strict significance level.\n")}
    dynnet[which(abs(dynnet_z) > cutoff, arr.ind = TRUE)] <- 1
    remove_ind <- which(colSums(dynnet) == 0 & rowSums(dynnet) == 0)
    if(length(remove_ind) >= dim(dynnet)[1]-1)
    { stop("Please select clarify = 'FALSE' to have network contents visualized
           and/or select a less strict significance level.\n")}
    if (clarify == "TRUE" & length(remove_ind) > 0) 
    {dynnet <- dynnet[-remove_ind, -remove_ind]}
    dyngraph <- graph.adjacency(dynnet, mode = c("directed"), 
                                add.rownames = TRUE)
    ident = unlist(strsplit(as.character(V(dyngraph)$name),"_"))[seq(2, length(V(dyngraph)$name)*2, by = 2)]
    V(dyngraph)$color = rep("grey", times = length(V(dyngraph)))
    V(dyngraph)$color[which(ident == "p")] = rep("red", times = length(which(ident == "p")))
    V(dyngraph)$color[which(ident == "g")] = rep("green", times = length(which(ident == "g")))
    V(dyngraph)$name = unlist(strsplit(V(dyngraph)$name, "_"))[seq(1, length(ident)*2, by = 2)]
    V(dyngraph)$label.cex = 0.6
    V(dyngraph)$label.color = "black"
    E(dyngraph)$color = rep("grey", times = length(E(dyngraph)))
    E(dyngraph)$color = "red4"
    E(dyngraph)$color[which(E(dyngraph)$weight==1)] = "green4"
        plot(dyngraph, main = "Dynamic consensus net", edge.arrow.size=0.5, layout = layout,...)
    legend("topleft" , legend = c("consensus proteins", "consensus genes"), 
           fill = c("red", "green"), cex = 0.7, ...)
    dev.off()
}





#' Plot consensus graph(s) from static analysis.
#' 
#' Consensus graph(s) plotted to pdf file stored in current working directory.
#'  
#' @param consensusGraphs result from static analysis: consensus graph generated 
#' by staticConsensusNet function.
#' @param data_omics OmicsData object.
#' @param ... further plotting/legend parameters.
#' @return pdf file in current working directory.
#' 
#' @keywords manip
#' @export
#' @examples
#' \dontrun{
#' data(OmicsExampleData)
#' data_omics = readOmics(tp_prots = c(0.25, 1, 4, 8, 13, 18, 24), 
#' tp_genes = c(1, 4, 8, 13, 18, 24), OmicsExampleData,
#' PWdatabase = c("biocarta", "kegg", "nci", "reactome"), 
#' TFtargetdatabase = c("userspec"))
#' data_omics = readPhosphodata(data_omics, 
#' phosphoreg = system.file("extdata", "phospho_reg_table.txt", 
#' package = "pwOmics")) 
#' data_omics = readTFdata(data_omics, 
#' TF_target_path = system.file("extdata", "TF_targets.txt", 
#' package = "pwOmics"))
#' data_omics_plus = readPWdata(data_omics,  
#' loadgenelists = system.file("extdata/Genelists", package = "pwOmics")) 
#' 
#' data_omics_plus = identifyPR(data_omics_plus)
#' setwd(system.file("extdata/Genelists", package = "pwOmics"))
#' data_omics = identifyPWs(data_omics_plus)
#' data_omics = identifyTFs(data_omics)
#' data_omics = identifyRsofTFs(data_omics, 
#' noTFs_inPW = 1, order_neighbors = 10)
#' data_omics = identifyPWTFTGs(data_omics)
#' statConsNet = staticConsensusNet(data_omics)
#' plot(ConsensusGraph, data_omics)
#' }
plotConsensusGraph <- function(consensusGraphs, data_omics, ...) {
    
    if(class(data_omics) != "OmicsData")
    {stop("Parameter 'data_omics' is not an OmicsData object.")}
    
    if(TRUE %in% !(names(consensusGraphs) %in% 
                       as.character(data_omics[[1]][[1]][[1]][[1]]) &
                       names(consensusGraphs) %in% 
                       as.character(data_omics[[1]][[1]][[1]][[2]])) ) 
    {warning("Warning: ConsensusGraphs do not match OmicsData object.")}
    
    pdf(paste(getwd(), "/", "ConsensusGraphs.pdf", sep = ""))
    same_tps = data_omics[[1]][[1]][[1]][[1]][which(data_omics[[1]][[1]][[1]][[1]] 
                                                    %in% data_omics[[1]][[1]][[1]][[2]])]
    for(k in 1: length(consensusGraphs))
    {
        V(consensusGraphs[[k]])$label.cex = 0.6
        plot.igraph(consensusGraphs[[k]], 
                    main = paste("Consensus graph\n time ", 
                                 same_tps[k], sep = ""), vertex.size = 7, 
                    vertex.label.color = "black", vertex.label.dist = 0.4, 
                    vertex.label.font = 2, ...)
        legend("topleft" , legend = c("consensus proteins", 
                                      "steiner node proteins", "consensus TFs", "consensus target genes"), 
               fill = c("red", "yellow", "lightblue", "green"), cex = 0.7, ...)
    }
    dev.off()
}


#' Plot consensus graph profiles of static consensus molecules.
#' 
#' Consensus graph profiles of static consensus molecules 
#' plotted as heatmap to pdf file stored in current working directory.
#'  
#' @param consensusGraphs result from static analysis: consensus graph generated 
#' by staticConsensusNet function.
#' @param data_omics OmicsData object.
#' @param subsel character vector of selected consensus molecules for plotting; 
#' if TRUE all consensus molecules are plotted
#' @param ... further plotting/legend parameters.
#' @return pdf file in current working directory.
#' 
#' @keywords manip
#' @export
#' @examples
#' \dontrun{
#' data(OmicsExampleData)
#' data_omics = readOmics(tp_prots = c(0.25, 1, 4, 8, 13, 18, 24), 
#' tp_genes = c(1, 4, 8, 13, 18, 24), OmicsExampleData,
#' PWdatabase = c("biocarta", "kegg", "nci", "reactome"), 
#' TFtargetdatabase = c("userspec"))
#' data_omics = readPhosphodata(data_omics, 
#' phosphoreg = system.file("extdata", "phospho_reg_table.txt", 
#' package = "pwOmics")) 
#' data_omics = readTFdata(data_omics, 
#' TF_target_path = system.file("extdata", "TF_targets.txt", 
#' package = "pwOmics"))
#' data_omics_plus = readPWdata(data_omics,  
#' loadgenelists = system.file("extdata/Genelists", package = "pwOmics")) 
#' 
#' data_omics_plus = identifyPR(data_omics_plus)
#' setwd(system.file("extdata/Genelists", package = "pwOmics"))
#' data_omics = identifyPWs(data_omics_plus)
#' data_omics = identifyTFs(data_omics)
#' data_omics = identifyRsofTFs(data_omics, 
#' noTFs_inPW = 1, order_neighbors = 10)
#' data_omics = identifyPWTFTGs(data_omics)
#' statConsNet = staticConsensusNet(data_omics)
#' plotConsensusProfiles(statConsNet, data_omics, subsel = TRUE)
#' }
plotConsensusProfiles <- function(consensusGraphs, data_omics, subsel = TRUE, ...) {
    
    if(class(data_omics) != "OmicsData")
    {stop("Parameter 'data_omics' is not an OmicsData object.")}
    
    if(TRUE %in% !(names(consensusGraphs) %in% 
                   as.character(data_omics[[1]][[1]][[1]][[1]]) &
                   names(consensusGraphs) %in% 
                   as.character(data_omics[[1]][[1]][[1]][[2]])) ) 
    {warning("Warning: ConsensusGraphs do not match OmicsData object.")}
    
    same_tps = data_omics[[1]][[1]][[1]][[1]][which(data_omics[[1]][[1]][[1]][[1]] 
                                                    %in% data_omics[[1]][[1]][[1]][[2]])]
    Cons_graph_members = vector()
    for(k in 1: length(same_tps))
    {Cons_graph_members = c(Cons_graph_members, unique(c(V(consensusGraphs[[k]])$name))) }
    Cons_graph_members = unique(Cons_graph_members)
    
    heat_CGM = matrix(nrow = length(Cons_graph_members), ncol = length(same_tps),
                      dimnames = list(Cons_graph_members, same_tps))
    for(s in 1: length(Cons_graph_members))
    {  for(h in 1: length(same_tps))
    {heat_CGM[s,h] = if(Cons_graph_members[s] %in% V(consensusGraphs[[h]])$name)
    {V(consensusGraphs[[h]])$color[which(V(consensusGraphs[[h]])$name == Cons_graph_members[s])]
    }else{"white"}
    }
    }
    colvec = c("white", "green", "lightblue", "red", "yellow")
    numvec = c(0,1,2,3,4)
    for(s in 1:5)
    {heat_CGM[which(heat_CGM == colvec[s])] = numvec[s]}
    heat_CGM = apply(heat_CGM,2, as.numeric)
    rownames(heat_CGM) = Cons_graph_members
    
    if((subsel == TRUE)[1])
    {ind_sel = seq(1, length(Cons_graph_members))
    }else{
        ind_sel = which(rownames(heat_CGM) %in% subsel)  
        if(!5 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1) || 
           (!5 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1) & !4 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1)) ||
           (!5 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1) &  !4 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1) &
            !3 %in% ((unique(as.vector(heat_CGM[ind_sel,])))+1) ))
        {colvec = colvec[sort(unique(as.vector(heat_CGM[ind_sel,])))+1]}
    }
    pdf(paste(getwd(), "/", "ConsensusGraphs.pdf", sep = ""))
    heatmap.2(heat_CGM[ind_sel,], Colv = NA, dendrogram = "row", col = colvec,
              scale = "none", main = "Consensus graph profiles", trace = "none", key = FALSE, ...)
    legend("topleft", fill = c("red", "yellow", "lightblue", "green"), 
           legend = c("consensus proteins", "steiner node proteins", "consensus TFs", "consensus target genes"), cex = 0.7)
    dev.off()
}



#' Plot temporal correlations of phosphoprotein expression levels and affected 
#' downstream transcripts.
#' 
#' This function plots an overview of consensus phosphoprotein expression level 
#' correlations with those transcripts that are affected downstream in a newly
#' generated folder in the working directory. The following 
#' additional information needs to be provided: Phosphorylation information amino
#' acid, position and multiplicity with the expression levels for each
#' of the time points the correlations should be plotted for. Furthermore 
#' data tables for fold changes/ratios of phosphoproteins and transcripts that 
#' are part of the data_omics object.
#' 
#' @param ConsensusGraph result from static analysis: consensus graph generated 
#' by staticConsensusNet function.
#' @param timepointsprot numeric vector with measurement time points in 
#' phosphoproteome data set, which should be used for correlation plots
#' @param timepointstrans numeric vector with measurement time points in
#' transcriptome data set, which should be used for correlation plots
#' @param foldername character vector specifying the name of the folder that
#' will be generated for the temporal correlations
#' @param trans_sign character vector specifying a tab-delimited file with the
#' transcriptome expression levels for all time points in timepointstrans
#' @param trans_sign_names character vector specifying column names in the
#' transcriptome file corresponding to the expression levels at different time
#' points.  
#' @param phospho_sign character vector specifying a tab-delimited file with the
#' phosphoproteome information (columns 'Gene.names', 'Amino.acid', 'Position',
#' 'Multiplicity') and expression levels for all time points in timepointstrans
#' @param phospho_sign_names character vector specifying column names in the
#' phosphoproteome file corresponding to the expression levels at different time
#' points. 
#' @param ... further plotting/legend parameters.
#' @return pdf file in current working directory.
#' 
#' @keywords manip
#' @export
#' @examples
#' \dontrun{
#' data(OmicsExampleData)
#' data_omics = readOmics(tp_prots = c(0.25, 1, 4, 8, 13, 18, 24), 
#' tp_genes = c(1, 4, 8, 13, 18, 24), OmicsExampleData,
#' PWdatabase = c("biocarta", "kegg", "nci", "reactome"), 
#' TFtargetdatabase = c("userspec"))
#' data_omics = readPhosphodata(data_omics, 
#' phosphoreg = system.file("extdata", "phospho_reg_table.txt", 
#' package = "pwOmics")) 
#' data_omics = readTFdata(data_omics, 
#' TF_target_path = system.file("extdata", "TF_targets.txt", 
#' package = "pwOmics"))
#' data_omics_plus = readPWdata(data_omics,  
#' loadgenelists = system.file("extdata/Genelists", package = "pwOmics")) 
#' 
#' data_omics_plus = identifyPR(data_omics_plus)
#' setwd(system.file("extdata/Genelists", package = "pwOmics"))
#' data_omics = identifyPWs(data_omics_plus)
#' data_omics = identifyTFs(data_omics)
#' data_omics = identifyRsofTFs(data_omics, 
#' noTFs_inPW = 1, order_neighbors = 10)
#' data_omics = identifyPWTFTGs(data_omics)
#' statConsNet = staticConsensusNet(data_omics)
#' temp_correlations(statConsNet, timepointsprot = c(1,4,8,13,18,24),
#' timepointstrans = c(1,4,8,13,18,24), foldername = "ProtCons_", 
#' trans_sign = system.file("extdata", "signif_single.csv", package = "pwOmics") 
#' trans_sign_names = c("FC_1",  "FC_2", "FC_3", "FC_4"), 
#' phospho_sign = system.file("extdata", "mat_phospho.csv", package = "pwOmics") 
#' phospho_sign_names = c("Rat1", "Rat2", "Rat3", "Rat4")) )
#' }
#' 
temp_correlations <- function (ConsensusGraph, timepointsprot, 
                                   timepointstrans, foldername = "ProtCons_", 
                                   trans_sign = "signif_single.csv", 
                                   trans_sign_names, 
                                   phospho_sign = "mat_phospho.csv", 
                                   phospho_sign_names) 
{
    trans_sign = read.csv(trans_sign, sep = "\t")
    phospho_sign = read.csv(phospho_sign, sep = ",")
    proteins = names(V(ConsensusGraph)[which(V(ConsensusGraph)$color == "red")])
    filenames = list.files(getwd(), pattern = "_downstream")
    filenames_prot = filenames[which(gsub("_downstream", "", filenames) %in% proteins)]
    for (g in 1:length(filenames_prot)) {
        print(filenames_prot[g])
        pdf(paste(foldername, filenames_prot[g], ".pdf", sep = ""))
        temp = NULL
        filenames = list.files(path = paste(getwd(), filenames_prot[g], 
                                            sep = "/"), pattern = "transcripts.csv")
        filenamesR = list.files(path = paste(getwd(), filenames_prot[g], 
                                             sep = "/"), pattern = ".RData")
        ind_vec = vector()
        for (k in 1:length(timepointsprot)) {
            ind_vec[k] = which(grepl(paste("tp", as.character(timepointsprot[k]), 
                                           "_", sep = ""), filenames) == TRUE)
        }
        times_cov = vector(length = length(timepointsprot))
        for (f in 1:length(timepointsprot)) {
            load(paste(getwd(), filenames_prot[g], filenamesR[ind_vec[f]], 
                       sep = "/"))
            if (length(axis) != 0) {
                temp = read.csv(paste(getwd(), filenames_prot[g], 
                                      filenames[ind_vec[f]], sep = "/"))
                times_cov[f] = TRUE
            }
        }
        names(times_cov) = paste("tp_", timepointsprot, sep = "")
        if (length(temp) > 0) {
            temp[, 1] = NULL
            transcripts_to_corr = as.character(temp[, 1])
            data_transcripts = trans_sign[match(transcripts_to_corr, 
                                            as.character(trans_sign[, 1])), ]
            data_phosphoprot = phospho_sign[grep(paste(gsub("_downstream", "", 
                                    filenames_prot[g]), "$", sep = ""), 
                                    phospho_sign$T..Gene.names), 1:14]
            
            colfunc <- "blue"
            par(mfrow = c(3, 2), oma = c(1, 1, 1, 1), mar = c(5, 4, 4, 2) + 0.1)
            plot(1, type="n", axes=FALSE, xlab="", ylab="")
            legend("center", legend = paste(timepointsprot, " min / ", 
                                            timepointstrans, " min", sep = ""), col = "black", 
                   pch = c(19, 18, 17, 15), title = "phosphoproteome / transcriptome time points", 
                   bty = "n")
            if(dim(data_phosphoprot)[1]>1)
            {data_phosphoprot_new = apply(data_phosphoprot[,phospho_sign_names],2,as.numeric)
            rownames(data_phosphoprot_new) = paste(data_phosphoprot[,"T..Gene.names"], 
                                                   "_", data_phosphoprot[,"C..Amino.acid"], 
                                                   data_phosphoprot[,"N..Position"], "_",
                                                   gsub("___", "M", data_phosphoprot[,"C..Multiplicity"]), sep = "")
            colh = colorRampPalette(c("blue", "orange"))(dim(data_phosphoprot_new)[1])
            for (j in 1:dim(data_transcripts)[1]) {
                plot(NA, NA, xlim = c(min(data_transcripts[, trans_sign_names], 
                                    na.rm = TRUE) - 0.4, max(data_transcripts[, 
                                    trans_sign_names], na.rm = TRUE) + 2), 
                     ylim = c(min(data_phosphoprot_new, na.rm = TRUE) - 0.1, 
                              max(data_phosphoprot_new, na.rm = TRUE) + 0.1), xlab = "FC downstream transcript", 
                     ylab = paste("log2 ratio ", gsub("_downstream", "", 
                             filenames_prot[g]), sep = ""), main = data_transcripts[j, 1])
                
                for(h in 1: dim(data_phosphoprot_new)[1])
                {   for (s in 1:length(timepointsprot) - 1) {
                    lines(as.numeric(data_transcripts[j, trans_sign_names[s:(s + 1)]]), 
                          data_phosphoprot_new[h,s:(s + 1)], col = colh[h], lwd = 1)}
                    points(data_transcripts[j, trans_sign_names], 
                           data_phosphoprot_new[h,], col = colh[h], 
                           pch = c(19, 18, 17, 15), cex = 1.2)
                }
                legend("bottomright", fill = colh, legend = rownames(data_phosphoprot_new), cex = 0.7)
            }
            }else{
                data_phosphoprot_new = apply(data_phosphoprot[,phospho_sign_names],2,as.numeric)
                nam = paste(data_phosphoprot[,"T..Gene.names"], "_", 
                            data_phosphoprot[,"C..Amino.acid"], 
                            data_phosphoprot[,"N..Position"], "_",
                            gsub("___", "M", data_phosphoprot[,"C..Multiplicity"]), sep = "")
                for (j in 1:dim(data_transcripts)[1]) {
                    plot(NA, NA, xlim = c(min(data_transcripts[, trans_sign_names], 
                                              na.rm = TRUE) - 0.4, 
                                          max(data_transcripts[, trans_sign_names], 
                                              na.rm = TRUE) + 2), 
                         ylim = c(min(data_phosphoprot_new, na.rm = TRUE) - 0.1, 
                                  max(data_phosphoprot_new, na.rm = TRUE) + 0.1), 
                         xlab = "FC downstream transcript", 
                         ylab = paste("log2 ratio ", gsub("_downstream", "", 
                                filenames_prot[g]), sep = ""), main = data_transcripts[j, 1])
                    for (s in 1:length(timepointsprot) - 1) {
                        lines(as.numeric(data_transcripts[j, trans_sign_names[s:(s + 1)]]), 
                              data_phosphoprot_new[s:(s + 1)], col = colfunc, lwd = 1)
                    }
                    points(data_transcripts[j, trans_sign_names],  data_phosphoprot_new, 
                           col = colfunc, pch = c(19, 18, 17, 15), cex = 1.2)
                    legend("bottomright", fill = colfunc, legend = nam, cex = 0.7)
                }
            }
        }
        dev.off()
    }
}