Finish R script linting

workflow/apps/r/CompileGTF.R

@@ -52,4 +52,3 @@ output_file <- file.path(outDir, "GTFCombined.xlsx")
 
 # Call the function to combine the files into a workbook with named sheets
 combineFilesToWorkbook(file_list, output_file)
-

workflow/apps/r/TSHeatmaps5dev2.R

@@ -36,11 +36,15 @@ base_dir <- args[6]
 output_dir <- args[7]
 study_nums <- args[8:length(args)]
 
-#import standard tables used in Sean's code That should be copied to each ExpStudy
+# Import standard tables used in Sean's code That should be copied to each ExpStudy
 labels <- read.csv(file = study_info_file, stringsAsFactors = FALSE)
 Ontology <- get_ontology(file = ontology_file, propagate_relationships = "is_a", extract_tags = "minimal")
 GO2ALLORFs <- as.list(org.Sc.sgdGO2ALLORFS) # all ORFs associated with GO term
-Terms <- read.delim(file=sgd_terms_tfile,header=FALSE,quote = "",col.names = c("GO_ID","GO_Term","GO_Aspect","GO_Term_Definition"))
+Terms <- read.delim(file = sgd_terms_tfile,
+  header = FALSE,
+  quote = "",
+  col.names = c("GO_ID", "GO_Term", "GO_Aspect", "GO_Term_Definition")
+)
 XX3 <- read.csv(file = all_sgd_terms_csv, stringsAsFactors = FALSE, header = TRUE)
 XX3[, 1] <- paste("GO:", formatC(XX3[, 1], width = 7, flag = "0"), sep = "")
 XX3[, 2] <- gsub(pattern = " ", replacement = "_", x = XX3[, 2])
@@ -48,7 +52,7 @@ XX3[,2] <- gsub(pattern = "/",replacement = "_",x = XX3[,2])
 
 # Load input files
 for (study_num in study_nums) {
-  input_file <- file.path(base_dir,paste('Exp',study_num),zscores_file)
+  input_file <- file.path(base_dir, paste("Exp", study_num), zscores_file)
   if (file.exists(input_file)) {
     assign(paste(X, study_num), read.csv(file = input_file, stringsAsFactors = FALSE, header = TRUE))
     assign(paste(Name, study_num), labels[study_num, 2])
@@ -76,12 +80,11 @@ if (length(study_nums) > 0) {
   try(X1[!is.na(X1$Z_lm_K) & X1$Z_lm_K >= 2, ]$Score_K <- "Deletion Suppressor")
   try(X1[!is.na(X1$Z_lm_K) & X1$Z_lm_K <= -2, ]$Score_K <- "Deletion Enhancer")
   
-  #express the na data as 0.001 in X1 for K and L
+  # Express the na data as 0.001 in X1 for K and L
   X1[is.na(X1$Z_lm_L), ]$Z_lm_L <- 0.001
   X1[is.na(X1$Z_lm_K), ]$Z_lm_K <- 0.001
   X1$Rank_L <- rank(X1$Z_lm_L)
   X1$Rank_K <- rank(X1$Z_lm_K)
-  
   X1 <- X1[order(X1$OrfRep, decreasing = FALSE), ]
   colnames(X1) <- paste(colnames(X1), "_X1", sep = "")
 }
@@ -130,13 +133,12 @@ if (length(study_nums) > 2) {
   try(X3[is.na(X3$Z_lm_K), ]$Score_K <- "No Growth")
   try(X3[!is.na(X3$Z_lm_K) & X3$Z_lm_K >= 2, ]$Score_K <- "Deletion Suppressor")
   try(X3[!is.na(X3$Z_lm_K) & X3$Z_lm_K <= -2, ]$Score_K <- "Deletion Enhancer")
-  #express the na data as 0.001 in X3
+  
+  # Express the na data as 0.001 in X3
   X3[is.na(X3$Z_lm_L), ]$Z_lm_L <- 0.001
   X3[is.na(X3$Z_lm_K), ]$Z_lm_K <- 0.001
-  
   X3$Rank_L <- rank(X3$Z_lm_L)
   X3$Rank_K <- rank(X3$Z_lm_K)
-  
   X3 <- X3[order(X3$OrfRep, decreasing = FALSE), ]
   colnames(X3) <- paste(colnames(X3), "_X3", sep = "")
   X <- cbind(X, X3)
@@ -158,13 +160,12 @@ if (length(study_nums) > 3) {
   try(X4[is.na(X4$Z_lm_K), ]$Score_K <- "No Growth")
   try(X4[!is.na(X4$Z_lm_K) & X4$Z_lm_K >= 2, ]$Score_K <- "Deletion Suppressor")
   try(X4[!is.na(X4$Z_lm_K) & X4$Z_lm_K <= -2, ]$Score_K <- "Deletion Enhancer")
-  #express the na data as 0.001 in X4
+  
+  # Express the na data as 0.001 in X4
   X4[is.na(X4$Z_lm_L), ]$Z_lm_L <- 0.001
   X4[is.na(X4$Z_lm_K), ]$Z_lm_K <- 0.001
-  
   X4$Rank_L <- rank(X4$Z_lm_L)
   X4$Rank_K <- rank(X4$Z_lm_K)
-  
   X4 <- X4[order(X4$OrfRep, decreasing = FALSE), ]
   colnames(X4) <- paste(colnames(X4), "_X4", sep = "")
   X <- cbind(X, X4)
@@ -186,13 +187,12 @@ if (length(study_nums) > 4) {
   try(X5[is.na(X5$Z_lm_K), ]$Score_K <- "No Growth")
   try(X5[!is.na(X5$Z_lm_K) & X5$Z_lm_K >= 2, ]$Score_K <- "Deletion Suppressor")
   try(X5[!is.na(X5$Z_lm_K) & X5$Z_lm_K <= -2, ]$Score_K <- "Deletion Enhancer")
-  #express the na data as 0.001 in X5
+  
+  # Express the na data as 0.001 in X5
   X5[is.na(X5$Z_lm_L), ]$Z_lm_L <- 0.001
   X5[is.na(X5$Z_lm_K), ]$Z_lm_K <- 0.001
-  
   X5$Rank_L <- rank(X5$Z_lm_L)
   X5$Rank_K <- rank(X5$Z_lm_K)
-  
   X5 <- X5[order(X5$OrfRep, decreasing = FALSE), ]
   colnames(X5) <- paste(colnames(X5), "_X5", sep = "")
   X <- cbind(X, X5)
@@ -201,12 +201,11 @@ if (length(study_nums) > 4) {
 X$ORF <- X$OrfRep_X1
 
 if (length(study_nums) > 1) {
-
   X$ORF <- gsub("_1", "", x = X$ORF)
-  
   try(X[X$Gene_X1 == "", ]$Gene_X1 <- X[X$Gene_X1 == "", ]$OrfRep_X1)
   try(X[X$Gene_X2 == "", ]$Gene_X2 <- X[X$Gene_X2 == "", ]$OrfRep_X2)
-  X_heatmap <- X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
+  X_heatmap <-
+    X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
     colnames(X) == "Z_Shift_K_X1" | colnames(X) == "Z_lm_K_X1" |
     colnames(X) == "Z_Shift_K_X2" | colnames(X) == "Z_lm_K_X2" |
     colnames(X) == "Z_Shift_L_X1" | colnames(X) == "Z_lm_L_X1" |
@@ -218,24 +217,24 @@ if (length(study_nums) > 1) {
   colnames(X_heatmap)[2] <- "Gene"
 }
 
-
 if (length(study_nums) > 2) {
   X$ORF <- gsub("_1", "", x = X$ORF)
   X$ORF <- gsub("_2", "", x = X$ORF)
-  
   try(X[X$Gene_X1 == "", ]$Gene_X1 <- X[X$Gene_X1 == "", ]$OrfRep_X1)
   try(X[X$Gene_X2 == "", ]$Gene_X2 <- X[X$Gene_X2 == "", ]$OrfRep_X2)
   try(X[X$Gene_X3 == "", ]$Gene_X3 <- X[X$Gene_X3 == "", ]$OrfRep_X3)
   
-  X_heatmap <- X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
+  X_heatmap <-
+    X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
     colnames(X) == "Z_Shift_K_X1" | colnames(X) == "Z_lm_K_X1" |
     colnames(X) == "Z_Shift_K_X2" | colnames(X) == "Z_lm_K_X2" |
     colnames(X) == "Z_Shift_K_X3" | colnames(X) == "Z_lm_K_X3" |
     colnames(X) == "Z_Shift_L_X1" | colnames(X) == "Z_lm_L_X1" |
     colnames(X) == "Z_Shift_L_X2" | colnames(X) == "Z_lm_L_X2" |
     colnames(X) == "Z_Shift_L_X3" | colnames(X) == "Z_lm_L_X3"]
+
   # Reorder columns
-  X_heatmap <- X_heatmap[,c(14,1,4,5,8,9,12,13,2,3,6,7,10,11)]  #Three
+  X_heatmap <- X_heatmap[, c(14, 1, 4, 5, 8, 9, 12, 13, 2, 3, 6, 7, 10, 11)]
   colnames(X_heatmap) <- gsub(pattern = "X1", replacement = Name1, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X2", replacement = Name2, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X3", replacement = Name3, colnames(X_heatmap))
@@ -246,13 +245,13 @@ if (length(study_nums) > 3) {
   X$ORF <- gsub("_1", "", x = X$ORF)
   X$ORF <- gsub("_2", "", x = X$ORF)
   X$ORF <- gsub("_3", "", x = X$ORF)
-  
   try(X[X$Gene_X1 == "", ]$Gene_X1 <- X[X$Gene_X1 == "", ]$OrfRep_X1)
   try(X[X$Gene_X2 == "", ]$Gene_X2 <- X[X$Gene_X2 == "", ]$OrfRep_X2)
   try(X[X$Gene_X3 == "", ]$Gene_X3 <- X[X$Gene_X3 == "", ]$OrfRep_X3)
   try(X[X$Gene_X4 == "", ]$Gene_X4 <- X[X$Gene_X4 == "", ]$OrfRep_X4)
   
-  X_heatmap <- X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
+  X_heatmap <-
+    X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
     colnames(X) == "Z_Shift_K_X1" | colnames(X) == "Z_lm_K_X1" |
     colnames(X) == "Z_Shift_K_X2" | colnames(X) == "Z_lm_K_X2" |
     colnames(X) == "Z_Shift_K_X3" | colnames(X) == "Z_lm_K_X3" |
@@ -261,8 +260,9 @@ if (length(study_nums) > 3) {
     colnames(X) == "Z_Shift_L_X2" | colnames(X) == "Z_lm_L_X2" |
     colnames(X) == "Z_Shift_L_X3" | colnames(X) == "Z_lm_L_X3" |
     colnames(X) == "Z_Shift_L_X4" | colnames(X) == "Z_lm_L_X4"]
+
   # Reorder columns
-  X_heatmap <- X_heatmap[,c(18,1,4,5,8,9,12,13,16,17,2,3,6,7,10,11,14,15)]  #Four
+  X_heatmap <- X_heatmap[, c(18, 1, 4, 5, 8, 9, 12, 13, 16, 17, 2, 3, 6, 7, 10, 11, 14, 15)]
   colnames(X_heatmap) <- gsub(pattern = "X1", replacement = Name1, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X2", replacement = Name2, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X3", replacement = Name3, colnames(X_heatmap))
@@ -270,20 +270,19 @@ if (length(study_nums) > 3) {
   colnames(X_heatmap)[2] <- "Gene"
 }
 
-
 if (length(study_nums) > 4) {
   X$ORF <- gsub("_1", "", x = X$ORF)
   X$ORF <- gsub("_2", "", x = X$ORF)
   X$ORF <- gsub("_3", "", x = X$ORF)
   X$ORF <- gsub("_4", "", x = X$ORF)
-  
   try(X[X$Gene_X1 == "", ]$Gene_X1 <- X[X$Gene_X1 == "", ]$OrfRep_X1)
   try(X[X$Gene_X2 == "", ]$Gene_X2 <- X[X$Gene_X2 == "", ]$OrfRep_X2)
   try(X[X$Gene_X3 == "", ]$Gene_X3 <- X[X$Gene_X3 == "", ]$OrfRep_X3)
   try(X[X$Gene_X4 == "", ]$Gene_X4 <- X[X$Gene_X4 == "", ]$OrfRep_X4)
   try(X[X$Gene_X5 == "", ]$Gene_X5 <- X[X$Gene_X5 == "", ]$OrfRep_X5)
   
-  X_heatmap <- X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
+  X_heatmap <-
+    X[colnames(X) == "ORF" | colnames(X) == "Gene_X1" |
     colnames(X) == "Z_Shift_K_X1" | colnames(X) == "Z_lm_K_X1" |
     colnames(X) == "Z_Shift_K_X2" | colnames(X) == "Z_lm_K_X2" |
     colnames(X) == "Z_Shift_K_X3" | colnames(X) == "Z_lm_K_X3" |
@@ -294,6 +293,7 @@ if (length(study_nums) > 4) {
     colnames(X) == "Z_Shift_L_X3" | colnames(X) == "Z_lm_L_X3" |
     colnames(X) == "Z_Shift_L_X4" | colnames(X) == "Z_lm_L_X4" |
     colnames(X) == "Z_Shift_L_X5" | colnames(X) == "Z_lm_L_X5"]
+
   # Reorder columns
   X_heatmap <- X_heatmap[, c(22, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 2, 3, 6, 7, 10, 11, 14, 15, 18, 19)]
   colnames(X_heatmap) <- gsub(pattern = "X1", replacement = Name1, colnames(X_heatmap))
@@ -301,16 +301,14 @@ if (length(study_nums) > 4) {
   colnames(X_heatmap) <- gsub(pattern = "X3", replacement = Name3, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X4", replacement = Name4, colnames(X_heatmap))
   colnames(X_heatmap) <- gsub(pattern = "X5", replacement = Name5, colnames(X_heatmap))
-  
   colnames(X_heatmap)[2] <- "Gene"
 }
 
-
-#theme elements for plots
+# Theme elements for plots
 theme_Publication <- function(base_size = 14, base_family = "sans") {
-  (theme_foundation(base_size=base_size, base_family=base_family)
-   + theme(plot.title = element_text(face = "bold",
-                                     size = rel(1.2), hjust = 0.5),
+  (theme_foundation(base_size = base_size, base_family = base_family) +
+    theme(
+      plot.title = element_text(face = "bold", size = rel(1.2), hjust = 0.5),
       text = element_text(),
       panel.background = element_rect(colour = NA),
       plot.background = element_rect(colour = NA),
@@ -332,26 +330,33 @@ theme_Publication <- function(base_size=14, base_family="sans") {
       plot.margin = unit(c(10, 5, 5, 5), "mm"),
       strip.background = element_rect(colour = "#f0f0f0", fill = "#f0f0f0"),
       strip.text = element_text(face = "bold")
-   ))
-  
+    )
+  )
 }
 
 scale_fill_Publication <- function(...) {
   library(scales)
-  discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
-  
+  discrete_scale(
+    "fill",
+    "Publication",
+    manual_pal(values = c("#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506", "#a6cee3", "#fb9a99", "#984ea3", "#ffff33")),
+    ...
+  )
 }
 
 scale_colour_Publication <- function(...) {
-  discrete_scale("colour","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
-  
+  discrete_scale(
+    "colour",
+    "Publication",
+    manual_pal(values = c("#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506", "#a6cee3", "#fb9a99", "#984ea3", "#ffff33")),
+    ...
+  )
 }
 
-
 theme_Publication_legend_right <- function(base_size = 14, base_family = "sans") {
-  (theme_foundation(base_size=base_size, base_family=base_family)
-   + theme(plot.title = element_text(face = "bold",
-                                     size = rel(1.2), hjust = 0.5),
+  (theme_foundation(base_size = base_size, base_family = base_family) +
+    theme(
+      plot.title = element_text(face = "bold", size = rel(1.2), hjust = 0.5),
       text = element_text(),
       panel.background = element_rect(colour = NA),
       plot.background = element_rect(colour = NA),
@@ -373,51 +378,68 @@ theme_Publication_legend_right <- function(base_size=14, base_family="sans") {
       plot.margin = unit(c(10, 5, 5, 5), "mm"),
       strip.background = element_rect(colour = "#f0f0f0", fill = "#f0f0f0"),
       strip.text = element_text(face = "bold")
-   ))
-  
+    )
+  )
 }
 
 scale_fill_Publication <- function(...) {
-  discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
-  
+  discrete_scale(
+    "fill",
+    "Publication",
+    manual_pal(values = c("#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506", "#a6cee3", "#fb9a99", "#984ea3", "#ffff33")),
+    ...
+  )
 }
 
 scale_colour_Publication <- function(...) {
-  discrete_scale("colour","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
-  
+  discrete_scale(
+    "colour",
+    "Publication",
+    manual_pal(values = c("#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506", "#a6cee3", "#fb9a99", "#984ea3", "#ffff33")),
+    ...
+  )
 }
 
-
 Ontology <- get_ontology(file = ontology_file, propagate_relationships = "is_a", extract_tags = "minimal")
 print(Ontology)
-#all ORFs associated with GO term
+
+# All ORFs associated with GO term
 GO2ALLORFs <- as.list(org.Sc.sgdGO2ALLORFS)
 
 # Terms is the GO term list jwr moved up to TAABLES
-Terms <- read.delim(file=sgd_terms_tfile,header=FALSE,quote = "",col.names = c("GO_ID","GO_Term","GO_Aspect","GO_Term_Definition"))
+Terms <- read.delim(file = sgd_terms_tfile,
+  header = FALSE,
+  quote = "",
+  col.names = c("GO_ID", "GO_Term", "GO_Aspect", "GO_Term_Definition")
+)
 
-
-#BIG LOOP BIG LOOP ------------------------------------------------------
 colormapbreaks <- c(-12, -10, -8, -6, -4, -2, 2, 4, 6, 8, 10, 12)
+
 for (s in 1:dim(XX3)[1]) {
-  #Ontology <- get_ontology(file="Documents/Hartman_Lab/SGD_Downloads/gene_ontology_edit.obo",propagate_relationships = "is_a",extract_tags = "minimal")
-  #Ontology_Everything <- get_ontology(file="Documents/Hartman_Lab/SGD_Downloads/gene_ontology_edit.obo",propagate_relationships = "is_a",extract_tags = "everything")
+  # Ontology <-
+  #   get_ontology(file = "Documents/Hartman_Lab/SGD_Downloads/gene_ontology_edit.obo",
+  #     propagate_relationships = "is_a", extract_tags = "minimal")
+  # Ontology_Everything <-
+  #   get_ontology(file = "Documents/Hartman_Lab/SGD_Downloads/gene_ontology_edit.obo",
+  #     propagate_relationships = "is_a", extract_tags = "everything")
+  #
   # GO_ID_Arg <- "GO:0006325"
   GO_ID_Arg_loop <- as.character(XX3[s, 1])
   GOTerm_parent <- get_descendants(Ontology, roots = GO_ID_Arg_loop)
   # GOTerm_parent <- get_descendants(Ontology,roots = "GO:0006325")
-  #only make plots if parent term has fewer than 500 children
- 
+  # Only make plots if parent term has fewer than 500 children
   Parent_Size <- length(as.vector(GO2ALLORFs[GO_ID_Arg_loop][[1]]))
   
   if (length(GOTerm_parent) > 100) {
     next()
   }
   
-  Parent_Size <- length(as.vector(GO2ALLORFs[GO_ID_Arg_loop][[1]]))
+  Parent_Size <- length(as.vector(GO2ALLORFs[GO_ID_Arg_loop][[1]])) # TODO why twice?
+
   if (Parent_Size < 2) {
     next()
   }
+
   if (Parent_Size > 2000) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 45, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -429,7 +451,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.5, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -439,7 +462,8 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
@@ -456,7 +480,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.6, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -466,11 +491,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size >= 500 && Parent_Size <= 1000) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 30, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -482,7 +509,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.6, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -492,11 +520,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size >= 200 && Parent_Size <= 500) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 25, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -508,7 +538,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -518,11 +549,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size >= 100 && Parent_Size <= 200) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 20, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -534,7 +567,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -544,11 +578,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size >= 60 && Parent_Size <= 100) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 15, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -560,7 +596,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -570,11 +607,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size >= 30 && Parent_Size <= 60) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 10, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -586,7 +625,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -596,10 +636,12 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           # ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
       if (dim(Genes_Annotated_to_Term)[1] <= 2 && dim(Genes_Annotated_to_Term)[1] > 0) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "none", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -609,12 +651,14 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           #ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     
     dev.off()
   }
+
   if (Parent_Size >= 3 && Parent_Size <= 30) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 7, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -626,7 +670,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -636,10 +681,12 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           # ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
       if (dim(Genes_Annotated_to_Term)[1] <= 2 && dim(Genes_Annotated_to_Term)[1] > 0) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "none", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -649,11 +696,13 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           # ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()
   }
+
   if (Parent_Size == 2) {
     pdf(file = paste(output_dir, XX3[s, 2], ".pdf", sep = ""), width = 12, height = 7, onefile = TRUE)
     for (i in 1:length(GOTerm_parent)) {
@@ -665,7 +714,8 @@ for(s in 1:dim(XX3)[1]){
       Genes_Annotated_to_Term <- X_heatmap[X_heatmap$ORF %in% All_Genes_Annotated_to_Term, ]
       X0 <- as.matrix(Genes_Annotated_to_Term[, 3:dim(Genes_Annotated_to_Term)[2]])
       if (dim(Genes_Annotated_to_Term)[1] > 2) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "row", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -675,10 +725,12 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           # ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
       if (dim(Genes_Annotated_to_Term)[1] <= 2 && dim(Genes_Annotated_to_Term)[1] > 0) {
-        try(heatmap.2(x=X0,
+        try(heatmap.2(
+          x = X0,
           Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
           dendrogram = "none", cexCol = 0.7, cexRow = 0.7, scale = "none",
           breaks = colormapbreaks, symbreaks = FALSE, colsep = c(2, 4, 6), sepcolor = "white", offsetCol = 0.1,
@@ -688,7 +740,8 @@ for(s in 1:dim(XX3)[1]){
           na.color = "red", col = brewer.pal(11, "PuOr"),
           main = GO_Term_Name,
           # ColSideColors = ev_repeat,
-                      labRow=as.character(Genes_Annotated_to_Term$Gene)))
+          labRow = as.character(Genes_Annotated_to_Term$Gene)
+        ))
       }
     }
     dev.off()

workflow/apps/r/gtaTemplate.R

@@ -40,16 +40,14 @@ if (length(args) >= 5) {
   output_dir <- "../../out/gta" # https://downloads.yeastgenome.org/curation/chromosomal_feature/gene_association.sgd
 }
 
-
-
 # # Set SGDgeneList file path
 # if (length(args) > 4) {
 #   SGDgeneList <- args[4]
 # } else {
 #   SGDgeneList <- "../Code/SGD_features.tab"
 
-
-#Begin for loop for experiments in this study-----------------ZScores_Interaction.csv
+# Begin for loop for experiments in this study
+# ZScores_Interaction.csv
 for (m in 1:length(zscores_file)) {
   
   #zscores_file <- paste(Wstudy,"/",expName[m],'/ZScores/ZScores_Interaction.csv',sep="")  #ArgsScore[1]
@@ -60,35 +58,41 @@ for(m in 1:length(zscores_file)){
   }
   
   #Terms is the GO term list
-  Terms <- read.delim(file = sgd_terms_file,header=FALSE,quote = "",col.names = c("GO_ID","GO_Term","GO_Aspect","GO_Term_Definition"))
+  Terms <- read.delim(file = sgd_terms_file, header = FALSE, quote = "",
+    col.names = c("GO_ID", "GO_Term", "GO_Aspect", "GO_Term_Definition"))
   #all ORFs associated with GO term
   GO2ALLORFs <- as.list(org.Sc.sgdGO2ALLORFS)
   #Gene_Association is the gene association to GO term file
-  Gene_Association <- read.delim(sgd_features_file,skip=8,header=FALSE,quote="",col.names = c("Database","Database_Object_ID","Database_Object_Symbol","NOT","GO_ID","Database_Reference","Evidence","With_or_From","Aspect","Database_Object_Name","Database_Object_Synonym","Database_Object_Type","taxon","Date","Assigned_By","OtherInfo","Empty"))
+  Gene_Association <- read.delim(sgd_features_file, skip = 8, header = FALSE, quote = "",
+    col.names = c("Database", "Database_Object_ID", "Database_Object_Symbol", "NOT", "GO_ID",
+      "Database_Reference", "Evidence", "With_or_From", "Aspect", "Database_Object_Name",
+      "Database_Object_Synonym", "Database_Object_Type", "taxon", "Date", "Assigned_By", "OtherInfo", "Empty"
+    )
+  )
   #Get the ORF names associated with each gene/GO term
   Gene_Association$ORF <- str_split_fixed(as.character(Gene_Association$Database_Object_Synonym), "\\|", 2)[, 1]
   #Get the numeric GO ID for matching
   Gene_Association$GO_ID_Numeric <- as.integer(str_split_fixed(as.character(Gene_Association$GO_ID), "\\:", 2)[, 2])
-  #get all unique GO terms
+  #Get all unique GO terms
   GO_Terms <- unique(Gene_Association$GO_ID)
-  #create a character vector with just the ColNames of the input file to store the scores for each GO term
+  #Create a character vector with just the ColNames of the input file to store the scores for each GO term
   Col_Names_X <- colnames(X)
-  #create a data_frame with header from input_file
+  #Create a data_frame with header from input_file
   GO_Term_Averages <- X[0, ]
-  #fill table with NAs same length as number of GO terms
+  #Fill table with NAs same length as number of GO terms
   GO_Term_Averages[1:length(GO_Terms), ] <- NA
-  #change the first and second col names to GO_ID and Term
+  #Change the first and second col names to GO_ID and Term
   colnames(GO_Term_Averages)[1] <- "GO_ID"
   colnames(GO_Term_Averages)[2] <- "Term"
   
-  #create new columns for Ontology, number genes (used to calculate the avg score), all possible genes in the GO term, and print genes/ORFs used
+  # Create new columns
   GO_Term_Averages$Ontology <- NA
   GO_Term_Averages$NumGenes <- NA
   GO_Term_Averages$AllPossibleGenes <- NA
   GO_Term_Averages$Genes <- NA
   GO_Term_Averages$ORFs <- NA
   
-  #create a data.frame for the standard deviation info 
+  # Create a data.frame for the standard deviation info
   GO_Term_SD <- X[0, ]
   GO_Term_SD[1:length(GO_Terms), ] <- NA
   colnames(GO_Term_SD)[1] <- "GO_ID"
@@ -96,72 +100,72 @@ for(m in 1:length(zscores_file)){
   
   # Loop for each GO term to get an average L and K Z score
   for (i in 1:length(GO_Terms)) {
-    #get the GO_Term
+    # Get the GO_Term
     ID <- GO_Terms[i]
     # Get data.frame for all genes associated to the GO Term
     ID_AllGenes <- Gene_Association[Gene_Association$GO_ID == ID, ]
-    #get a vector of just the gene names
+    # Get a vector of just the gene names
     ID_AllGenes_vector <- as.vector(GO2ALLORFs[as.character(ID)][[1]])
     if (length(unique(ID_AllGenes_vector)) > 4000) {
       next()
     }
-    #get the GO term character description where numeric Terms ID matches GO_Term's ID
+    # Get the GO term character description where numeric Terms ID matches GO_Term's ID
     GO_Description_Term <- as.character(Terms[Terms$GO_ID %in% ID_AllGenes$GO_ID_Numeric, ]$GO_Term[1])
-    #get the Z scores for all genes in the GO_ID
+    # Get the Z scores for all genes in the GO_ID
     Zscores_For_ID <- X[X$ORF %in% ID_AllGenes_vector, ]
-    #get the Gene names and ORFs for the term
+    # Get the Gene names and ORFs for the term
     GO_Term_Averages$Genes[i] <- paste(unique(Zscores_For_ID$Gene), collapse = " | ")
     GO_Term_Averages$ORFs[i] <- paste(unique(Zscores_For_ID$ORF), collapse = " | ")
     
-    #dataframe to report the averages for a GO term
-    #get the GO ID
+    # Dataframe to report the averages for a GO term
+    # Get the GO ID
     GO_Term_Averages$GO_ID[i] <- as.character(ID)
-    #get the term name
+    # Get the term name
     GO_Term_Averages$Term[i] <- GO_Description_Term
-    #get total number of genes annotated to the Term that we have in our library
+    # Get total number of genes annotated to the Term that we have in our library
     GO_Term_Averages$NumGenes[i] <- length(unique(Zscores_For_ID$ORF))
-    #get total number of genes annotated to the Term in SGD
+    # Get total number of genes annotated to the Term in SGD
     GO_Term_Averages$AllPossibleGenes[i] <- length(unique(ID_AllGenes_vector))
-    #get the ontology of the term
+    # Get the ontology of the term
     GO_Term_Averages$Ontology[i] <- as.character(ID_AllGenes$Aspect[1])
     
-    #calculate the average score for every column
+    # Calculate the average score for every column
     for (j in 3:length(X[1, ])) {
       GO_Term_Averages[i, j] <- mean(Zscores_For_ID[, j], na.rm = TRUE)
       # GO_Scores <- colMeans(Zscores_For_ID[,3:length(X[1,])])
     }
     
-    #also calculate same values for the SD
+    # Also calculate same values for the SD
     GO_Term_SD$GO_ID[i] <- as.character(ID)
-    #get the term name
+    # Get the term name
     GO_Term_SD$Term[i] <- GO_Description_Term
     
-    #calculate column scores for SD
+    # Calculate column scores for SD
     for (j in 3:length(X[1, ])) {
       GO_Term_SD[i, j] <- sd(Zscores_For_ID[, j], na.rm = TRUE)
       # GO_Scores <- colMeans(Zscores_For_ID[,3:length(X[1,])])
     }
   }
-  #add either _Avg or _SD depending on if the calculated score is an average or SD
+  # Add either _Avg or _SD depending on if the calculated score is an average or SD
   colnames(GO_Term_Averages) <- paste(colnames(GO_Term_Averages), "Avg", sep = "_")
   colnames(GO_Term_SD) <- paste(colnames(GO_Term_SD), "SD", sep = "_")
-  #combine the averages with the SDs to make one big data.frame
+  # Combine the averages with the SDs to make one big data.frame
   X2 <- cbind(GO_Term_Averages, GO_Term_SD)
-  #test[ , order(names(test))]
+  # Test[ , order(names(test))]
   X2 <- X2[, order(names(X2))]
   X2 <- X2[!is.na(X2$Z_lm_L_Avg), ]
-  #create output file
-  write.csv(X2,file=paste(output_dir,"/",expName[m],"/Average_GOTerms_All.csv",sep=""),row.names=FALSE)
-  #remove NAs
+  # Create output file
+  write.csv(X2, file.path(output_dir, expName[m], "Average_GOTerms_All.csv"), row.names = FALSE)
+  # Remove NAs
   X3 <- X2[!is.na(X2$Z_lm_L_Avg), ]
-  #identify redundant GO terms
   
+  # Identify redundant GO terms
   for (i in 1:length(X3[, 1])) {
-    #loop through each GO term - get term
+    # Loop through each GO term - get term
     GO_term_ID <- as.character(X3$GO_ID_Avg[i])
-    #get term in the X3
+    # Get term in the X3
     X3_Temp <- X3[X3$GO_ID_Avg == GO_term_ID, ]
-    #get anywhere that has the same number K_Avg value
+    # Get anywhere that has the same number K_Avg value
     X3_Temp2 <- X3[X3$Z_lm_K_Avg %in% X3_Temp, ]
     if (length(X3_Temp2[, 1]) > 1) {
       if (length(unique(X3_Temp2$Genes_Avg)) == 1) {
@@ -176,24 +180,17 @@ for(m in 1:length(zscores_file)){
       Y <- rbind(Y, X3_Temp2)
     }
   }
+
   Y1 <- unique(Y)
-  
-  write.csv(Y1,file=paste(output_dir,"/",exp_name,"/Average_GOTerms_All_NonRedundantTerms.csv",sep=""),row.names = FALSE)
-  
+  write.csv(Y1, file.path(output_dir, exp_name, "Average_GOTerms_All_NonRedundantTerms.csv"), row.names = FALSE)
   Y2 <- Y1[Y1$Z_lm_L_Avg >= 2 | Y1$Z_lm_L_Avg <= -2, ]
   Y2 <- Y2[!is.na(Y2$Z_lm_L_Avg), ]
-  write.csv(Y2,file=paste(output_dir,"/",exp_name,"/Average_GOTerms_NonRedundantTerms_Above2SD_L.csv",sep=""),row.names = FALSE)
-  
+  write.csv(Y2, file.path(output_dir, exp_name, "Average_GOTerms_NonRedundantTerms_Above2SD_L.csv"), row.names = FALSE)
   Y3 <- Y2[Y2$NumGenes_Avg > 2, ]
-  write.csv(Y3,file=paste(output_dir,"/",exp_name,"/Average_GOTerms_NonRedundantTerms_Above2SD_L_Above2Genes.csv",sep=""),row.names = FALSE)
-  
+  write.csv(Y3, file.path(output_dir, exp_name, "Average_GOTerms_NonRedundantTerms_Above2SD_L_Above2Genes.csv"), row.names = FALSE)
   Y4 <- Y1[Y1$Z_lm_K_Avg >= 2 | Y1$Z_lm_K_Avg <= -2, ]
   Y4 <- Y4[!is.na(Y4$Z_lm_K_Avg), ]
-  write.csv(Y4,file=paste(output_dir,"/",exp_name,"/Average_GOTerms_NonRedundantTerms_Above2SD_K.csv",sep=""),row.names = FALSE)
-  
+  write.csv(Y4, file.path(output_dir, exp_name, "Average_GOTerms_NonRedundantTerms_Above2SD_K.csv"), row.names = FALSE)
   Y5 <- Y4[Y4$NumGenes_Avg > 2, ]
-  write.csv(Y5,file=paste(output_dir,"/",exp_name,"/Average_GOTerms_NonRedundantTerms_Above2SD_K_Above2Genes.csv",sep=""),row.names = FALSE)
-  
-  #End of 'for loop'
+  write.csv(Y5, file.path(output_dir, exp_name, "Average_GOTerms_NonRedundantTerms_Above2SD_K_Above2Genes.csv"), row.names = FALSE)
 }
-

workflow/apps/r/joinInteractExps.R

@@ -9,7 +9,7 @@ args <- commandArgs(TRUE)
 
 # Set output dir
 if (length(args) >= 1) {
-  outDir <- args[1]
+  outDir <- file.path(args[1])
 } else {
   outDir <- "./" # for legacy workflow
 }
@@ -24,7 +24,7 @@ print(paste("SD=",sd))
 
 # Set studyInfo file
 if (length(args) >= 3) {
-  studyInfo <- args[3]
+  studyInfo <- file.path(args[3])
 } else {
   studyInfo <- "../Code/StudyInfo.csv" # for legacy workflow
 }
@@ -32,7 +32,7 @@ if (length(args) >= 3) {
 studies <- args[3:length(args)]
 inputFiles <- c()
 for (study in 1:length(studies)) {
-  zsFile <- file.path(study, 'zscores', 'zscores_interaction.csv')
+  zsFile <- file.path(study, "zscores", "zscores_interaction.csv")
   if (file.exists(zsFile)) {
     inputFiles[study] <- zsFile
   }
@@ -43,29 +43,31 @@ print(length(inputFiles))
 # TODO this is better handled in a loop in case you want to compare more experiments?
 # The input is already  designed for this
 # Read in the files for your experiment and
-# Join the two files at a time as a function of how many inputFile, list the larger file first ? in this example X2 has the larger number of genes.
+# Join the two files at a time as a function of how many inputFile
+# list the larger file first ? in this example X2 has the larger number of genes
 # If X1 has a larger number of genes, switch the order of X1 and X2
 if (length(inputFiles) == 2) {
   X1 <- read.csv(file = inputFiles[1], stringsAsFactors = FALSE)
   X2 <- read.csv(file = inputFiles[2], stringsAsFactors = FALSE)
   X <- join(X1, X2, by = "OrfRep")
-  OBH=X[,order(colnames(X))]  #OrderByHeader
-  headSel=select(OBH,  contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
-  headSel=select(headSel, -'Gene.1') #remove 'Gene.1 column
-  headSel2=select(OBH,  contains('OrfRep'), matches('Gene')) #Frame for interleaving Z_lm with Shift colums
-  headSel2=select(headSel2, -'Gene.1') #remove 'Gene.1 column   #Frame for interleaving Z_lm with Shift colums
-  
+  OBH <- X[, order(colnames(X))]  # OrderByHeader
+  headSel <- select(OBH, contains("OrfRep"), matches("Gene"),
+    contains("Z_lm_K"), contains("Z_Shift_K"), contains("Z_lm_L"), contains("Z_Shift_L"))
+  headSel <- select(headSel, -"Gene.1") #remove "Gene.1 column
+  headSel2 <- select(OBH,  contains("OrfRep"), matches("Gene")) #Frame for interleaving Z_lm with Shift colums
+  headSel2 <- select(headSel2, -"Gene.1") #remove "Gene.1 column   #Frame for interleaving Z_lm with Shift colums
 } else if (length(inputFiles) == 3) {
   X1 <- read.csv(file = inputFiles[1], stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
   X2 <- read.csv(file = inputFiles[2], stringsAsFactors = FALSE) #exp2File,stringsAsFactors = FALSE)
   X3 <- read.csv(file = inputFiles[3], stringsAsFactors = FALSE) #exp3File,stringsAsFactors = FALSE)
   X <- join(X1, X2, by = "OrfRep")
   X <- join(X, X3, by = "OrfRep")
-  OBH=X[,order(colnames(X))]  #OrderByHeader
-  headSel=select(OBH,  contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
-  headSel=select(headSel, -'Gene.1',-'Gene.2') 
-  headSel2=select(OBH,  contains('OrfRep'), matches('Gene')) 
-  headSel2=select(headSel2, -'Gene.1',-'Gene.2') 
+  OBH <- X[, order(colnames(X))]  #OrderByHeader
+  headSel <- select(OBH,  contains("OrfRep"), matches("Gene"),
+    contains("Z_lm_K"), contains("Z_Shift_K"), contains("Z_lm_L"), contains("Z_Shift_L"))
+  headSel <- select(headSel, -"Gene.1", -"Gene.2")
+  headSel2 <- select(OBH,  contains("OrfRep"), matches("Gene"))
+  headSel2 <- select(headSel2, -"Gene.1", -"Gene.2")
   
 } else if (length(inputFiles) == 4) {
   X1 <- read.csv(file = inputFiles[1], stringsAsFactors = FALSE) #exp1File,stringsAsFactors = FALSE)
@@ -75,46 +77,46 @@ if(length(inputFiles)==2) {
   X <- join(X1, X2, by = "OrfRep")
   X <- join(X, X3, by = "OrfRep")
   X <- join(X, X4, by = "OrfRep")
-  OBH=X[,order(colnames(X))]  #OrderByHeader
-  headSel=select(OBH,  contains('OrfRep'), matches('Gene'), contains('Z_lm_K'), contains('Z_Shift_K'),contains('Z_lm_L'), contains('Z_Shift_L'))
-  headSel=select(headSel, -'Gene.1',-'Gene.2',-'Gene.3') 
-  headSel2=select(OBH,  contains('OrfRep'), matches('Gene')) 
-  headSel2=select(headSel2, -'Gene.1',-'Gene.2',-'Gene.3') 
+  OBH <- X[, order(colnames(X))]  #OrderByHeader
+  headSel <- select(OBH,  contains("OrfRep"), matches("Gene"),
+    contains("Z_lm_K"), contains("Z_Shift_K"), contains("Z_lm_L"), contains("Z_Shift_L"))
+  headSel <- select(headSel, -"Gene.1", -"Gene.2", -"Gene.3")
+  headSel2 <- select(OBH,  contains("OrfRep"), matches("Gene"))
+  headSel2 <- select(headSel2, -"Gene.1", -"Gene.2", -"Gene.3")
 }
 
-#headSel$contains('Z_Shift') %>% replace_na(0.001)
+# headSel$contains("Z_Shift") %>% replace_na(0.001)
 headers <- colnames(headSel)
-i=0
+i <- 0
 for (i in 1:length(headers)) {
-  
   if (grepl("Shift", headers[i])) {
-    headSel[headers[i]][is.na(headSel[headers[i]])]=0.001
+    headSel[headers[i]][is.na(headSel[headers[i]])] <- 0.001
   }
   if (grepl("Z_lm_", headers[i])) {
-    headSel[headers[i]][is.na(headSel[headers[i]])]=0.0001
+    headSel[headers[i]][is.na(headSel[headers[i]])] <- 0.0001
   }
 }
 
 # 2SD option code to exclude Z_lm values less than 2 standard Deviations
-
-REMcRdy=select(headSel,  contains('OrfRep'), matches('Gene'), contains('Z_lm_'))
-shiftOnly=select(headSel,  contains('OrfRep'), matches('Gene'), contains('Z_Shift'))
+REMcRdy <- select(headSel,  contains("OrfRep"), matches("Gene"), contains("Z_lm_"))
+shiftOnly <- select(headSel,  contains("OrfRep"), matches("Gene"), contains("Z_Shift"))
 
 # Code to replace the numeric (.1 .2 .3) headers with experiment names from StudyInfo.txt
 Labels <- read.csv(file = "../Code/StudyInfo.csv", stringsAsFactors = FALSE, sep = ",")
 
 # Using Text search grepl to relabel headers
-REMcRdyHdr=colnames(REMcRdy)
-REMcRdyLabels='asdf'
-shftHdr=colnames(shiftOnly)
-shiftLabels='asdf'
+REMcRdyHdr <- colnames(REMcRdy)
+REMcRdyLabels <- "asdf"
+shftHdr <- colnames(shiftOnly)
+shiftLabels <- "asdf"
 shiftLabels[1:2] <- shftHdr[1:2]
 REMcRdyLabels[1:2] <- REMcRdyHdr[1:2]
   
 for (i in 3:(length(shftHdr))) {
   if (i == 3) {
     shiftLabels[3] <- paste0(Labels[1, 2], ".", shftHdr[3])
-    REMcRdyLabels[3]<-paste0(Labels[1,2],".",REMcRdyHdr[3]) }
+    REMcRdyLabels[3] <- paste0(Labels[1, 2], ".", REMcRdyHdr[3])
+  }
   if (i == 5) {
     shiftLabels[5] <- paste0(Labels[1, 2], ".", shftHdr[5])
     REMcRdyLabels[5] <- paste0(Labels[1, 2], ".", REMcRdyHdr[5])
@@ -125,17 +127,20 @@ for(i in 3:(length(shftHdr))){
   }
   if (grepl(".1", shftHdr[i], fixed = true)) {
     shiftLabels[i] <- paste0(Labels[2, 2], ".", shftHdr[i])
-    REMcRdyLabels[i]<-paste0(Labels[2,2],".",REMcRdyHdr[i])}
+    REMcRdyLabels[i] <- paste0(Labels[2, 2], ".", REMcRdyHdr[i])
+  }
   if (grepl(".2", shftHdr[i], fixed = true)) {
     shiftLabels[i] < -paste0(Labels[3, 2], ".", shftHdr[i])
-    REMcRdyLabels[i]<-paste0(Labels[3,2],".",REMcRdyHdr[i])}
+    REMcRdyLabels[i] <- paste0(Labels[3, 2], ".", REMcRdyHdr[i])
+  }
   if (grepl(".3", shftHdr[i], fixed = true)) {
     shiftLabels[i] <- paste0(Labels[4, 2], ".", shftHdr[i])
-    REMcRdyLabels[i]<-paste0(Labels[4,2],".",REMcRdyHdr[i])}
+    REMcRdyLabels[i] <- paste0(Labels[4, 2], ".", REMcRdyHdr[i])
+  }
 }
 
 for (i in 3:(length(REMcRdyLabels))) {
-  j=as.integer(i)
+  j <- as.integer(i)
   REMcRdyLabels[j] <- gsub("[.]", "_", REMcRdyLabels[j])
   shiftLabels[j] <- gsub("[.]", "_", shiftLabels[j])
 }
@@ -143,86 +148,86 @@ for(i in 3:(length(REMcRdyLabels))){
 colnames(shiftOnly) <- shiftLabels
 colnames(REMcRdy) <- REMcRdyLabels
 
-combI=headSel2 #Starting Template orf, Genename columns
+combI <- headSel2 # starting Template orf, Genename columns
 
 # headersRemc<-colnames(REMcRdy)
 # Reorder columns to produce an interleaved set of Z_lm and Shift data for all the cpps.
 for (i in 3:length(colnames(REMcRdy))) {
-  combI=cbind.data.frame(combI, shiftOnly[i])
-  combI=cbind.data.frame(combI, REMcRdy[i])
+  combI <- cbind.data.frame(combI, shiftOnly[i])
+  combI <- cbind.data.frame(combI, REMcRdy[i])
 }
 
-Vec1=NA
-Vec2=NA
-Vec3=NA
-Vec4=NA
-Vec5=NA
-Vec6=NA
-Vec7=NA
-Vec8=NA
+Vec1 <- NA
+Vec2 <- NA
+Vec3 <- NA
+Vec4 <- NA
+Vec5 <- NA
+Vec6 <- NA
+Vec7 <- NA
+Vec8 <- NA
 
 if (length(REMcRdy) == 6) {
-  Vec1=abs(REMcRdy[,3])>=std
-  Vec2=abs(REMcRdy[,4])>=std
-  Vec3=abs(REMcRdy[,5])>=std
-  Vec4=abs(REMcRdy[,6])>=std
-  bolVec=Vec1 | Vec2 |Vec3 |Vec4
-  REMcRdyGT2=REMcRdy[bolVec,1:2]
-  REMcRdyGT2[ ,3:6]=REMcRdy[bolVec,3:6]
-  shiftOnlyGT2=shiftOnly[bolVec,1:2]
-  shiftOnlyGT2[ ,3:6]=shiftOnly[bolVec,3:6]
+  Vec1 <- abs(REMcRdy[, 3]) >= std
+  Vec2 <- abs(REMcRdy[, 4]) >= std
+  Vec3 <- abs(REMcRdy[, 5]) >= std
+  Vec4 <- abs(REMcRdy[, 6]) >= std
+  bolVec <- Vec1 | Vec2 | Vec3 | Vec4
+  REMcRdyGT2 <- REMcRdy[bolVec, 1:2]
+  REMcRdyGT2[, 3:6] <- REMcRdy[bolVec, 3:6]
+  shiftOnlyGT2 <- shiftOnly[bolVec, 1:2]
+  shiftOnlyGT2[, 3:6] <- shiftOnly[bolVec, 3:6]
 }
 
 if (length(REMcRdy) == 8) {
-  Vec1=abs(REMcRdy[,3])>=std
-  Vec2=abs(REMcRdy[,4])>=std
-  Vec3=abs(REMcRdy[,5])>=std
-  Vec4=abs(REMcRdy[,6])>=std
-  Vec5=abs(REMcRdy[,7])>=std
-  Vec6=abs(REMcRdy[,8])>=std
-  bolVec=Vec1 | Vec2 |Vec3 | Vec4 |Vec5 |Vec6
-  REMcRdyGT2=REMcRdy[bolVec,1:2]
-  REMcRdyGT2[ ,3:8]=REMcRdy[bolVec,3:8]
-  shiftOnlyGT2=shiftOnly[bolVec,1:2]
-  shiftOnlyGT2[ ,3:8]=shiftOnly[bolVec,3:8]
+  Vec1 <- abs(REMcRdy[, 3]) >= std
+  Vec2 <- abs(REMcRdy[, 4]) >= std
+  Vec3 <- abs(REMcRdy[, 5]) >= std
+  Vec4 <- abs(REMcRdy[, 6]) >= std
+  Vec5 <- abs(REMcRdy[, 7]) >= std
+  Vec6 <- abs(REMcRdy[, 8]) >= std
+  bolVec <- Vec1 | Vec2 | Vec3 | Vec4 | Vec5 | Vec6
+  REMcRdyGT2 <- REMcRdy[bolVec, 1:2]
+  REMcRdyGT2[, 3:8] <- REMcRdy[bolVec, 3:8]
+  shiftOnlyGT2 <- shiftOnly[bolVec, 1:2]
+  shiftOnlyGT2[, 3:8] <- shiftOnly[bolVec, 3:8]
 }
 
 if (length(REMcRdy) == 10) {
-  Vec1=abs(REMcRdy[,3])>=std
-  Vec2=abs(REMcRdy[,4])>=std
-  Vec3=abs(REMcRdy[,5])>=std
-  Vec4=abs(REMcRdy[,6])>=std
-  Vec5=abs(REMcRdy[,7])>=std
-  Vec6=abs(REMcRdy[,8])>=std
-  Vec7=abs(REMcRdy[,9])>=std
-  Vec8=abs(REMcRdy[,10])>=std
-  bolVec=Vec1 | Vec2 |Vec3 |Vec4|Vec5|Vec6|Vec7|Vec8 
-  REMcRdyGT2=REMcRdy[bolVec,1:2]
-  REMcRdyGT2[ ,3:10]=REMcRdy[bolVec,3:10]
-  shiftOnlyGT2=shiftOnly[bolVec,1:2]
-  shiftOnlyGT2[ ,3:10]=shiftOnly[bolVec,3:10]
+  Vec1 <- abs(REMcRdy[, 3]) >= std
+  Vec2 <- abs(REMcRdy[, 4]) >= std
+  Vec3 <- abs(REMcRdy[, 5]) >= std
+  Vec4 <- abs(REMcRdy[, 6]) >= std
+  Vec5 <- abs(REMcRdy[, 7]) >= std
+  Vec6 <- abs(REMcRdy[, 8]) >= std
+  Vec7 <- abs(REMcRdy[, 9]) >= std
+  Vec8 <- abs(REMcRdy[, 10]) >= std
+  bolVec <- Vec1 | Vec2 | Vec3 | Vec4 | Vec5 | Vec6 | Vec7 | Vec8
+  REMcRdyGT2 <- REMcRdy[bolVec, 1:2]
+  REMcRdyGT2[, 3:10] <- REMcRdy[bolVec, 3:10]
+  shiftOnlyGT2 <- shiftOnly[bolVec, 1:2]
+  shiftOnlyGT2[, 3:10] <- shiftOnly[bolVec, 3:10]
 }
 
 if (std != 0) {
-  REMcRdy=REMcRdyGT2  # [,2:length(REMcRdyGT2)]
-  shiftOnly=shiftOnlyGT2  # [,2:length(shiftOnlyGT2)]
+  REMcRdy <- REMcRdyGT2  # [,2:length(REMcRdyGT2)]
+  shiftOnly <- shiftOnlyGT2  # [,2:length(shiftOnlyGT2)]
 }
 
 if (std == 0) {
-  REMcRdy=REMcRdy # [,2:length(REMcRdy)]
-  shiftOnly=shiftOnly # [,2:length(shiftOnly)]
+  REMcRdy <- REMcRdy # [,2:length(REMcRdy)]
+  shiftOnly <- shiftOnly # [,2:length(shiftOnly)]
 }
 
 # R places hidden "" around the header names. The following
 # is intended to remove those quote so that the "" do not blow up the Java REMc.
 # Use ,quote=F in the write.csv statement to fix R output file.
-#write.csv(combI,file=file.path(outDir,"CombinedKLzscores.csv"),row.names=FALSE)
-write.csv(REMcRdy,file=file.path(outDir,"REMcRdy_lm_only.csv"),row.names=FALSE, quote=F)
-write.csv(shiftOnly,file=file.path(outDir,"Shift_only.csv"),row.names=FALSE, quote=F)
+# write.csv(combI,file.path(outDir,"CombinedKLzscores.csv"), row.names = FALSE)
+write.csv(REMcRdy, file.path(outDir, "REMcRdy_lm_only.csv"), row.names = FALSE, quote = FALSE)
+write.csv(shiftOnly, file.path(outDir, "Shift_only.csv"), row.names = FALSE, quote = FALSE)
 #LabelStd <- read.table(file="./parameters.csv",stringsAsFactors = FALSE,sep = ",")
 
 LabelStd <- read.csv(file = studyInfo, stringsAsFactors = FALSE)
 print(std)
-LabelStd[,4]=as.numeric(std)
+LabelStd[, 4] <- as.numeric(std)
 write.csv(LabelStd, file = file.path(outDir, "parameters.csv"), row.names = FALSE)
 write.csv(LabelStd, file = studyInfo, row.names = FALSE)

workflow/qhtcp-workflow

@@ -1620,7 +1620,12 @@ gta() {
   for combo in "${study_combos[@]}"; do
     # Split on comma and assign to array
     IFS=',' read -ra studies <<< "$combo"
-    r_gta_pairwiselk "${studies[0]}" "${studies[1]}" "$STUDY_INFO_FILE" "$gta_out_dir"
+    r_gta_pairwiselk \
+      "${studies[0]}" \
+      "${studies[1]}" \
+      "$STUDY_INFO_FILE" \
+      "Average_GOTerms_All.csv" \
+      "$gta_out_dir"
   done
 
   # All studies
@@ -1663,14 +1668,15 @@ wrapper r_gta
 # * Is GTAtemplate.R actually a template?
 # * Do we need to allow user customization?
 #
-# Files
+# INPUT
 #
 # * [gene_association.sgd](https://downloads.yeastgenome.org/curation/chromosomal_feature/gene_association.sgd)
 # * go_terms.tab
 #
-# Output
+# OUTPUT
+#
+# * Average_GOTerms_All.csv
 #
-#   *
 #
 # @arg $1 string Exp# name
 # @arg $2 string ZScores_Interaction.csv file
@@ -1683,7 +1689,7 @@ r_gta() {
   cat <<-EOF
 
 	EOF
-  script="$APPS_DIR/r/GTAtemplate.R"
+  script="$APPS_DIR/r/gtaTemplate.R"
   [[ -d $5 ]] || mkdir -p "$5"
   debug "$RSCRIPT $script $*"
   "$RSCRIPT" "$script" \
@@ -1701,11 +1707,13 @@ wrapper r_gta_pairwiselk
 #
 # TODO
 #
-#   * Should move directory creation from PairwiseLK.R to gta module
+# * Move directory creation from PairwiseLK.R to gta module
+# * Needs better output filenames and directory organization
+# * Needs more for looping to reduce verbosity
 #
-# Files
+# INPUT
 #
-#   * 
+# * Average_GOTerms_All.csv
 # * 
 #
 # Output
@@ -1724,7 +1732,7 @@ wrapper r_gta_pairwiselk
 # @arg $4 string output directory
 #
 r_gta_pairwiselk() {
-  debug "Running: ${FUNCNAME[0]}"
+  debug "Running: ${FUNCNAME[0]} $*"
   cat <<-EOF
 
 	EOF
@@ -1734,7 +1742,12 @@ r_gta_pairwiselk() {
   [[ -d $4 ]] || mkdir -p "$4"
 
   debug "$RSCRIPT $script $*"
-  "$RSCRIPT" "$script" "$@"
+  "$RSCRIPT" "$script" \
+    "$1" \
+    "$2" \
+    "$3" \
+    "$4" \
+    "${@:5}"
 }
 
 
@@ -1842,6 +1855,10 @@ r_interactions() {
 wrapper r_join_interactions
 # @description JoinInteractExps3dev.R creates REMcRdy_lm_only.csv and Shift_only.csv
 #
+# TODO
+# * Needs more loops to reduce verbosity
+#
+#
 # Output
 #
 # * REMcRdy_lm_only.csv