R/plotEmap.R

In BaderLab/POPPATHR: Population-based pathway analysis of SNP-SNP coevolution

#' Plots an EnrichmentMap in Cytoscape to visualize selection-enriched pathways
#'
#' @description Create a network where nodes are pathway gene sets significantly
#'		enriched in two independent population-based GSEA (i.e., validated
#' 		selection-enriched gene sets.) Pathways are then clustered to identify
#'		main pathways and themes of the selection-enriched gene sets.
#' @param gmt_file (char) file path to GMT file (this needs to be absolute path name).
#' @param EM_file (char) file path to results_EM.txt (generated by writeEmapFile.R).
#' @param EM_group_file (char) file path to write pathway groupings as determined
#' 		by EnrichmentMap.
#' @param output_folder (char) file path to write gene set grouping file. This will be
#' 		used for between-pathway SNP association analysis.
#' @param net_name (char) name for network in Cytoscape.
#' @param image_format (char) one of PNG, PDF, SVG, or JPEG.
#'
#' # Not run because requires Cytoscape to be installed and open.
#' # plotEmap(gmt_file=gmt_file, EM_file=EM_file, output_folder=enrich_folder,
#'		net_name="generic", image_format="png").
#' @return Filename of image to which EnrichmentMap is exported. Also side
#' 		effect of plotting the EnrichmentMap in an open session of Cytoscape.
#' @export
#'

plotEmap <- function(gmt_file, EM_file, EM_group_file, output_folder,
										 net_name="generic", image_format="png") {

	# Confirm that Cytoscape is installed and opened
	cytoscapePing()

  # Create EM using given parameters
	if (net_name %in% getNetworkList()) {
		deleteNetwork(net_name)
	}

	cat("* Building the network\n")
	em_command <- paste('enrichmentmap build analysisType="generic"',
			'gmtFile=', gmt_file,
			'enrichmentsDataset1=', EM_file,
			'pvalue=', 1,
			'qvalue=', 1,
			'similaritycutoff=', 0.375,
			'coefficients=', "COMBINED",
		  'combinedConstant=', 0.5)
  response <- commandsGET(em_command)
	renameNetwork(net_name, getNetworkSuid())

 	# Annotate the network using AutoAnnotate app
	cat("* Annotating the network using AutoAnnotate\n")
  aa_command <- paste("autoannotate annotate-clusterBoosted",
		"clusterAlgorithm=MCL",
		"maxWords=3",
		"network=", net_name)
	#print(aa_command)
	response <- commandsGET(aa_command)

	# Apply style
	cat("* Creating or applying style\n")
  styleName <- "EMapStyle"
	defaults <- list("NODE_SHAPE"="ellipse",
  			"NODE_SIZE"=30,
  			"EDGE_TRANSPARENCY"=200,
  			"NODE_TRANSPARENCY"=255,
				"EDGE_STROKE_UNSELECTED_PAINT"="#999999") # get node colour in paper

	createVisualStyle(styleName, defaults)
	setVisualStyle(styleName)
	layoutNetwork("attributes-layout NodeAttribute=__mclCLuster")

	redraw_command <- sprintf("autoannotate redraw network=%s", getNetworkSuid())
  response <- commandsGET(redraw_command)
	fitContent()

	# Pull node names from AutoAnnotated clusters
	# Get node and edge table that have the cluster annotations
	node_table <- getTableColumns(table="node", network=getNetworkSuid())
	edge_table <- getTableColumns(table="edge", network=getNetworkSuid())

	# Get the correct attribute names
	descr_attrib <- colnames(node_table)[grep(colnames(node_table), pattern="GS_DESCR")]

	# Get all the cluster numbers
	cluster_num <- node_table$'__mclCluster'

	# Get the unique set of clusters
	set_clusters <- unique(cluster_num)
	set_clusters <- set_clusters[which(set_clusters != 0)]

	# Calculate the cluster names using AutoAnnotate and collate all the nodes
	# associated with each cluster
	cluster_names <- c()
	for (i in 1:length(set_clusters)) {
	  current_cluster <- set_clusters[i]
	  gs_cluster <- node_table$name[which(node_table$'__mclCluster' == current_cluster)]

	  # For this cluster of gs get the gs descr to use in defining in AutoAnnotate
	  gs_cluster_suid <- node_table$SUID[which(node_table$name %in% gs_cluster)]
	  suids_aa <- paste("SUID", gs_cluster_suid, sep=":")

	  # Get the annotation for the given cluster
	  curernt_name = NULL
	  aa_label <- paste("autoannotate label-clusterBoosted",
				"labelColumn=", descr_attrib,
				"maxWords=3",
				"nodeList=", paste(suids_aa, collapse=","))

	  current_name <- commandsGET(aa_label)
	  cluster_names <- rbind(cluster_names,
			c(current_cluster, current_name, length(gs_cluster), paste(gs_cluster, collapse=";"))
		)
	}

	# Create list of cluster names and associated pathways
	EM_group_list <- list()
	for (k in 1:nrow(cluster_names)) {
		EM_group_list[[k]] <- unlist(strsplit(cluster_names[k,4], ";"))
		names(EM_group_list)[k] <- gsub(" ", "_", toupper(cluster_names[k,2]))
	}
	# Get single unclustered gene sets and add to list
	single_path <- node_table[-which(node_table$name %in% unlist(EM_group_list)),]$name

	# Remove annotation info from gene set names
	single_path_name <- gsub("\\%.*", "", single_path)
	single_path_name <- gsub(" ", "_", single_path_name)
	names(single_path) <- single_path_name

	# Add to existing list
	EM_group_list <- c(EM_group_list, single_path)

	# Write out rda file with pathway EM groupings
	save(EM_group_list, file=EM_group_file)

	# Write out EnrichmentMap to EM_file directory
	out_file <- gsub("\\..*", "", EM_file)
	cat(sprintf("* Writing out enrichment map image to %s\n", out_file))
	exportImage(out_file, type=image_format)
}