cdhit_wrapper.R

#! /bin/Rscript
library("optparse")

option_list <- list(
  make_option(opt_str = c("-i", "--input"), default = NULL, help = "Input bed-file. Fourth column is expected to contain names, last column must be sequences.", metavar = "character"),
  make_option(opt_str = c("-c", "--identity"), default = 0.8, help = "Identity threshold. Default = %default", metavar = "double >= 0.8"),
  make_option(opt_str = c("-A", "--coverage"), default = 8, help = "Minimal alignment length for both sequences in nucelotides. Default = %default", metavar = "integer"),
  make_option(opt_str = c("-o", "--output"), default = "cluster.bed", help = "Output file same as input but with appended column of cluster numbers. Default = %default", metavar = "character"),
  make_option(opt_str = c("--clean"), default = TRUE, help = "Delete all temporary files. Default = %default", metavar = "logical"),
  # make_option(opt_str = c("-G", "--global"), default = FALSE),
  # make_option(opt_str = c("-b", "--band_width"), default = 20),
  # make_option(opt_str = c("-M", "--memory"), default = 800),
  make_option(opt_str = c("-T", "--threads"), default = 1)#,
  # make_option(opt_str = c("-n", "--word_length"), default = 3),
  # make_option(opt_str = c("-l", "--throw_away_sequences"), default = 5),
  # make_option(opt_str = c("-d", "--description")), # can not produce bed if this is != 0
  # make_option(opt_str = c("-s", "--length_dif_cutoff_shorter_p"), default = 0.0),
  # make_option(opt_str = c("-S", "--length_dif_cutoff_shorter_n"), default = 999999),
  # make_option(opt_str = c("-aL", "--alignment_coverage_longer_p"), default = 0.0),
  # make_option(opt_str = c("-AL", "--alignment_coverage_longer_n"), default = 99999999),
  # make_option(opt_str = c("-aS", "--alignment_coverage_shorter_p"), default = 0.0),
  # make_option(opt_str = c("-AS", "--alignment_coverage_shorter_n"), default = 99999999),
  # make_option(opt_str = c("-uL", "--max_unmatched_longer_p"), default = 1.0),
  # make_option(opt_str = c("-uS", "--max_unmatched_shorter_p"), default = 99999999),
  # make_option(opt_str = c("-U", "--max_unmatched_both_n"), default = 99999999),
  # make_option(opt_str = c("-g", "--fast_cluster"), default = FALSE),
  # make_option(opt_str = c("-r", "--ignore_strand"), default = TRUE),
  # make_option(opt_str = c("--match"), default = 2),
  # make_option(opt_str = c("--mismatch"), default = -2),
  # make_option(opt_str = c("--gap"), default = -6),
  # make_option(opt_str = c("--gap_ext"), default = -1),
  # make_option(opt_str = c("-sc", "--sort_cluster_by_size"), default = TRUE)
  # TODO more args
)

opt_parser <- OptionParser(option_list = option_list, 
                           description = "CD-HIT-EST Wrapper function.")

opt <- parse_args(opt_parser)

#' cd-hit wrapper
#' 
#' @param input Data.table or file in bed format (requires names in fourth column and sequences in last column).
#' @param identity Identity threshold.
#' @param coverage Minimal alignment length for both sequences in nucelotides.
#' @param output Clustered bedfile. Adds cluster number in last column (lower number = bigger cluster).
#' @param clean Clean up after run.
#' @param threads Number of threads to use in cd-hit (0 = all cpus).
#' 
#' TODO add all cdhit parameter (necessary ?)
#' TODO check whether cdhit is installed
cdhitest <- function(input, identity = 0.8, coverage = 8, output = "cluster.bed", clean = TRUE, threads = 1) {
  message("Loading bed.")
  # load bed if neccessary
  if (!data.table::is.data.table(input)) {
    bed_table <- data.table::fread(input = input, header = FALSE)
  } else {
    bed_table <- input
  }
  
  message("Convert to fasta.")
  ### convert bed to fasta
  # 4th column = name
  # last column = sequence
  fasta_file <- "converted_bed.fasta"
  seqinr::write.fasta(sequences = as.list(bed_table[[ncol(bed_table)]]), names = bed_table[[4]], as.string = TRUE, file.out = fasta_file)
  
  message("Clustering.")
  ### cd-hit-est
  cdhit_output <- "cdhit_output"
  cdhit_call <- paste("cd-hit-est -i", fasta_file, "-o", cdhit_output, "-c", identity, "-A", coverage, "-T ", threads, "-G 0 -n 3 -g 1 -r 0 -l 5 -sc 1 -d 0")
  
  system(command = cdhit_call, wait = TRUE)
  
  message("Formatting/ writing results.")
  # reformat cluster file
  # columns: id, clstr, clstr_size, length, clstr_rep, clstr_iden, clstr_cov
  cluster_file <- "reformated.clstr"
  cluster_call <- paste("clstr2txt.pl", paste0(cdhit_output, ".clstr"), ">", cluster_file)
  
  system(command = cluster_call, wait = TRUE)
  
  # load reformated file
  cluster <- data.table::fread(cluster_file)
  
  ### add cluster to bed_table
  result <- merge(x = bed_table, y = cluster[, c("id", "clstr")], by.x = "V4", by.y = "id", sort = FALSE)[, union(names(bed_table), names(cluster)[2]), with = FALSE]
  
  # delete files
  if (clean) {
    file.remove(fasta_file, paste0(cdhit_output, ".clstr"), cdhit_output, cluster_file)
  }
  
  data.table::fwrite(x = result, file = output, sep = "\t", col.names = FALSE)
}

# call function with given parameter if not in interactive context (e.g. run from shell)
if (!interactive()) {
  # remove last parameter (help param)
  params <- opt[-length(opt)]
  do.call(cdhitest, args = params)
}