Lint R scripts

This commit is contained in:
2024-08-13 15:59:18 -04:00
parent f190967383
commit 79862ddab4
4 changed files with 185 additions and 170 deletions

1
.gitignore vendored
View File

@@ -3,3 +3,4 @@ manual.odt
mwe mwe
centos-upgrade-plan.txt centos-upgrade-plan.txt
workflow/out/ workflow/out/
workflow/scans/

View File

@@ -10,25 +10,25 @@ library(sos)
args <- commandArgs(TRUE) args <- commandArgs(TRUE)
if (length(args) >= 1) { if (length(args) >= 1) {
finalTable <- args[1] finalTable <- file.path(args[1])
} else { } else {
finalTable <- "REMcRdy_lm_only.csv-finalTable.csv" # for legacy workflow finalTable <- "REMcRdy_lm_only.csv-finalTable.csv" # for legacy workflow
} }
if (length(args) >= 2) { if (length(args) >= 2) {
shiftFile <- args[2] shiftFile <- file.path(args[2])
} else { } else {
shiftFile <- "Shift_only.csv" # for legacy workflow shiftFile <- "Shift_only.csv" # for legacy workflow
} }
if (length(args) >= 3) { if (length(args) >= 3) {
studyInfo <- args[3] studyInfo <- file.path(args[3])
} else { } else {
studyInfo <- "../Code/StudyInfo.csv" # for legacy workflow studyInfo <- "../Code/StudyInfo.csv" # for legacy workflow
} }
if (length(args) >= 4) { if (length(args) >= 4) {
output <- args[4] output <- file.path(args[4])
} else { } else {
output <- "REMcHeatmaps/REMcWithShift.csv" # for legacy workflow output <- "REMcHeatmaps/REMcWithShift.csv" # for legacy workflow
} }

View File

@@ -171,7 +171,7 @@ for (i in 1:num_unique_clusts) {
if (cluster_length != 1) { if (cluster_length != 1) {
X0 <- as.matrix(cluster_data[, 4:(length(hmapfile[1, ]) - 2)]) X0 <- as.matrix(cluster_data[, 4:(length(hmapfile[1, ]) - 2)])
if (cluster_length >= 2001) { if (cluster_length >= 2001) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 20, width = 15) pdf(file = mypath, height = 20, width = 15)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -191,7 +191,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length >= 201 && cluster_length <= 2000) { if (cluster_length >= 201 && cluster_length <= 2000) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 15, width = 12) pdf(file = mypath, height = 15, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -210,7 +210,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length >= 150 && cluster_length <= 200) { if (cluster_length >= 150 && cluster_length <= 200) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 12, width = 12) pdf(file = mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -228,7 +228,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length >= 101 && cluster_length <= 149) { if (cluster_length >= 101 && cluster_length <= 149) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, mypath, height = 12, width = 12) pdf(file = mypath, mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -246,7 +246,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length >= 60 && cluster_length <= 100) { if (cluster_length >= 60 && cluster_length <= 100) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 12, width = 12) pdf(file = mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -264,7 +264,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length <= 59 && cluster_length >= 30) { if (cluster_length <= 59 && cluster_length >= 30) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 9, width = 12) pdf(file = mypath, height = 9, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,
@@ -282,7 +282,7 @@ for (i in 1:num_unique_clusts) {
dev.off() dev.off()
} }
if (cluster_length <= 29) { if (cluster_length <= 29) {
mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), sep = ""), ".pdf") mypath <- file.path(outDir, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file = mypath, height = 7, width = 12) pdf(file = mypath, height = 7, width = 12)
heatmap.2( heatmap.2(
x = X0, x = X0,

View File

@@ -8,57 +8,58 @@ library(tidyverse)
args <- commandArgs(TRUE) args <- commandArgs(TRUE)
# Need to give the input "finalTable.csv" file after running REMc generated by eclipse # Need to give the input "finalTable.csv" file after running REMc generated by eclipse
inputFinalTable <- args[1] inputFinalTable <- file.path(args[1])
# Give the DAmP_list.txt as the third argument - will color the gene names differently # Give the DAmP_list.txt as the third argument - will color the gene names differently
DAmPs <- Args[2] DAmPs <- file.path(Args[2])
DAmP_list <- read.delim(file=DAmPs,header=F,stringsAsFactors = F) DAmP_list <- read.delim(file = DAmPs, header = FALSE, stringsAsFactors = FALSE)
# Give the yeast human homology mapping as the fourth argument - will add the genes to the finalTable and use info for heatmaps # Give the yeast human homology mapping as the fourth argument - will add the genes to the finalTable and use info for heatmaps
mapFile <- Args[3] mapFile <- file.path(Args[3])
mapping <- read.csv(file=mapFile,stringsAsFactors = F) mapping <- read.csv(file = mapFile, stringsAsFactors = FALSE)
# Define the output path for the heatmaps - create this folder first - in linux terminal in the working folder use > mkdir filename_heatmaps # Define the output path for the heatmaps - create this folder first - in linux terminal in the working folder use > mkdir filename_heatmaps
outputPath <- Args[4] outputPath <- file.path(Args[4])
# Read in finalTablewithShift # Read in finalTablewithShift
hmapfile <- data.frame(read.csv(file=inputFinalTable,header=TRUE,sep=",",stringsAsFactors = FALSE)) hmapfile <- data.frame(read.csv(file = inputFinalTable, header = TRUE, sep = ",", stringsAsFactors = FALSE))
# Map the finalTable to the human homolog file # Map the finalTable to the human homolog file
hmapfile_map <- hmapfile hmapfile_map <- hmapfile
# Match using OrfRep after dropping the _1 _2 _3 _4 # Match using OrfRep after dropping the _1 _2 _3 _4
# But need to also account for some older files have ORF as column name rather than OrfRep in finalTable file # But need to also account for some older files have ORF as column name rather than OrfRep in finalTable file
if(colnames(hmapfile_map)[2] == "OrfRep"){ if (colnames(hmapfile_map)[2] == "OrfRep") {
try(hmapfile_map$ORFMatch <- hmapfile_map$OrfRep) try(hmapfile_map$ORFMatch <- hmapfile_map$OrfRep)
} }
if(colnames(hmapfile_map)[2] == "ORF"){ if (colnames(hmapfile_map)[2] == "ORF") {
try(hmapfile_map$ORFMatch <- hmapfile_map$ORF) try(hmapfile_map$ORFMatch <- hmapfile_map$ORF)
} }
hmapfile_map$ORFMatch <- gsub("_1","",x=hmapfile_map$ORFMatch) hmapfile_map$ORFMatch <- gsub("_1", "", x = hmapfile_map$ORFMatch)
hmapfile_map$ORFMatch <- gsub("_2","",x=hmapfile_map$ORFMatch) hmapfile_map$ORFMatch <- gsub("_2", "", x = hmapfile_map$ORFMatch)
hmapfile_map$ORFMatch <- gsub("_3","",x=hmapfile_map$ORFMatch) hmapfile_map$ORFMatch <- gsub("_3", "", x = hmapfile_map$ORFMatch)
hmapfile_map$ORFMatch <- gsub("_4","",x=hmapfile_map$ORFMatch) hmapfile_map$ORFMatch <- gsub("_4", "", x = hmapfile_map$ORFMatch)
# Join the hmapfile using # Join the hmapfile using
hmapfile_w_homolog <- full_join(hmapfile_map,mapping,by=c("ORFMatch"="ensembl_gene_id")) hmapfile_w_homolog <- full_join(hmapfile_map, mapping, by = c("ORFMatch" = "ensembl_gene_id"))
# Remove matches that are not from the finalTable # Remove matches that are not from the finalTable
hmapfile_w_homolog <- hmapfile_w_homolog[is.na(hmapfile_w_homolog$likelihood) == F,] hmapfile_w_homolog <- hmapfile_w_homolog[is.na(hmapfile_w_homolog$likelihood) == FASLE, ]
# Write csv with all info from mapping file # Write csv with all info from mapping file
write.csv(hmapfile_w_homolog,file=paste(outputPath,"/",inputFinalTable,"_WithHomologAll.csv",sep=""),row.names = F) write.csv(hmapfile_w_homolog, file.path(outputPath, paste(inputFinalTable, "_WithHomologAll.csv", sep = "")), row.names = FALSE)
# Remove the non matches and output another mapping file - this is also one used to make heatmaps # Remove the non matches and output another mapping file - this is also one used to make heatmaps
hmapfile_w_homolog <- hmapfile_w_homolog[is.na(hmapfile_w_homolog$external_gene_name_Human) == F,] hmapfile_w_homolog <- hmapfile_w_homolog[is.na(hmapfile_w_homolog$external_gene_name_Human) == FALSE, ]
write.csv(hmapfile_w_homolog,file=paste(outputPath,"/",inputFinalTable,"_WithHomologMatchesOnly.csv",sep=""),row.names = F) write.csv(hmapfile_w_homolog, file.path(outputPath, paste(inputFinalTable, "_WithHomologMatchesOnly.csv", sep = ""), row.names = FALSE))
# Add human gene name to the Gene column # Add human gene name to the Gene column
hmapfile_w_homolog$Gene <- paste(hmapfile_w_homolog$Gene,hmapfile_w_homolog$external_gene_name_Human,sep="/") hmapfile_w_homolog$Gene <- paste(hmapfile_w_homolog$Gene, hmapfile_w_homolog$external_gene_name_Human, sep = "/")
# Only keep the finalTable file columns and the homology info # Only keep the finalTable file columns and the homology info
hmap_len <- dim(hmapfile)[2] hmap_len <- dim(hmapfile)[2]
hmapfile_w_homolog_remake <- cbind(hmapfile_w_homolog[,1:hmap_len], hsapiens_homolog_orthology_type=hmapfile_w_homolog$hsapiens_homolog_orthology_type) hmapfile_w_homolog_remake <-
cbind(hmapfile_w_homolog[, 1:hmap_len], hsapiens_homolog_orthology_type = hmapfile_w_homolog$hsapiens_homolog_orthology_type)
hmapfile <- hmapfile_w_homolog_remake hmapfile <- hmapfile_w_homolog_remake
# Set NAs to NA # Set NAs to NA
@@ -68,19 +69,20 @@ hmapfile[hmapfile == 0.001] <- NA
hmapfile[hmapfile == -0.001] <- NA hmapfile[hmapfile == -0.001] <- NA
# Select the number of rows based on the number of genes # Select the number of rows based on the number of genes
num_total_genes <- length(hmapfile[,1]) num_total_genes <- length(hmapfile[, 1])
# break out the cluster names so each part of the cluster origin can be accessed # Break out the cluster names so each part of the cluster origin can be accessed
# line below removed because it adds to many genes to clusters when going past 1-0-10 since it cannot differentiate between 1-0-1 and 1-0-10 when using grepl. # Line below removed because it adds to many genes to clusters when going past 1-0-10
# hmapfile$cluster.origin = gsub(" ","",x=hmapfile$cluster.origin) # since it cannot differentiate between 1-0-1 and 1-0-10 when using grepl.
# hmapfile$cluster.origin = gsub(" ","",x = hmapfile$cluster.origin)
hmapfile$cluster.origin = gsub(";"," ;",x=hmapfile$cluster.origin) hmapfile$cluster.origin <- gsub(";", " ;", x = hmapfile$cluster.origin)
hmapfile$cluster.origin = strsplit(hmapfile$cluster.origin,';') hmapfile$cluster.origin <- strsplit(hmapfile$cluster.origin, ";")
# use tail(x,n) for accessing the outward most cluster # use tail(x,n) for accessing the outward most cluster
clust_rounds <- 0 clust_rounds <- 0
for(i in 1:num_total_genes){ for (i in 1:num_total_genes) {
if(length(hmapfile$cluster.origin[[i]]) > clust_rounds){ if (length(hmapfile$cluster.origin[[i]]) > clust_rounds) {
clust_rounds <- length(hmapfile$cluster.origin[[i]]) clust_rounds <- length(hmapfile$cluster.origin[[i]])
} }
} }
@@ -89,12 +91,13 @@ unique_clusts <- unique(hmapfile$cluster.origin[1:num_total_genes])
unique_clusts <- unique_clusts[unique_clusts != " "] unique_clusts <- unique_clusts[unique_clusts != " "]
#select only the unique cluster names #select only the unique cluster names
unique_clusts <- sort(unique(unlist(unique_clusts,use.names= FALSE)),decreasing=FALSE) unique_clusts <- sort(unique(unlist(unique_clusts, use.names = FALSE)), decreasing = FALSE)
num_unique_clusts <- length(unique_clusts) num_unique_clusts <- length(unique_clusts)
# Base the color key on a statistical analysis of the L and K data # Base the color key on a statistical analysis of the L and K data
# need to create "breaks" to set the color key, need to have 12 different breaks (for 11 colors) # need to create "breaks" to set the color key, need to have 12 different breaks (for 11 colors)
# scale() will calculate the mean and standard deviation of the entire vector, then "scale" each element by those values by subtracting the mean and dividing by the sd. # scale() will calculate the mean and standard deviation of the entire vector
# then "scale" each element by those values by subtracting the mean and dividing by the sd
# hmapfile[,4:(length(hmapfile[1,]) - 2)] <- scale(hmapfile[,4:(length(hmapfile[1,]) - 2)]) # hmapfile[,4:(length(hmapfile[1,]) - 2)] <- scale(hmapfile[,4:(length(hmapfile[1,]) - 2)])
@@ -107,19 +110,19 @@ L_MAX <- 0
KcolumnValues <- vector() KcolumnValues <- vector()
LcolumnValues <- vector() LcolumnValues <- vector()
for(i in 4:(length(hmapfile[1,]) - 3)){ for (i in 4:(length(hmapfile[1, ]) - 3)){
if(grepl("_Z_lm_K",colnames(hmapfile)[i],fixed=TRUE) == TRUE){ if (grepl("_Z_lm_K", colnames(hmapfile)[i], fixed = TRUE) == TRUE) {
KcolumnValues <- append(KcolumnValues,i) KcolumnValues <- append(KcolumnValues, i)
} }
if(grepl("_Z_lm_L",colnames(hmapfile)[i],fixed=TRUE) == TRUE){ if (grepl("_Z_lm_L", colnames(hmapfile)[i], fixed = TRUE) == TRUE) {
LcolumnValues <- append(LcolumnValues,i) LcolumnValues <- append(LcolumnValues, i)
} }
} }
# L_MAX <- quantile(hmapfile[,LcolumnValues],c(0,.01,.5,.99,1),na.rm=TRUE)[4] # L_MAX <- quantile(hmapfile[,LcolumnValues],c(0,.01,.5,.99,1),na.rm = TRUE)[4]
# K_MIN <- quantile(hmapfile[,KcolumnValues],c(0,.01,.5,.99,1),na.rm=TRUE)[2] # K_MIN <- quantile(hmapfile[,KcolumnValues],c(0,.01,.5,.99,1),na.rm = TRUE)[2]
# L_MAX <- quantile(hmapfile[,LcolumnValues],c(0,.01,.5,.975,1),na.rm=TRUE)[4] # L_MAX <- quantile(hmapfile[,LcolumnValues],c(0,.01,.5,.975,1),na.rm = TRUE)[4]
# K_MIN <- quantile(hmapfile[,KcolumnValues],c(0,.025,.5,.99,1),na.rm=TRUE)[2] # K_MIN <- quantile(hmapfile[,KcolumnValues],c(0,.025,.5,.99,1),na.rm = TRUE)[2]
# Z scores are # Z scores are
L_MAX <- 12 L_MAX <- 12
@@ -128,7 +131,7 @@ K_MIN <- -12
# L_Multiplier <- as.numeric(abs(K_MIN/L_MAX)) # L_Multiplier <- as.numeric(abs(K_MIN/L_MAX))
# hmapfile[,LcolumnValues] <- hmapfile[,LcolumnValues] * L_Multiplier # hmapfile[,LcolumnValues] <- hmapfile[,LcolumnValues] * L_Multiplier
# if(grepl("SHIFT",colnames(hmapfile)[4],fixed=TRUE) == TRUE){ # if(grepl("SHIFT",colnames(hmapfile)[4],fixed = TRUE) == TRUE){
# print("FOUND SHIFT VALUES") # print("FOUND SHIFT VALUES")
# hmapfile[,(LcolumnValues - 1)] <- hmapfile[,(LcolumnValues-1)] * L_Multiplier # hmapfile[,(LcolumnValues - 1)] <- hmapfile[,(LcolumnValues-1)] * L_Multiplier
# } # }
@@ -142,15 +145,19 @@ print(KEY_MIN)
print(L_MAX) print(L_MAX)
#print(L_Multiplier) #print(L_Multiplier)
colormapbreaks <- c(KEY_MIN,KEY_MIN*(5/6),KEY_MIN*(4/6),KEY_MIN*(3/6),KEY_MIN*(2/6),KEY_MIN*(1/6),KEY_MAX*(1/6),KEY_MAX*(2/6),KEY_MAX*(3/6),KEY_MAX*(4/6),KEY_MAX*(5/6),KEY_MAX) colormapbreaks <- c(KEY_MIN, KEY_MIN * (5 / 6), KEY_MIN * (4 / 6), KEY_MIN * (3 / 6),
#print(colormapbreaks) KEY_MIN * (2 / 6), KEY_MIN * (1 / 6), KEY_MAX * (1 / 6), KEY_MAX * (2 / 6),
KEY_MAX * (3 / 6), KEY_MAX * (4 / 6), KEY_MAX * (5 / 6), KEY_MAX)
# print(colormapbreaks)
# Probably should give a way to detect shift in case that is is not in the first row... (maybe just grepl for the whole column name?) # Probably should give a way to detect shift in case that is is not in the first row... (maybe just grepl for the whole column name?)
# However since also using this to amend the first part. Could possibly identify all the ones that contain the word shift and then create an object containing just those numbers # However since also using this to amend the first part.
# then could just use these values and create spaces only between interaction values - possibly could get rid of redundant shift values if we don't want to view these # Could possibly identify all the ones that contain the word shift and then create an object containing just those numbers
# then could just use these values and create spaces only between interaction values
# possibly could get rid of redundant shift values if we don't want to view these
# could we pool all the shift data/average it? # could we pool all the shift data/average it?
if(grepl("Shift",colnames(hmapfile)[4],fixed=TRUE) == TRUE){ if (grepl("Shift", colnames(hmapfile)[4], fixed = TRUE) == TRUE) {
even_columns <- seq(from= 2, to= (length(hmapfile[1,]) - 7),by=2) even_columns <- seq(from = 2, to = (length(hmapfile[1, ]) - 7), by = 2)
# ev_repeat = rep("white",length(even_columns)) # ev_repeat = rep("white",length(even_columns))
# ev_repeat = rep("red",(length(hmapfile[1,]) - 5)) # ev_repeat = rep("red",(length(hmapfile[1,]) - 5))
# middle_col <- (length(hmapfile[1,]) - 5)/2 # middle_col <- (length(hmapfile[1,]) - 5)/2
@@ -158,8 +165,8 @@ if(grepl("Shift",colnames(hmapfile)[4],fixed=TRUE) == TRUE){
# print(ev_repeat) # print(ev_repeat)
} }
if(grepl("Shift",colnames(hmapfile)[4],fixed=TRUE) == FALSE){ if (grepl("Shift", colnames(hmapfile)[4], fixed = TRUE) == FALSE) {
even_columns <- seq(from= 2, to= (length(hmapfile[1,]) - 7),by=1) even_columns <- seq(from = 2, to = (length(hmapfile[1, ]) - 7), by = 1)
print("NO SHIFT VALS FOUND") print("NO SHIFT VALS FOUND")
} }
@@ -167,43 +174,43 @@ if(grepl("Shift",colnames(hmapfile)[4],fixed=TRUE) == FALSE){
# even_columns <- c(2,5,7,10,12,15,17) # even_columns <- c(2,5,7,10,12,15,17)
# m <- 0 # m <- 0
colnames_edit <- as.character(colnames(hmapfile)[4:(length(hmapfile[1,]) - 3)]) colnames_edit <- as.character(colnames(hmapfile)[4:(length(hmapfile[1, ]) - 3)])
colnames(DAmP_list)[1] <- "ORF" colnames(DAmP_list)[1] <- "ORF"
hmapfile$DAmPs <- "YKO" hmapfile$DAmPs <- "YKO"
colnames(hmapfile)[2] <- "ORF" colnames(hmapfile)[2] <- "ORF"
try(hmapfile[hmapfile$ORF %in% DAmP_list$ORF,]$DAmPs <- "YKD") try(hmapfile[hmapfile$ORF %in% DAmP_list$ORF, ]$DAmPs <- "YKD")
# X <- X[order(X$DAmPs,decreasing = TRUE),] # X <- X[order(X$DAmPs,decreasing = TRUE),]
hmapfile$color2 <- NA hmapfile$color2 <- NA
try(hmapfile[hmapfile$DAmPs == "YKO",]$color2 <- "black") try(hmapfile[hmapfile$DAmPs == "YKO", ]$color2 <- "black")
try(hmapfile[hmapfile$DAmPs == "YKD",]$color2 <- "red") try(hmapfile[hmapfile$DAmPs == "YKD", ]$color2 <- "red")
hmapfile$color <- NA hmapfile$color <- NA
try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_many2many",]$color <- "#F8766D") try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_many2many", ]$color <- "#F8766D")
try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_one2many",]$color <- "#00BA38") try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_one2many", ]$color <- "#00BA38")
try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_one2one",]$color <- "#619CFF") try(hmapfile[hmapfile$hsapiens_homolog_orthology_type == "ortholog_one2one", ]$color <- "#619CFF")
# print(colnames_edit) # print(colnames_edit)
for(i in 1:length(colnames_edit)){ for (i in 1:length(colnames_edit)) {
if(grepl("Shift",colnames_edit[i],fixed=TRUE) == TRUE){ if (grepl("Shift", colnames_edit[i], fixed = TRUE) == TRUE) {
colnames_edit[i] <- "" colnames_edit[i] <- ""
colnames_edit[i+1] <- gsub(pattern = "_Z_lm_",replacement = " ",x = colnames_edit[i+1]) colnames_edit[i + 1] <- gsub(pattern = "_Z_lm_", replacement = " ", x = colnames_edit[i + 1])
try(colnames_edit[i+1] <- gsub(pattern = "_",replacement = " ",x = colnames_edit[i+1])) try(colnames_edit[i + 1] <- gsub(pattern = "_", replacement = " ", x = colnames_edit[i + 1]))
# INT_store <- strsplit(colnames_edit[i+1], "Z_lm") # INT_store <- strsplit(colnames_edit[i+1], "Z_lm")
# print(length(unlist(INT_store))) # print(length(unlist(INT_store)))
# if(length(unlist(INT_store)) == 4){ # if(length(unlist(INT_store)) == 4){
# colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[3],sep=" ") # colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[3],sep = " ")
# } # }
# if(length(unlist(INT_store)) == 3){ # if(length(unlist(INT_store)) == 3){
# #
# colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],sep=" ") # colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],sep = " ")
# } # }
# if(length(unlist(INT_store)) == 5){ # if(length(unlist(INT_store)) == 5){
# colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[3],unlist(INT_store)[4],sep=" ") # colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[3],unlist(INT_store)[4],sep = " ")
# } # }
# if(length(unlist(INT_store)) == 6){ # if(length(unlist(INT_store)) == 6){
# colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[6],sep=" ") # colnames_edit[i+1] <- paste(unlist(INT_store)[1],unlist(INT_store)[2],unlist(INT_store)[6],sep = " ")
# } # }
} }
@@ -217,135 +224,142 @@ print(colnames_edit)
# colnames_edit[20] <- "TEM HL L" # colnames_edit[20] <- "TEM HL L"
# Create the heatmaps # Create the heatmaps
for(i in 1:num_unique_clusts){ for (i in 1:num_unique_clusts) {
cluster <- unique_clusts[i] cluster <- unique_clusts[i]
cluster_data <- subset(hmapfile,grepl(cluster,cluster.origin)) cluster_data <- subset(hmapfile, grepl(cluster, cluster.origin))
cluster_length <- length(cluster_data[,1]) cluster_length <- length(cluster_data[, 1])
if(cluster_length != 1){ if (cluster_length != 1) {
X0 <- as.matrix(cluster_data[,4:(length(hmapfile[1,]) - 6)]) X0 <- as.matrix(cluster_data[, 4:(length(hmapfile[1, ]) - 6)])
if(cluster_length >= 2001){ if (cluster_length >= 2001) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=20,width=15) pdf(file = mypath, height = 20, width = 15)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor= "white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
# zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
# cellnote = round(X0,digits=0), notecex = 0.1, key=TRUE, # cellnote = round(X0,digits = 0), notecex = 0.1, key = TRUE,
keysize=0.7, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 0.7, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
# ColSideColors=ev_repeat, # ColSideColors = ev_repeat,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
# abline(v=0.5467,col="black") )
# abline(v = 0.5467,col = "black")
dev.off() dev.off()
} }
if(cluster_length >= 201 && cluster_length <= 2000){ if (cluster_length >= 201 && cluster_length <= 2000) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=15,width=12) pdf(file = mypath, height = 15, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
# zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.1, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.1, key = TRUE,
keysize=0.7, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 0.7, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
# abline(v=0.5316,col="black") )
# abline(v = 0.5316,col = "black")
dev.off() dev.off()
} }
if(cluster_length >= 150 && cluster_length <= 200){ if (cluster_length >= 150 && cluster_length <= 200) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=12,width=12) pdf(file = mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.1, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
# zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.2, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.2, key = TRUE,
keysize=1, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 1, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
)
dev.off() dev.off()
} }
if(cluster_length >= 101 && cluster_length <= 149){ if (cluster_length >= 101 && cluster_length <= 149) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,mypath,height=12,width=12) pdf(file = mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.2, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.2, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
# zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.3, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.3, key = TRUE,
keysize=1, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 1, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
)
dev.off() dev.off()
} }
if(cluster_length >= 60 && cluster_length <= 100){ if (cluster_length >= 60 && cluster_length <= 100) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=12,width=12) pdf(file = mypath, height = 12, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.4, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.4, scale = "none",
breaks=colormapbreaks,symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
#zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.3, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.3, key = TRUE,
keysize=1, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 1, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
)
dev.off() dev.off()
} }
if(cluster_length <= 59 && cluster_length >= 30){ if (cluster_length <= 59 && cluster_length >= 30) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=9,width=12) pdf(file = mypath, height = 9, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, distfun = dist, hclustfun = hclust, Rowv = TRUE, Colv = NA, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.6, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.6, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
#zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.4, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.4, key = TRUE,
keysize=1, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 1, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
)
dev.off() dev.off()
} }
if(cluster_length <= 29){ if (cluster_length <= 29) {
mypath = file.path(outputPath,paste("cluster_",gsub(" ","",cluster), ".pdf",sep="")) mypath <- file.path(outputPath, paste("cluster_", gsub(" ", "", cluster), ".pdf", sep = ""))
pdf(file=mypath,height=7,width=12) pdf(file = mypath, height = 7, width = 12)
heatmap.2( heatmap.2(
x=X0, x = X0,
Rowv=TRUE, Colv=NA, Rowv = TRUE, Colv = NA,
distfun = dist, hclustfun = hclust, distfun = dist, hclustfun = hclust,
dendrogram = "row", cexCol = 0.8, cexRow = 0.9, scale = "none", dendrogram = "row", cexCol = 0.8, cexRow = 0.9, scale = "none",
breaks=colormapbreaks, symbreaks=FALSE, colsep = even_columns, sepcolor="white", offsetCol = 0.1, breaks = colormapbreaks, symbreaks = FALSE, colsep = even_columns, sepcolor = "white", offsetCol = 0.1,
#zlim=c(-132,132), # zlim = c(-132,132),
xlab = "Type of Media", ylab = "Gene Name", xlab = "Type of Media", ylab = "Gene Name",
cellnote = round(X0,digits=0), notecex = 0.4, key=TRUE, cellnote = round(X0, digits = 0), notecex = 0.4, key = TRUE,
keysize=1, trace="none", density.info=c("none"), margins=c(10, 8), keysize = 1, trace = "none", density.info = c("none"), margins = c(10, 8),
na.color="red", col=brewer.pal(11,"PuOr"), na.color = "red", col = brewer.pal(11, "PuOr"),
main=cluster, main = cluster,
labRow=as.character(cluster_data$Gene), labCol=colnames_edit, colRow=cluster_data$color2,RowSideColors=cluster_data$color) labRow = as.character(cluster_data$Gene), labCol = colnames_edit, colRow = cluster_data$color2, RowSideColors = cluster_data$color
)
dev.off() dev.off()
} }
} }
# print(paste("FINISHED", "CLUSTER",cluster,sep=" ")) # print(paste("FINISHED", "CLUSTER",cluster,sep = " "))
} }