merge_similar_clusters.R

#!/usr/bin/env Rscript
if (!require(optparse)) install.packages("optparse"); library(optparse)
if (!require(igraph)) install.packages("igraph"); library(igraph)

option_list <- list(
  make_option(opt_str = c("-i", "--input"), default = NULL, help = "Input TSV-file. Output from merged tomtom results", metavar = "character"),
  make_option(opt_str = c("-l", "--list"), default = NULL, help = "Numerically sorted whitespace separated list of absolute fasta-file paths", metavar = "character"),
  make_option(opt_str = c("-w", "--min"), default = NULL, help = "Minimum weight of edge allowed in graph clusters.", metavar = "character")
)

opt_parser <- OptionParser(option_list = option_list,
                           description = "Adding Cluster ID to Query_ID Column",
                           epilogue = "Author: Rene Wiegandt <Rene.Wiegandt@mpi-bn.mpg.de>")

opt <- parse_args(opt_parser)

#' Merging FASTA-files, which motifs are similar.
#'
#' @parameter tsv_in <string> Path to TSV file generated by Tomtom.
#'                            The input for Tomtom is a from all clusters merged meme-file.
#' @parameter file_list <string> Numerically sorted comma separated list of absolute fasta-file paths
#' @parameter min_weight <INT> Minimum weight of edge allowed in graph clusters.
#'
#' @author René Wiegandt
#' @contact rene.wiegandt(at)mpi-bn.mpg.de
merge_similar <- function(tsv_in, file_list, min_weight){

  files <- unlist(strsplit(file_list, ","))

  # split the string on the character "." in the first to columns and safe the last value each, to get the number of the cluster.
  tsv <- data.table::fread(tsv_in, header = TRUE, sep = "\t",colClasses = 'character')
  query_cluster <- vapply(strsplit(tsv[["Query_ID"]],"\\."), function(l){
    tail(l, n = 1)
  },"")
  target_cluster <- vapply(strsplit(tsv[["Target_ID"]],"\\."), function(l){
    tail(l, n = 1)
  },"")

  # create data.table with only the cluster-numbers
  sim_not_unique <- data.table::data.table(as.numeric(query_cluster),as.numeric(target_cluster))

  # remove rows if column 1 is identical to column 2
  edgelist <- sim_not_unique[query_cluster != target_cluster]
  unique_edgelist <- unique(edgelist)
  
  # create graph from data.frame
  g <- graph_from_edgelist(as.matrix(unique(edgelist)), directed = F)
  g_tidy <- delete.vertices(simplify(g), degree(g) == 0)
  g_tidy <- set_vertex_attr(g_tidy,"name",V(g_tidy), value = (sort(unique(c(unique_edgelist$V1,unique_edgelist$V2)))) - 1)

  # plotting
  png('motif_clusters.png')
  plot(g_tidy)
  dev.off()
  clust <- clusters(g)
  if (clust$no < 1) {
    b <- lapply(files, function(f){
      system(paste("cat",f,">",basename(f)))
    })
  }
  clust_out <- file("cluster.txt")
  # merge FASTA-files depending on the clustered graphs
  a <- vapply(seq(from = 1, to = clust$no, by = 1), function(i){
    #get vector with cluster ids, which are clustered together in "motif cluster" i
    cl <- as.vector(which(clust$membership %in% c(i)))
    #create string, which represents a list, containing all FASTA-files to be merged
    fasta_list <- paste(files[cl], collapse = " ")
    #create name for merged FASTA-file
    name <- paste0("Cluster_",i - 1 ,".fasta")
    #merge fasta files
    system(paste("cat",fasta_list,">",name))
    paste("Cluster",i - 1,"is generated by merging cluster:",paste(cl - 1,collapse = ","))
  }, "")
  writeLines(a[a != ""], con = clust_out)
}
  


# run function merge_similar with given parameteres if not in interactive context (e.g. run from shell)
if (!interactive()) {
  merge_similar(opt$input, opt$list, opt$min)
}