| Previous changeset 24:64b822045467 (2016-01-29) Next changeset 26:8371404ab073 (2016-01-29) |
|
Commit message:
Uploaded |
|
added:
bubbles_v9_NSAF_natural_log.R |
| b |
| diff -r 64b822045467 -r ab602bbf4ac5 bubbles_v9_NSAF_natural_log.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bubbles_v9_NSAF_natural_log.R Fri Jan 29 09:33:58 2016 -0500 |
| [ |
| b'@@ -0,0 +1,268 @@\n+rm(list = ls()) \r\n+###################################################################################################\r\n+# R-code: Multi-bubble graph generation from SAINTexpress output\r\n+# Author: Brent Kuenzi\r\n+###################################################################################################\r\n+# This Script generates the bubble graphs based upon Saint output. \r\n+###################################################################################################\r\n+# Copyright (C) Brent Kuenzi.\r\n+# Permission is granted to copy, distribute and/or modify this document\r\n+# under the terms of the GNU Free Documentation License, Version 1.3\r\n+# or any later version published by the Free Software Foundation;\r\n+# with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.\r\n+# A copy of the license is included in the section entitled "GNU\r\n+# Free Documentation License".\r\n+###################################################################################################\r\n+## REQUIRED INPUT ##\r\n+\r\n+# 1) listfile: SAINTexpress generated "list.txt" file\r\n+# 2) preyfile: SAINT pre-processing generated "prey.txt" file used to run SAINTexpress\r\n+\r\n+## OPTIONAL INPUT ##\r\n+\r\n+# 3) crapome: raw output from crapome Workflow 1 query (http://www.crapome.org)\r\n+# 4) color: bubble color (default = "red")\r\n+# - color = "crapome": color bubbles based on Crapome(%)\r\n+# - Also recognizes any color within R\'s built-in colors() vector\r\n+# 5) label: Adds gene name labels to bubbles within the "zoomed in" graphs (default = FALSE)\r\n+# 6) cutoff: Saintscore cutoff to be assigned for filtering the "zoomed in" graphs (default = 0.8)\r\n+# 7) type: Specifies if the data is MaxQuant (MQ) or Scaffold (SC) data (default = "SC")\r\n+# 8) inc_file: Selects only the uniprot ids in the provided list (default ="None")\r\n+# 9) exc_file: Removes the proteins in the list (default = "None")\r\n+###################################################################################################\r\n+\r\n+\r\n+ins_check_run <- function(){\r\n+ if (\'dplyr\' %in% rownames(installed.packages())){}\r\n+ else {\r\n+ install.packages(\'dplyr\', repos = \'http://cran.us.r-project.org\')\r\n+ }\r\n+ if (\'tidyr\' %in% rownames(installed.packages())){}\r\n+ else {\r\n+ install.packages(\'tidyr\', repos = \'http://cran.us.r-project.org\')\r\n+ }\r\n+ if (\'ggplot2\' %in% rownames(installed.packages())){}\r\n+ else {\r\n+ install.packages(\'ggplot2\', repos = \'http://cran.us.r-project.org\')\r\n+ }\r\n+}\r\n+\r\n+ins_check_run()\r\n+library(dplyr); library(tidyr); library(ggplot2)\r\n+\r\n+\r\n+main <- function(listfile, preyfile, crapome = FALSE, color = "red", label = FALSE, cutoff = 0.8, type = "SC", inc_file = "None", exc_file = "None" ) {\r\n+ cutoff_check(cutoff)\r\n+ listfile <- list_type(listfile, inc_file, exc_file)\r\n+ if(type == "SC") {\r\n+ df <- merge_files_sc(listfile, preyfile, crapome)\r\n+ }\r\n+ if(type == "MQ") {\r\n+ df <- merge_files_mq(listfile, preyfile, crapome)\r\n+ }\r\n+ bubble_NSAF(df, color)\r\n+ bubble_SAINT(df, color)\r\n+ bubble_zoom_SAINT(df, color, label, cutoff)\r\n+ bubble_zoom_NSAF(df, color, label, cutoff)\r\n+ write.table(df, "output.txt", sep = "\\t", quote = FALSE, row.names = FALSE)\r\n+}\r\n+###################################################################################################\r\n+# Include and Exclude list filtering\r\n+###################################################################################################\r\n+list_type <- function(df, inc_file, exc_file) {\r\n+ Saint <- read.delim(df, stringsAsFactors = FALSE)\r\n+ if (inc_file != "None") {\r\n+ if (exc_file == "None") {\r\n+ inc_prots <- read.delim(inc_file, sep = \'\\t\', header = FALSE, stringsAsFactors = FALSE)\r\n+ print(inc_prots[, 1])\r\n+ print(Saint$Prey)\r\n+ filtered_df = subset(Saint, Saint$Prey == inc_prots[, 1])\r\n+ }\r\n+ else {\r\n+ inc_prots <- read.delim(inc_file, sep = \'\\t\', header = FALSE, stringsAsFactors = FALSE)\r\n+ exc_prots <- read'..b'yGene))\r\n+ p <- qplot(x = SaintScore, y = log2(FoldChange), data = a, colour = I(color), size = SpecSum) +\r\n+ scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + \r\n+ geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)\r\n+ if(label == TRUE & length(a$NSAF != 0)) {\r\n+ p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)\r\n+ }\r\n+ if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}\r\n+ return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_SAINT.png"))\r\n+ }\r\n+}\r\n+###################################################################################################\r\n+# Filter proteins on Saintscore cutoff and plot for each bait x = log(NSAF), y = Log2(FoldChange)\r\n+###################################################################################################\r\n+bubble_zoom_NSAF <- function(data, color, label = FALSE, cutoff = 0.8) {\r\n+ if(color == "crapome") {\r\n+ a <- subset(data, CrapomePCT < 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))\r\n+ b <- subset(data, CrapomePCT >= 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))\r\n+ p <- qplot(x = log(NSAF), y = log2(FoldChange), data = a, colour = I("tan"), size = SpecSum) + \r\n+ scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + \r\n+ geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)\r\n+ if(label == TRUE & length(a$NSAF != 0)) {\r\n+ p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black")\r\n+ }\r\n+ if(length(levels(a$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}\r\n+ p <- p + geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum, color = CrapomePCT), data = b) + \r\n+ scale_colour_gradient(limits = c(80, 100), low = "tan", high = "red") + \r\n+ labs(colour = "CRAPome Probability \\nof Specific Interaction (%)", x = "ln(NSAF)") + \r\n+ geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = b)\r\n+ if(label == TRUE & length(b$NSAF != 0)) {\r\n+ p <- p + geom_text(data = b, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)\r\n+ }\r\n+ return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_NSAF.png"))\r\n+ }\r\n+ if(color != "crapome") {\r\n+ a <- subset(data, SaintScore >= cutoff, select = c(NSAF, SpecSum, FoldChange, SaintScore, Bait, PreyGene))\r\n+ p <- qplot(x = log(NSAF), y = log2(FoldChange), data = a, colour = I(color), size = SpecSum) +\r\n+ scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + \r\n+ geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a) + \r\n+ labs(x = "ln(NSAF)")\r\n+ if(label == TRUE & length(a$NSAF != 0)) {\r\n+ p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)\r\n+ }\r\n+ if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}\r\n+ return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_NSAF.png"))\r\n+ }\r\n+}\r\n+###################################################################################################\r\n+# Check Saintscore cutoff and stop program if not between 0 and 1\r\n+###################################################################################################\r\n+cutoff_check <- function(cutoff){\r\n+ if( any(cutoff < 0 | cutoff > 1) ) stop(\'SAINT score cutoff not between 0 and 1. Please correct and try again\')\r\n+}\r\n+\r\n+args <- commandArgs(trailingOnly = TRUE)\r\n+main(args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9])\r\n' |