annotate bubbles_v9_NSAF_natural_log.R @ 27:ab9ace68e3b5 draft

Uploaded

1 rm(list = ls()) 
2 ###################################################################################################
3 # R-code: Multi-bubble graph generation from SAINTexpress output
4 # Author: Brent Kuenzi
5 ###################################################################################################
6 # This Script generates the bubble graphs based upon Saint output. 
7 ###################################################################################################
8 # Copyright (C)  Brent Kuenzi.
9 # Permission is granted to copy, distribute and/or modify this document
10 # under the terms of the GNU Free Documentation License, Version 1.3
11 # or any later version published by the Free Software Foundation;
12 # with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
13 # A copy of the license is included in the section entitled "GNU
14 # Free Documentation License".
15 ###################################################################################################
16 ## REQUIRED INPUT ##
17 
18 # 1) listfile: SAINTexpress generated "list.txt" file
19 # 2) preyfile: SAINT pre-processing generated "prey.txt" file used to run SAINTexpress
20 
21 ## OPTIONAL INPUT ##
22 
23 # 3) crapome: raw output from crapome Workflow 1 query (http://www.crapome.org)
24 # 4) color: bubble color (default = "red")
25 #     - color = "crapome": color bubbles based on Crapome(%)
26 #     - Also recognizes any color within R's built-in colors() vector
27 # 5) label: Adds gene name labels to bubbles within the "zoomed in" graphs (default = FALSE)
28 # 6) cutoff: Saintscore cutoff to be assigned for filtering the "zoomed in" graphs (default = 0.8)
29 # 7) type: Specifies if the data is MaxQuant (MQ) or Scaffold (SC) data (default = "SC")
30 # 8) inc_file: Selects only the uniprot ids in the provided list (default ="None")
31 # 9) exc_file: Removes the proteins in the list (default = "None")
32 ###################################################################################################
33 
34 
35 ins_check_run <- function(){
36   if ('dplyr' %in% rownames(installed.packages())){}
37   else {
38     install.packages('dplyr', repos = 'http://cran.us.r-project.org')
39   }
40   if ('tidyr' %in% rownames(installed.packages())){}
41   else {
42     install.packages('tidyr', repos = 'http://cran.us.r-project.org')
43   }
44   if ('ggplot2' %in% rownames(installed.packages())){}
45   else {
46     install.packages('ggplot2', repos = 'http://cran.us.r-project.org')
47   }
48 }
49 
50 ins_check_run()
51 library(dplyr); library(tidyr); library(ggplot2)
52 
53 
54 main <- function(listfile, preyfile, crapome = FALSE, color = "red", label = FALSE, cutoff = 0.8, type = "SC", inc_file = "None", exc_file = "None" ) {
55   cutoff_check(cutoff)
56   listfile <- list_type(listfile, inc_file, exc_file)
57   if(type == "SC") {
58     df <- merge_files_sc(listfile, preyfile, crapome)
59   }
60   if(type == "MQ") {
61     df <- merge_files_mq(listfile, preyfile, crapome)
62   }
63   bubble_NSAF(df, color)
64   bubble_SAINT(df, color)
65   bubble_zoom_SAINT(df, color, label, cutoff)
66   bubble_zoom_NSAF(df, color, label, cutoff)
67   write.table(df, "output.txt", sep = "\t", quote = FALSE, row.names = FALSE)
68 }
69 ###################################################################################################
70 # Include and Exclude list filtering
71 ###################################################################################################
72 list_type <- function(df, inc_file, exc_file) {
73   Saint <- read.delim(df, stringsAsFactors = FALSE)
74   if (inc_file != "None") {
75     if (exc_file == "None") {
76       inc_prots <- read.delim(inc_file, sep = '\t', header = FALSE, stringsAsFactors = FALSE)
77       print(inc_prots[, 1])
78       print(Saint$Prey)
79       filtered_df = subset(Saint, Saint$Prey == inc_prots[, 1])
80     }
81     else {
82       inc_prots <- read.delim(inc_file, sep = '\t', header = FALSE, stringsAsFactors = FALSE)
83       exc_prots <- read.delim(exc_file, sep = '\t', header = FALSE, stringsAsFactors = FALSE)
84       filtered_df = subset(Saint, Saint$Prey == inc_prots[, 1])
85       filtered_df = subset(filtered_df, filtered_df$Prey != exc_prots[, 1])
86     }
87   }
88   else if (exc_file != "None") {
89     exc_prots <- read.delim(exc_file, sep = '\t', header = FALSE, stringsAsFactors = FALSE)
90     filtered_df = subset(Saint, Saint$Prey != exc_prots[, 1])
91   }
92   else {
93     filtered_df = Saint
94   }
95   return(filtered_df)
96   
97 }
98 ###################################################################################################
99 # Merge input files and caculate Crapome(%) and NSAF for each protein for each bait
00 ###################################################################################################
01 merge_files_mq <- function(SAINT, prey_DF, crapome = FALSE) {
02   #SAINT <- read.table(SAINT_DF, sep = '\t', header = TRUE)
03   #Some of these read.table()'s don't use stringsAsFactors = FALSE. Is this on purpose? Factors give rise to some really weird and unpredictable behavior; suggest always using stringsAsFactors = FALSE
04   prey <- read.table(prey_DF, sep = '\t', header = FALSE); colnames(prey) <- c("Prey", "Length", "PreyGene")
05   DF <- merge(SAINT, prey)
06   DF$SpecSum <- log2(DF$SpecSum)
07   
08   if(crapome != FALSE) {
09     crapome <- read.table(crapome, sep = '\t', header = TRUE)
10     colnames(crapome) <- c("Prey", "Symbol", "Num.of.Exp", "Ave.SC", "Max.SC")
11     DF1 <- merge(DF, crapome); as.character(DF1$Num.of.Exp); DF1$Symbol <- NULL;
12                     DF1$Ave.SC <- NULL; DF1$Max.SC <- NULL # Removes unnecessary columns.
13     DF1$Num.of.Exp <- sub("^$", "0 / 1", DF1$Num.of.Exp ) # Replace blank values with 0 / 1.
14     DF <- DF1 %>% separate(Num.of.Exp, c("NumExp", "TotalExp"), " / ") # Split into 2 columns.
15     DF$CrapomePCT <- 100 - (as.integer(DF$NumExp) / as.integer(DF$TotalExp) * 100) # Calculate the crapome %.
16   }
17   DF$SAF <- DF$AvgSpec / DF$Length
18   DF2 = DF %>% group_by(Bait) %>% mutate(NSAF = SAF/sum(SAF))
19   DF$NSAF = DF2$NSAF
20   return(DF)
21 }
22 
23 merge_files_sc <- function(SAINT, prey_DF, crapome = FALSE) {
24   #SAINT <- read.table(SAINT_DF, sep = '\t', header = TRUE)
25   prey <- read.table(prey_DF, sep = '\t', header = FALSE); colnames(prey) <- c("Prey", "Length", "PreyGene")
26   DF <- merge(SAINT, prey)
27   
28   if(crapome != FALSE) {
29     crapome <- read.table(crapome, sep = '\t', header = TRUE)
30     colnames(crapome) <- c("Prey", "Symbol", "Num.of.Exp", "Ave.SC", "Max.SC")
31     DF1 <- merge(DF, crapome); as.character(DF1$Num.of.Exp); DF1$Symbol <- NULL;
32                     DF1$Ave.SC <- NULL; DF1$Max.SC <- NULL # Removes unnecessary columns.
33     DF1$Num.of.Exp <- sub("^$", "0 / 1", DF1$Num.of.Exp ) # Replace blank values with 0 / 1.
34     DF <- DF1 %>% separate(Num.of.Exp, c("NumExp", "TotalExp"), " / ") # Split into 2 columns.
35     DF$CrapomePCT <- 100 - (as.integer(DF$NumExp) / as.integer(DF$TotalExp) * 100) # Calculate the crapome %.
36   }
37   DF$SAF <- DF$AvgSpec / DF$Length
38   DF2 = DF %>% group_by(Bait) %>% mutate(NSAF = SAF/sum(SAF))
39   DF$NSAF = DF2$NSAF
40   return(DF)
41 }
42 ###################################################################################################
43 # Plot all proteins for each bait by x = ln(NSAF), y = Log2(FoldChange)
44 ###################################################################################################
45 bubble_NSAF <- function(data, color) {
46     if(color == "crapome") {
47       a <- subset(data, CrapomePCT < 80, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait))
48       b <- subset(data, CrapomePCT >= 80, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait))
49       p <- qplot(x = log(NSAF), y = log2(FoldChange), data = a, colour = I("tan"), size = SpecSum) + scale_size(range = c(1, 10)) + 
50         geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)
51       if(length(levels(a$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")} # multiple graphs if multiple baits
52       #The text says ln() which is log base e, but the code uses log base 10. Fix code or the axis label. 
53       p <- p + geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum, color = CrapomePCT), data = b) + 
54         scale_colour_gradient(limits = c(80, 100), low = "tan", high = "red") + 
55         labs(colour = "CRAPome Probability \nof Specific Interaction (%)", x = "ln(NSAF)") + 
56         geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = b)
57       return(ggsave(p, width = 8, height = 4, filename = "bubble_NSAF.png"))
58     }
59    if(color != "crapome") {
60       p <- qplot(x = log(NSAF), y = log2(FoldChange), data = data, colour = I(color), size = SpecSum) + scale_size(range = c(1, 10)) + 
61         geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = data) + # add bubble outlines
62           labs(x = "ln(NSAF)")
63         if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
64       return(ggsave(p, width = 8, height = 4, filename = "bubble_NSAF.png"))
65     }
66   }
67 ###################################################################################################
68 # Plot all proteins for each bait by x = Saintscore, y = Log2(FoldChange)
69 ###################################################################################################
70 bubble_SAINT <- function(data, color) {
71     if(color == "crapome") {
72       a <- subset(data, CrapomePCT < 80, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait)) #filter on CRAPome
73       b <- subset(data, CrapomePCT >= 80, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait))
74       p <- qplot(x = SaintScore, y = log2(FoldChange), data = a, colour = I("tan"), size = SpecSum) + 
75         scale_size(range = c(1, 10)) + geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)
76       if(length(levels(a$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
77       p <- p + geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum, color = CrapomePCT), data = b) + 
78         scale_colour_gradient(limits = c(80, 100), low = "tan", high = "red") + 
79         labs(colour = "CRAPome Probability \nof Specific Interaction (%)") +
80         geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = b)
81       return(ggsave(p, width = 8, height = 4, filename = "bubble_SAINT.png"))
82     }
83     if(color != "crapome") {
84       p <- qplot(x = SaintScore, y = log2(FoldChange), data = data, colour = I(color), size = SpecSum) +
85         scale_size(range = c(1, 10)) + geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = data)
86       if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
87       return(ggsave(p, width = 8, height = 4, filename = "bubble_SAINT.png"))
88     }
89   }
90 ###################################################################################################
91 # Filter proteins on Saintscore cutoff and plot for each bait x = Saintscore, y = Log2(FoldChange)
92 ###################################################################################################
93 bubble_zoom_SAINT <- function(data, color, label = FALSE, cutoff = 0.8) {
94   if(color == "crapome") {
95     a <- subset(data, CrapomePCT < 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))
96     b <- subset(data, CrapomePCT >= 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))
97     p <- qplot(x = SaintScore, y = log2(FoldChange), data = a, colour = I("tan"), size = SpecSum) + 
98       scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score")+geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)
99     if(label == TRUE & length(a$NSAF != 0)) {
00       p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black")
01     }
02     if(length(levels(a$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
03     p <- p + geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum, color = CrapomePCT), data = b) + 
04       scale_colour_gradient(limits = c(80, 100), low = "tan", high = "red") + 
05       labs(colour = "CRAPome Probability \nof Specific Interaction (%)") + 
06       geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = b)
07     if(label == TRUE & length(b$NSAF != 0)) {
08       p <- p + geom_text(data = b, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)
09     }
10     return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_SAINT.png"))
11   }
12   if(color != "crapome") {
13     a <- subset(data, SaintScore >= cutoff, select = c(NSAF, SpecSum, FoldChange, SaintScore, Bait, PreyGene))
14     p <- qplot(x = SaintScore, y = log2(FoldChange), data = a, colour = I(color), size = SpecSum) +
15       scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + 
16       geom_point(aes(x = SaintScore, y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)
17     if(label == TRUE & length(a$NSAF != 0)) {
18       p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)
19     }
20     if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
21     return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_SAINT.png"))
22   }
23 }
24 ###################################################################################################
25 # Filter proteins on Saintscore cutoff and plot for each bait x = log(NSAF), y = Log2(FoldChange)
26 ###################################################################################################
27 bubble_zoom_NSAF <- function(data, color, label = FALSE, cutoff = 0.8) {
28   if(color == "crapome") {
29     a <- subset(data, CrapomePCT < 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))
30     b <- subset(data, CrapomePCT >= 80 & SaintScore >= cutoff, select = c(NSAF, SpecSum, CrapomePCT, FoldChange, SaintScore, Bait, PreyGene))
31     p <- qplot(x = log(NSAF), y = log2(FoldChange), data = a, colour = I("tan"), size = SpecSum) + 
32       scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + 
33       geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a)
34     if(label == TRUE & length(a$NSAF != 0)) {
35       p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black")
36     }
37     if(length(levels(a$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
38     p <- p + geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum, color = CrapomePCT), data = b) + 
39       scale_colour_gradient(limits = c(80, 100), low = "tan", high = "red") + 
40       labs(colour = "CRAPome Probability \nof Specific Interaction (%)", x = "ln(NSAF)") + 
41       geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = b)
42     if(label == TRUE & length(b$NSAF != 0)) {
43       p <- p + geom_text(data = b, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)
44     }
45     return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_NSAF.png"))
46   }
47   if(color != "crapome") {
48     a <- subset(data, SaintScore >= cutoff, select = c(NSAF, SpecSum, FoldChange, SaintScore, Bait, PreyGene))
49     p <- qplot(x = log(NSAF), y = log2(FoldChange), data = a, colour = I(color), size = SpecSum) +
50       scale_size(range = c(1, 10)) + ggtitle("Filtered on SAINT score") + 
51       geom_point(aes(x = log(NSAF), y = log2(FoldChange), size = SpecSum), colour = "black", shape = 21, data = a) + 
52       labs(x = "ln(NSAF)")
53     if(label == TRUE & length(a$NSAF != 0)) {
54       p <- p + geom_text(data = a, aes(label = PreyGene, size = 10, vjust = 0, hjust = 0), colour = "black", show_guide = FALSE)
55     }
56     if(length(levels(data$Bait) > 1)) {p <- p + facet_wrap(~Bait, scales = "free_y")}
57     return(ggsave(p, width = 8, height = 4, filename = "bubble_zoom_NSAF.png"))
58   }
59 }
60 ###################################################################################################
61 # Check Saintscore cutoff and stop program if not between 0 and 1
62 ###################################################################################################
63 cutoff_check <- function(cutoff){
64   if( any(cutoff < 0 | cutoff > 1) ) stop('SAINT score cutoff not between 0 and 1. Please correct and try again')
65 }
66 
67 args <- commandArgs(trailingOnly = TRUE)
68 main(args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9])