Feature/Pass additional args to harmony

src/harmony/bin/run_harmony.R

@@ -1,8 +1,10 @@
 #!/usr/bin/env Rscript
 
-print("##################################################")
-print("# Harmony: Algorithm for single cell integration #")
-print("##################################################")
+print("#############################################################")
+print("# Harmony: Algorithm for single cell integration            #")
+print('# GitHub: https://github.com/immunogenomics/harmony         #')
+print('# Paper: https://www.nature.com/articles/s41592-019-0619-0  #')
+print("#############################################################")
 
 # Loading dependencies scripts
 library("argparse")
@@ -21,6 +23,13 @@ parser$add_argument(
     default = "foo",
   	help="Prefix path to save output files. [default %default]"
 )
+parser$add_argument(
+    '--seed',
+    type="character",
+    dest='seed',
+    default=617,
+    help='Seed. [default %default]'
+)
 parser$add_argument(
     "--vars-use",
     type="character",
@@ -29,13 +38,6 @@ parser$add_argument(
 	  default=NULL,
     help='If meta_data is dataframe, this defined which variable(s) to remove (character vector).'
 )
-parser$add_argument(
-    '--seed',
-    type="character",
-    dest='seed',
-    default=617,
-    help='Seed. [default %default]'
-)
 parser$add_argument(
 	'--do-pca',
 	type="logical",
@@ -44,9 +46,36 @@ parser$add_argument(
 	default=FALSE,
 	help='Whether to perform PCA on input matrix.'
 )
+parser$add_argument(
+    '--theta',
+    type="double",
+    dest='theta',
+    default=NULL,
+    help='Diversity clustering penalty parameter. Specify for each variable in vars_use Default theta=2. theta=0 does not encourage any diversity. Larger values of theta result in more diverse clusters. [default %default]'
+)
+parser$add_argument(
+    '--lambda',
+    type="double",
+    dest='lambda',
+    default=NULL,
+    help='Ridge regression penalty parameter. Specify for each variable in vars_use. Default lambda=1. Lambda must be strictly positive. Smaller values result in more aggressive correction. [default %default]'
+)
+parser$add_argument(
+    '--epsilon-harmony',
+    type="double",
+    dest='epsilon_harmony',
+    default=1e-04,
+    help='Convergence tolerance for Harmony. Set to -Inf to never stop early. [default %default]'
+)
+
 
 args <- parser$parse_args()
 
+if(args$epsilon_harmony < 0) {
+  args$epsilon_harmony <- -Inf
+  print("Setting epsilon.harmony argument to -Inf...")
+}
+
 cat("Parameters: \n")
 print(args)
 
@@ -61,6 +90,19 @@ if(!is.null(args$seed)) {
   warnings("No seed is set, this will likely give none reproducible results.")
 }
 
+# Required for reproducibility in case numeric parameters are passed (e.g.: theta, lambda)
+args <- lapply(X = args, FUN = function(arg) {
+  if(is.numeric(x = arg)) {
+    if(arg %% 1 == 0) {
+      return (as.integer(x = arg))
+    } else {
+      return (arg)
+    }
+  }
+  return (arg)
+})
+
+
 input_ext <- tools::file_ext(args$input)
 
 if(input_ext == "h5ad") {
@@ -84,7 +126,7 @@ if(input_ext == "h5ad") {
   stop(paste0("Unrecognized input file format: ", input_ext, "."))
 }
 
-print(paste0("PCA embeddings matrix has ", dim(x = data)[1], " rows, ", dim(x = data)[2], " columns."))
+print(paste0("PCA embeddings matrix has ", dim(x = pca_embeddings)[1], " rows, ", dim(x = pca_embeddings)[2], " columns."))
 
 if(sum(args$vars_use %in% colnames(x = metadata)) != length(x = args$vars_use)) {
 	stop("Some argument value from the parameter(s) --vars-use are not found in the metadata.")
@@ -99,6 +141,9 @@ harmony_embeddings <- harmony::HarmonyMatrix(
   meta_data = metadata,
   vars_use = args$vars_use,
   do_pca = args$do_pca,
+  theta = args$theta,
+  lambda = args$lambda,
+  epsilon.harmony = args$epsilon_harmony,
   verbose = FALSE
 )
 

src/harmony/harmony.config

@@ -4,6 +4,9 @@ params {
             container = 'vibsinglecellnf/harmony:1.0-1'
             report_ipynb = "${params.misc.test.enabled ? '../../..' : ''}/src/harmony/bin/reports/sc_harmony_report.ipynb"
             varsUse = ['batch']
+            // theta = ''
+            // lambda = ''
+            // epsilonHarmony = '' // Set to -1 for you need to set this parameter to -Inf in HarmonyMatrix R function
         }
     }
 }

src/harmony/processes/runHarmony.nf

@@ -20,13 +20,17 @@ process SC__HARMONY__HARMONY_MATRIX {
 
     script:
         def sampleParams = params.parseConfig(sampleId, params.global, params.sc.harmony)
-		processParams = sampleParams.local
-        varsUseAsArguments = processParams.varsUse.collect({ '--vars-use' + ' ' + it }).join(' ')
+		def processParams = sampleParams.local
+        // Arguments
+        def varsUseAsArguments = processParams.varsUse.collect({ '--vars-use' + ' ' + it }).join(' ')
         """
         ${binDir}run_harmony.R \
             ${f} \
             --seed ${params.global.seed} \
             ${varsUseAsArguments} \
+            ${processParams?.theta ? "--theta "+ processParams.theta : "" } \
+            ${processParams?.lambda ? "--lambda "+ processParams.lambda : "" } \
+            ${processParams?.epsilonHarmony ? "--epsilon-harmony "+ processParams.epsilonHarmony : "" } \
             --output-prefix "${sampleId}.SC__HARMONY__HARMONY_MATRIX"
         """