diff --git a/Jenkinsfile b/Jenkinsfile
index ccb1f52..1e8f4e4 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -18,6 +18,7 @@ pipeline {
                              steps {
                                  echo 'Building R package....'
                                  sh "cd R && R CMD build TARGET.QC/ && curl -v --user ${NEXUS_USER}:${NEXUS_PASS} --upload-file TARGET.QC_0.0.1.tar.gz http://zsibio.ii.pw.edu.pl/nexus/repository/r-zsibio/src/contrib/TARGET.QC_0.0.1.tar.gz"
+                                 sh "cd R && R CMD build REFERENCE.SAMPLE.SET.SELECTOR/ && curl -v --user ${NEXUS_USER}:${NEXUS_PASS} --upload-file REFERENCE.SAMPLE.SET.SELECTOR_0.0.1.tar.gz http://zsibio.ii.pw.edu.pl/nexus/repository/r-zsibio/src/contrib/REFERENCE.SAMPLE.SET.SELECTOR_0.0.1.tar.gz"
                                  sh "cd R && R CMD build CODEXCOV/ && curl -v --user ${NEXUS_USER}:${NEXUS_PASS} --upload-file CODEXCOV_0.0.1.tar.gz http://zsibio.ii.pw.edu.pl/nexus/repository/r-zsibio/src/contrib/CODEXCOV_0.0.1.tar.gz"
                                  sh "cd R && R CMD build EXOMEDEPTHCOV/ && curl -v --user ${NEXUS_USER}:${NEXUS_PASS} --upload-file EXOMEDEPTHCOV_0.0.1.tar.gz http://zsibio.ii.pw.edu.pl/nexus/repository/r-zsibio/src/contrib/EXOMEDEPTHCOV_0.0.1.tar.gz"
                                  sh "cd R && R CMD build CNVCALLER.RUNNER/ && curl -v --user ${NEXUS_USER}:${NEXUS_PASS} --upload-file CNVCALLER.RUNNER_0.0.1.tar.gz http://zsibio.ii.pw.edu.pl/nexus/repository/r-zsibio/src/contrib/CNVCALLER.RUNNER_0.0.1.tar.gz"
diff --git a/R/CNVCALLER.RUNNER/R/CODEXCOV_wrapper.R b/R/CNVCALLER.RUNNER/R/CODEXCOV_wrapper.R
index e4636ac..5fa7645 100644
--- a/R/CNVCALLER.RUNNER/R/CODEXCOV_wrapper.R
+++ b/R/CNVCALLER.RUNNER/R/CODEXCOV_wrapper.R
@@ -1,24 +1,16 @@
 library('CODEXCOV')
 
-#' Function Description
-#'
-#' Function description.
-#' @param K_from
-#' @param K_to
-#' @param lmax
-#' @param cov_table
-#' @keywords 
-#' @export
-#' @examples
-#' run_wrapper_CODEXCOV
 run_wrapper_CODEXCOV <- function(K_from,
                                  K_to,
                                  lmax,
+                                 reference_set_select_method,
+                                 num_of_samples_in_reference_set,
                                  cov_table){
   calls <- run_CODEXCOV(strtoi(K_from),
                         strtoi(K_to),
                         strtoi(lmax),
-                        cov_table
-  )
+                        reference_set_select_method,
+                        strtoi(num_of_samples_in_reference_set),
+                        cov_table)
   calls
 }
diff --git a/R/CNVCALLER.RUNNER/R/EXOMEDEPTHCOV_wrapper.R b/R/CNVCALLER.RUNNER/R/EXOMEDEPTHCOV_wrapper.R
index a247bfa..27736ae 100644
--- a/R/CNVCALLER.RUNNER/R/EXOMEDEPTHCOV_wrapper.R
+++ b/R/CNVCALLER.RUNNER/R/EXOMEDEPTHCOV_wrapper.R
@@ -1,7 +1,10 @@
 library('EXOMEDEPTHCOV')
 
-run_wrapper_EXOMEDEPTHCOV <- function(cov_table){
-  calls <- run_EXOMEDEPTHCOV(cov_table
-  )
+run_wrapper_EXOMEDEPTHCOV <- function(reference_set_select_method,
+                                      num_of_samples_in_reference_set,
+                                      cov_table){
+  calls <- run_EXOMEDEPTHCOV(reference_set_select_method,
+                             num_of_samples_in_reference_set,
+                             cov_table)
   calls
 }
diff --git a/R/CNVCALLER.RUNNER/inst/run_cnvcaller.R b/R/CNVCALLER.RUNNER/inst/run_cnvcaller.R
index 1183a0c..19c88a7 100755
--- a/R/CNVCALLER.RUNNER/inst/run_cnvcaller.R
+++ b/R/CNVCALLER.RUNNER/inst/run_cnvcaller.R
@@ -33,6 +33,8 @@ read_parameters <- function(tabName, id, conn){
   K_to <- parameters[1,'k_to']
   lmax <- parameters[1,'lmax']
   chr <- parameters[1,'chr']
+  reference_set_select_method <- parameters[1,'reference_set_select_method']
+  num_of_samples_in_reference_set <- parameters[1,'num_of_samples_in_reference_set']
   scenario_id  <- parameters[1,'scenario_id']
 
   return(list(caller=caller, 
@@ -48,6 +50,8 @@ read_parameters <- function(tabName, id, conn){
               K_to=K_to, 
               lmax=lmax,
               chr=chr,
+              reference_set_select_method=reference_set_select_method,
+              num_of_samples_in_reference_set=num_of_samples_in_reference_set,
               scenario_id=scenario_id))
 }
 
@@ -106,12 +110,15 @@ run_caller <- function(parameters, cov_table){
     calls <- run_wrapper_CODEXCOV(parameters$K_from,
                                   parameters$K_to,
                                   parameters$lmax,
+                                  parameters$reference_set_select_method,
+                                  parameters$num_of_samples_in_reference_set,
                                   cov_table
     )
     calls
   } else if (parameters$caller == "exomedepth"){
-    calls <- run_wrapper_EXOMEDEPTHCOV(cov_table
-    )
+    calls <- run_wrapper_EXOMEDEPTHCOV(parameters$reference_set_select_method,
+                                       parameters$num_of_samples_in_reference_set,
+                                       cov_table)
     calls
   } else if(parameters$caller == "xhmm") {
   }
diff --git a/R/CODEXCOV/DESCRIPTION b/R/CODEXCOV/DESCRIPTION
index f764af5..4f63e70 100755
--- a/R/CODEXCOV/DESCRIPTION
+++ b/R/CODEXCOV/DESCRIPTION
@@ -14,7 +14,8 @@ Depends:
     devtools (>= 1.13.2),
     DBI (== 0.7),
     optparse (== 1.4.4),
-    CODEX (>= 1.8.0)
+    CODEX (>= 1.8.0),
+    REFERENCE.SAMPLE.SET.SELECTOR (>= 0.0.1)
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
diff --git a/R/CODEXCOV/R/run_CODEXCOV.R b/R/CODEXCOV/R/run_CODEXCOV.R
index 3562acd..78929de 100644
--- a/R/CODEXCOV/R/run_CODEXCOV.R
+++ b/R/CODEXCOV/R/run_CODEXCOV.R
@@ -9,6 +9,8 @@
 run_CODEXCOV <- function(K_from,
                          K_to,
                          lmax,
+                         reference_set_select_method,
+                         num_of_samples_in_reference_set,
                          cov_table){
   
   sampname <- unique(cov_table[,"sample_name"])
diff --git a/R/EXOMEDEPTHCOV/DESCRIPTION b/R/EXOMEDEPTHCOV/DESCRIPTION
index e99ae52..2ef92a8 100644
--- a/R/EXOMEDEPTHCOV/DESCRIPTION
+++ b/R/EXOMEDEPTHCOV/DESCRIPTION
@@ -15,7 +15,8 @@ Depends:
     DBI (== 0.7),
     optparse (== 1.4.4),
     IRanges (>= 2.0.0),
-    ExomeDepth (>= 1.1.10)
+    ExomeDepth (>= 1.1.10),
+    REFERENCE.SAMPLE.SET.SELECTOR (>= 0.0.1)
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
diff --git a/R/EXOMEDEPTHCOV/R/run_EXOMEDEPTHCOV.R b/R/EXOMEDEPTHCOV/R/run_EXOMEDEPTHCOV.R
index 45fcf8c..c7fc967 100644
--- a/R/EXOMEDEPTHCOV/R/run_EXOMEDEPTHCOV.R
+++ b/R/EXOMEDEPTHCOV/R/run_EXOMEDEPTHCOV.R
@@ -1,8 +1,10 @@
 library(ExomeDepth)
 library(methods)
 
-run_EXOMEDEPTHCOV <- function(cov_table){
-  
+run_EXOMEDEPTHCOV <- function(reference_set_select_method,
+                              num_of_samples_in_reference_set,
+                              cov_table){
+
   sampname <- unique(cov_table[,"sample_name"])
   targets <- cov_table[,c("target_id", "chr", "pos_min", "pos_max")]
   targets <- targets[!duplicated(targets[,"target_id"]),]
@@ -26,13 +28,14 @@ run_EXOMEDEPTHCOV <- function(cov_table){
       for (i in 1:nrow(Y)) {
         target_length <- c(target_length, width(ref[i]))
       }
-      reference_set <- select.reference.set (test.counts = Y[,actual_sample_id],
-                                             reference.counts = Y[,-actual_sample_id],
-                                             bin.length = target_length,
-                                             n.bins.reduced = 10000)
+      reference_samples <- run_REFERENCE.SAMPLE.SET.SELECTOR(actual_sample,
+                                                             Y,
+                                                             reference_set_select_method,
+                                                             num_of_samples_in_reference_set,
+                                                             target_length)
 
       ## ----construct.ref-------------------------------------------------------
-      my.matrix <- as.matrix(Y[,reference_set$reference.choice])
+      my.matrix <- as.matrix(Y[,reference_samples])
       my.reference.selected <- apply(X = my.matrix, 
                                      MAR = 1, 
                                      FUN = sum)
diff --git a/R/REFERENCE.SAMPLE.SET.SELECTOR/DESCRIPTION b/R/REFERENCE.SAMPLE.SET.SELECTOR/DESCRIPTION
new file mode 100755
index 0000000..ffb4b8e
--- /dev/null
+++ b/R/REFERENCE.SAMPLE.SET.SELECTOR/DESCRIPTION
@@ -0,0 +1,16 @@
+Package: REFERENCE.SAMPLE.SET.SELECTOR
+Title: Package For Selecting Reference Sample Set
+Version: 0.0.1
+Authors@R: c(
+    person("Tomasz", "Gambin", email = "tgambin@gmail.com", role = c("aut", "cre")),
+    person("Marek", "Wiewiórka", email = "marek.wiewiorka@gmail.com", role = c("aut")),
+    person("Wiktor", "Kuśmirek", email = "kusmirekwiktor@gmail.com", role = c("aut")))
+Description: An implementation of the package with multiple function
+    for selecting reference sample set.
+Depends:
+    R (>= 3.2.3),
+    ExomeDepth (>= 1.1.10)
+License: GPL-3
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 6.0.1.9000
diff --git a/R/REFERENCE.SAMPLE.SET.SELECTOR/NAMESPACE b/R/REFERENCE.SAMPLE.SET.SELECTOR/NAMESPACE
new file mode 100755
index 0000000..884a631
--- /dev/null
+++ b/R/REFERENCE.SAMPLE.SET.SELECTOR/NAMESPACE
@@ -0,0 +1,2 @@
+# Generated by roxygen2: fake comment so roxygen2 overwrites silently.
+exportPattern("^[^\\.]")
diff --git a/R/REFERENCE.SAMPLE.SET.SELECTOR/R/functions_REFERENCE.SAMPLE.SET.SELECTOR.R b/R/REFERENCE.SAMPLE.SET.SELECTOR/R/functions_REFERENCE.SAMPLE.SET.SELECTOR.R
new file mode 100644
index 0000000..6ac3736
--- /dev/null
+++ b/R/REFERENCE.SAMPLE.SET.SELECTOR/R/functions_REFERENCE.SAMPLE.SET.SELECTOR.R
@@ -0,0 +1,32 @@
+library(ExomeDepth)
+
+canoes_method <- function(investigated_sample, Y, num_refs){
+  if (num_refs == 0) {
+    num_refs <- 30  # in CANOES application num_refs is default set to 30
+  }
+  samples <- colnames(Y)
+  cov <- cor(Y[, samples], Y[, samples])
+  reference_samples <- setdiff(samples, investigated_sample)
+  covariances <- cov[investigated_sample, reference_samples]
+  reference_samples <- names(sort(covariances, 
+          decreasing=T)[1:min(num_refs, length(covariances))])
+  return(list(reference_samples=reference_samples))
+}
+
+exomedepth_method <- function(investigated_sample, Y, num_refs, target_length){
+  samples <- colnames(Y)
+  reference_samples <- setdiff(samples, investigated_sample)
+  reference_set <- select.reference.set(test.counts = Y[,investigated_sample],
+                                        reference.counts = Y[,reference_samples],
+                                        bin.length = target_length,
+                                        n.bins.reduced = 10000)
+  if (num_refs == 0) {
+    reference_samples <- reference_set$reference.choice
+  } else {
+    reference <- reference_set$summary.stats[1:num_refs,'ref.samples']
+    reference_samples <- c()
+    for (s in reference)
+      reference_samples <- c(reference_samples, c(s))
+  } 
+  return(list(reference_samples=reference_samples))
+}
diff --git a/R/REFERENCE.SAMPLE.SET.SELECTOR/R/run_REFERENCE.SAMPLE.SET.SELECTOR.R b/R/REFERENCE.SAMPLE.SET.SELECTOR/R/run_REFERENCE.SAMPLE.SET.SELECTOR.R
new file mode 100644
index 0000000..f025b95
--- /dev/null
+++ b/R/REFERENCE.SAMPLE.SET.SELECTOR/R/run_REFERENCE.SAMPLE.SET.SELECTOR.R
@@ -0,0 +1,16 @@
+run_REFERENCE.SAMPLE.SET.SELECTOR <- function(investigated_sample,
+                                              Y,
+                                              select_method,
+                                              num_refs,
+                                              target_length){
+  if(select_method == "canoes") {
+    reference_samples <- canoes_method(investigated_sample, Y, num_refs)$reference_samples
+  } else if(select_method == "codex") {
+    #reference_samples <- codex_method(investigated_sample, Y, num_refs)$reference_samples
+  } else if(select_method == "exomedepth") {
+    reference_samples <- exomedepth_method(investigated_sample, Y, num_refs, target_length)$reference_samples
+  } else if(select_method == "clamms") {
+    #reference_samples <- clamms_method(investigated_sample, Y, num_refs)$reference_samples
+  }
+  reference_samples
+}
diff --git a/R/tests/run_tests.R b/R/tests/run_tests.R
index 3db6e0d..8b5d86f 100755
--- a/R/tests/run_tests.R
+++ b/R/tests/run_tests.R
@@ -7,6 +7,7 @@ if (length(which(installed.packages()[,1] == "CODEX")) == 0){
 if (length(which(installed.packages()[,1] == "ExomeDepth")) == 0){install.packages("ExomeDepth",repos="https://cloud.r-project.org/")}
 
 if (length(which(installed.packages()[,1] == "TARGET.QC")) > 0){remove.packages("TARGET.QC")}
+if (length(which(installed.packages()[,1] == "REFERENCE.SAMPLE.SET.SELECTOR")) > 0){remove.packages("REFERENCE.SAMPLE.SET.SELECTOR")}
 if (length(which(installed.packages()[,1] == "CODEXCOV")) > 0){remove.packages("CODEXCOV")}
 if (length(which(installed.packages()[,1] == "EXOMEDEPTHCOV")) > 0){remove.packages("EXOMEDEPTHCOV")}
 if (length(which(installed.packages()[,1] == "CNVCALLER.RUNNER")) > 0){remove.packages("CNVCALLER.RUNNER")}
@@ -14,6 +15,7 @@ if (length(which(installed.packages()[,1] == "CNVCALLER.EVALUATOR")) > 0){remove
 
 setwd('tests/')
 devtools::install('../TARGET.QC')
+devtools::install('../REFERENCE.SAMPLE.SET.SELECTOR')
 devtools::install('../CODEXCOV')
 devtools::install('../EXOMEDEPTHCOV')
 devtools::install('../CNVCALLER.RUNNER')
diff --git a/R/tests/test_CODEXCOV_wrapper.R b/R/tests/test_CODEXCOV_wrapper.R
index c9dadf2..9fb945e 100644
--- a/R/tests/test_CODEXCOV_wrapper.R
+++ b/R/tests/test_CODEXCOV_wrapper.R
@@ -7,6 +7,8 @@ test_that("basic test for run_wrapper_CODEXCOV function, without calls",{
   K_from <- "1"
   K_to <- "3"
   lmax <- "200"
+  reference_set_select_method <- "codex"
+  num_of_samples_in_reference_set <- "0"
   cov_table <- as.data.frame(matrix(data=c("NA12044",6,"2",17211,173310,151.636363636364,
                                            "NA11829",4,"2",16712,167190,100,
                                            "NA07051",7,"2",30275,304310,155.929936305732,
@@ -67,6 +69,8 @@ test_that("basic test for run_wrapper_CODEXCOV function, without calls",{
   calls <- run_wrapper_CODEXCOV(K_from,
                                 K_to,
                                 lmax,
+                                reference_set_select_method,
+                                num_of_samples_in_reference_set,
                                 cov_table)
   expect_equal(length(calls), 0)
 })
diff --git a/R/tests/test_EXOMEDEPTHCOV_wrapper.R b/R/tests/test_EXOMEDEPTHCOV_wrapper.R
index 9e61c22..15753e9 100644
--- a/R/tests/test_EXOMEDEPTHCOV_wrapper.R
+++ b/R/tests/test_EXOMEDEPTHCOV_wrapper.R
@@ -4,6 +4,8 @@ library(CNVCALLER.RUNNER)
 context("Testing run_wrapper_EXOMEDEPTHCOV function")
 
 test_that("basic test for run_wrapper_EXOMEDEPTHCOV function, without calls",{
+  reference_set_select_method <- "exomedepth"
+  num_of_samples_in_reference_set <- "0"
   cov_table <- as.data.frame(matrix(data=c("NA12044",6,"2",17211,173310,151.636363636364,
                                            "NA11829",4,"2",16712,167190,100,
                                            "NA07051",7,"2",30275,304310,155.929936305732,
@@ -61,7 +63,9 @@ test_that("basic test for run_wrapper_EXOMEDEPTHCOV function, without calls",{
   cov_table[,"pos_min"] <- as.integer(as.character(cov_table[,"pos_min"]))
   cov_table[,"pos_max"] <- as.integer(as.character(cov_table[,"pos_max"]))
   cov_table[,"read_count"] <- as.numeric(as.character(cov_table[,"read_count"]))
-  calls <- run_wrapper_EXOMEDEPTHCOV(cov_table)
+  calls <- run_wrapper_EXOMEDEPTHCOV(reference_set_select_method,
+                                     num_of_samples_in_reference_set,
+                                     cov_table)
   expect_equal(length(calls), 0)
 })
 
diff --git a/R/tests/test_functions_REFERENCE_SAMPLE_SET_SELECTOR.R b/R/tests/test_functions_REFERENCE_SAMPLE_SET_SELECTOR.R
new file mode 100644
index 0000000..424d923
--- /dev/null
+++ b/R/tests/test_functions_REFERENCE_SAMPLE_SET_SELECTOR.R
@@ -0,0 +1,111 @@
+library(testthat)
+library(REFERENCE.SAMPLE.SET.SELECTOR)
+
+context("Testing canoes_method function")
+
+test_that("basic test for canoes_method function with num_refs equal to 5",{
+  investigated_sample <- "sample_1"
+  num_refs <- 5
+  Y = matrix(as.integer(c(52, 53, 54, 56, 70, 147, 58, 51,
+                          25, 22, 39, 31, 27, 217, 17, 27,
+                          57, 61, 66, 65, 64, 177, 67, 59, 
+                          101, 100, 90, 85, 98, 182, 92, 103, 
+                          149, 154, 140, 160, 165, 213, 149, 162, 
+                          95, 90, 90, 98, 88, 262, 93, 103, 
+                          50, 55, 45, 52, 59, 135, 65, 55, 
+                          79, 57, 63, 90, 71, 216, 69, 66)), nrow=8, ncol=8, byrow = TRUE)
+  colnames(Y) <- c("sample_1", "sample_2", "sample_3", "sample_4", "sample_5", "sample_6", "sample_7", "sample_8")
+  reference_samples <- canoes_method(investigated_sample, 
+                                     Y,
+                                     num_refs)
+  expect_equal(length(reference_samples), 1)
+  expect_equal(length(reference_samples$reference_samples), num_refs)
+  expect_equal(reference_samples$reference_samples[1], 'sample_8')
+  expect_equal(reference_samples$reference_samples[2], 'sample_4')
+  expect_equal(reference_samples$reference_samples[3], 'sample_2')
+  expect_equal(reference_samples$reference_samples[4], 'sample_3')
+  expect_equal(reference_samples$reference_samples[5], 'sample_7')
+})
+
+test_that("basic test for canoes_method function with num_refs equal to 0",{
+  investigated_sample <- "sample_1"
+  num_refs <- 0
+  Y = matrix(as.integer(c(52, 53, 54, 56, 70, 147, 58, 51,
+                          25, 22, 39, 31, 27, 217, 17, 27,
+                          57, 61, 66, 65, 64, 177, 67, 59, 
+                          101, 100, 90, 85, 98, 182, 92, 103, 
+                          149, 154, 140, 160, 165, 213, 149, 162, 
+                          95, 90, 90, 98, 88, 262, 93, 103, 
+                          50, 55, 45, 52, 59, 135, 65, 55, 
+                          79, 57, 63, 90, 71, 216, 69, 66)), nrow=8, ncol=8, byrow = TRUE)
+  colnames(Y) <- c("sample_1", "sample_2", "sample_3", "sample_4", "sample_5", "sample_6", "sample_7", "sample_8")
+  reference_samples <- canoes_method(investigated_sample, 
+                                     Y,
+                                     num_refs)
+  expect_equal(length(reference_samples), 1)
+  expect_equal(length(reference_samples$reference_samples), 7)
+  expect_equal(reference_samples$reference_samples[1], 'sample_8')
+  expect_equal(reference_samples$reference_samples[2], 'sample_4')
+  expect_equal(reference_samples$reference_samples[3], 'sample_2')
+  expect_equal(reference_samples$reference_samples[4], 'sample_3')
+  expect_equal(reference_samples$reference_samples[5], 'sample_7')
+  expect_equal(reference_samples$reference_samples[6], 'sample_5')
+  expect_equal(reference_samples$reference_samples[7], 'sample_6')
+})
+
+
+context("Testing exomedepth_method function")
+
+test_that("basic test for exomedepth_method function with num_refs equal to 5",{
+  investigated_sample <- "sample_1"
+  num_refs <- 5
+  target_length <- c(1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000)
+  Y = matrix(as.integer(c(52, 53, 54, 56, 70, 147, 58, 51,
+                          25, 22, 39, 31, 27, 217, 17, 27,
+                          57, 61, 66, 65, 64, 177, 67, 59, 
+                          101, 100, 90, 85, 98, 182, 92, 103, 
+                          149, 154, 140, 160, 165, 213, 149, 162, 
+                          95, 90, 90, 98, 88, 262, 93, 103, 
+                          50, 55, 45, 52, 59, 135, 65, 55, 
+                          79, 57, 63, 90, 71, 216, 69, 66)), nrow=8, ncol=8, byrow = TRUE)
+  colnames(Y) <- c("sample_1", "sample_2", "sample_3", "sample_4", "sample_5", "sample_6", "sample_7", "sample_8")
+  reference_samples <- exomedepth_method(investigated_sample, 
+                                         Y,
+                                         num_refs,
+                                         target_length)
+  expect_equal(length(reference_samples), 1)
+  expect_equal(length(reference_samples$reference_samples), num_refs)
+  expect_equal(reference_samples$reference_samples[1], 'sample_8')
+  expect_equal(reference_samples$reference_samples[2], 'sample_4')
+  expect_equal(reference_samples$reference_samples[3], 'sample_5')
+  expect_equal(reference_samples$reference_samples[4], 'sample_2')
+  expect_equal(reference_samples$reference_samples[5], 'sample_3')
+})
+
+test_that("basic test for exomedepth_method function with num_refs equal to 0 (auto)",{
+  investigated_sample <- "sample_1"
+  num_refs <- 0
+  target_length <- c(1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000)
+  Y = matrix(as.integer(c(52, 53, 54, 56, 70, 147, 58, 51,
+                          25, 22, 39, 31, 27, 217, 17, 27,
+                          57, 61, 66, 65, 64, 177, 67, 59, 
+                          101, 100, 90, 85, 98, 182, 92, 103, 
+                          149, 154, 140, 160, 165, 213, 149, 162, 
+                          95, 90, 90, 98, 88, 262, 93, 103, 
+                          50, 55, 45, 52, 59, 135, 65, 55, 
+                          79, 57, 63, 90, 71, 216, 69, 66)), nrow=8, ncol=8, byrow = TRUE)
+  colnames(Y) <- c("sample_1", "sample_2", "sample_3", "sample_4", "sample_5", "sample_6", "sample_7", "sample_8")
+  reference_samples <- exomedepth_method(investigated_sample, 
+                                         Y,
+                                         num_refs,
+                                         target_length)
+  expect_equal(length(reference_samples), 1)
+  expect_equal(length(reference_samples$reference_samples), 6)
+  expect_equal(reference_samples$reference_samples[1], 'sample_8')
+  expect_equal(reference_samples$reference_samples[2], 'sample_4')
+  expect_equal(reference_samples$reference_samples[3], 'sample_5')
+  expect_equal(reference_samples$reference_samples[4], 'sample_2')
+  expect_equal(reference_samples$reference_samples[5], 'sample_3')
+  expect_equal(reference_samples$reference_samples[6], 'sample_7')
+})
+
diff --git a/R/tests/test_run_cnvcaller.R b/R/tests/test_run_cnvcaller.R
index 49e5b73..0d3e8fb 100644
--- a/R/tests/test_run_cnvcaller.R
+++ b/R/tests/test_run_cnvcaller.R
@@ -31,7 +31,7 @@ test_that("basic test for reading test parameters",{
   parameters <- dbGetQuery(conn_psql, "select * from test_parameters")
   expect_equal(nrow(parameters), 20)
   read_parameters_result <- read_parameters("test_parameters", parameters[1,"id"], conn_psql)
-  expect_equal(length(read_parameters_result), 14)
+  expect_equal(length(read_parameters_result), 16)
   expect_equal(read_parameters_result$caller, parameters[1,"caller"])
   expect_equal(read_parameters_result$cov_table, parameters[1,"cov_table"])
   expect_equal(read_parameters_result$mapp_thresh, parameters[1,"mapp_thresh"])
@@ -45,6 +45,8 @@ test_that("basic test for reading test parameters",{
   expect_equal(read_parameters_result$K_to, parameters[1,"k_to"])
   expect_equal(read_parameters_result$lmax, parameters[1,"lmax"])
   expect_equal(read_parameters_result$chr, parameters[1,"chr"])
+  expect_equal(read_parameters_result$reference_set_select_method, parameters[1,"reference_set_select_method"])
+  expect_equal(read_parameters_result$num_of_samples_in_reference_set, parameters[1,"num_of_samples_in_reference_set"])
 })
 
 context("Testing read_coverage_table function")
@@ -82,7 +84,9 @@ test_that("basic test for run_caller function, without calls",{
                      gc_thresh_to="100",
                      K_from="1",
                      K_to="3",
-                     lmax="200")
+                     lmax="200",
+                     reference_set_select_method="codex",
+                     num_of_samples_in_reference_set="0")
   cov_table <- as.data.frame(matrix(data=c("NA12044",6,"2",17211,173310,151.636363636364,
                                            "NA11829",4,"2",16712,167190,100,
                                            "NA07051",7,"2",30275,304310,155.929936305732,
diff --git a/sql/parameter_generator.sql b/sql/parameter_generator.sql
index 7815785..9e9dc96 100644
--- a/sql/parameter_generator.sql
+++ b/sql/parameter_generator.sql
@@ -17,7 +17,9 @@ select cast(row_number() over () as integer ) as id,
   k_from ,
   k_to ,
   lmax,
-  chr
+  chr,
+  reference_set_select_method,
+  num_of_samples_in_reference_set
 from test_params_bck par, (select cast (generate_series(1,22) as text ) as chr
 union
 select 'X'
diff --git a/sql/parameters_calls_model.sql b/sql/parameters_calls_model.sql
index 74f0b62..6fe4d45 100644
--- a/sql/parameters_calls_model.sql
+++ b/sql/parameters_calls_model.sql
@@ -15,7 +15,9 @@ CREATE TABLE IF NOT EXISTS CNV.PARAMETERS (
     gc_thresh_to TEXT,
     K_from TEXT,
     K_to TEXT,
-    lmax TEXT
+    lmax TEXT,
+    reference_set_select_method TEXT,
+    num_of_samples_in_reference_set TEXT
 );
 
 CREATE TABLE IF NOT EXISTS CNV.CALLS (
diff --git a/sql/parameters_calls_test_model.sql b/sql/parameters_calls_test_model.sql
index 82c2551..6fff3b4 100644
--- a/sql/parameters_calls_test_model.sql
+++ b/sql/parameters_calls_test_model.sql
@@ -11,39 +11,41 @@ CREATE TABLE IF NOT EXISTS TEST_PARAMETERS (
     gc_thresh_to TEXT,
     K_from TEXT,
     K_to TEXT,
-    lmax TEXT
+    lmax TEXT,
+    reference_set_select_method TEXT,
+    num_of_samples_in_reference_set TEXT
 );
 
-INSERT INTO TEST_PARAMETERS VALUES (1, 'codex', 'cnv.coverage_target', 2, 3, 4, 5, 6, 7, 8, 9, 10, 11);
-INSERT INTO TEST_PARAMETERS VALUES (2, 'codex', 'cnv.coverage_target', 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
-INSERT INTO TEST_PARAMETERS VALUES (3, 'codex', 'cnv.coverage_target', 4, 5, 6, 7, 8, 9, 10, 11, 12, 13);
-INSERT INTO TEST_PARAMETERS VALUES (4, 'codex', 'cnv.coverage_target', 5, 6, 7, 8, 9, 10, 11, 12, 13, 14);
-INSERT INTO TEST_PARAMETERS VALUES (5, 'codex', 'cnv.coverage_target', 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
-INSERT INTO TEST_PARAMETERS VALUES (6, 'codex', 'cnv.coverage_target', 7, 8, 9, 10, 11, 12, 13, 14, 15, 16);
-INSERT INTO TEST_PARAMETERS VALUES (7, 'codex', 'cnv.coverage_target', 8, 9, 10, 11, 12, 13, 14, 15, 16, 17);
-INSERT INTO TEST_PARAMETERS VALUES (8, 'codex', 'cnv.coverage_target', 9, 10, 11, 12, 13, 14, 15, 16, 17, 18);
-INSERT INTO TEST_PARAMETERS VALUES (9, 'codex', 'cnv.coverage_target', 10, 11, 12, 13, 14, 15, 16, 17, 18, 19);
-INSERT INTO TEST_PARAMETERS VALUES (10, 'xhmm', 'cnv.coverage_target', 11, 12, 13, 14, 15, 16, 17, 18, 19, 20);
-INSERT INTO TEST_PARAMETERS VALUES (11, 'xhmm', 'cnv.coverage_target', 12, 13, 14, 15, 16, 17, 18, 19, 20, 21);
-INSERT INTO TEST_PARAMETERS VALUES (12, 'xhmm', 'cnv.coverage_target', 13, 14, 15, 16, 17, 18, 19, 20, 21, 22);
-INSERT INTO TEST_PARAMETERS VALUES (13, 'xhmm', 'cnv.coverage_target', 14, 15, 16, 17, 18, 19, 20, 21, 22, 23);
-INSERT INTO TEST_PARAMETERS VALUES (14, 'xhmm', 'cnv.coverage_target', 15, 16, 17, 18, 19, 20, 21, 22, 23, 24);
-INSERT INTO TEST_PARAMETERS VALUES (15, 'xhmm', 'cnv.coverage_target', 16, 17, 18, 19, 20, 21, 22, 23, 24, 25);
-INSERT INTO TEST_PARAMETERS VALUES (16, 'xhmm', 'cnv.coverage_target', 17, 18, 19, 20, 21, 22, 23, 24, 25, 26);
-INSERT INTO TEST_PARAMETERS VALUES (17, 'xhmm', 'cnv.coverage_target', 18, 19, 20, 21, 22, 23, 24, 25, 26, 27);
-INSERT INTO TEST_PARAMETERS VALUES (18, 'xhmm', 'cnv.coverage_target', 19, 20, 21, 22, 23, 24, 25, 26, 27, 28);
-INSERT INTO TEST_PARAMETERS VALUES (19, 'xhmm', 'cnv.coverage_target', 20, 21, 22, 23, 24, 25, 26, 27, 28, 29);
-INSERT INTO TEST_PARAMETERS VALUES (20, 'xhmm', 'cnv.coverage_target', 21, 22, 23, 24, 25, 26, 27, 28, 29, 30);
-INSERT INTO TEST_PARAMETERS VALUES (21, 'exomedepth', 'cnv.coverage_target', 22, 23, 24, 25, 26, 27, 28, 29, 30, 31);
-INSERT INTO TEST_PARAMETERS VALUES (22, 'exomedepth', 'cnv.coverage_target', 23, 24, 25, 26, 27, 28, 29, 30, 31, 32);
-INSERT INTO TEST_PARAMETERS VALUES (23, 'exomedepth', 'cnv.coverage_target', 24, 25, 26, 27, 28, 29, 30, 31, 32, 33);
-INSERT INTO TEST_PARAMETERS VALUES (24, 'exomedepth', 'cnv.coverage_target', 25, 26, 27, 28, 29, 30, 31, 32, 33, 34);
-INSERT INTO TEST_PARAMETERS VALUES (25, 'exomedepth', 'cnv.coverage_target', 26, 27, 28, 29, 30, 31, 32, 33, 34, 35);
-INSERT INTO TEST_PARAMETERS VALUES (26, 'exomedepth', 'cnv.coverage_target', 27, 28, 29, 30, 31, 32, 33, 34, 35, 36);
-INSERT INTO TEST_PARAMETERS VALUES (27, 'exomedepth', 'cnv.coverage_target', 28, 29, 30, 31, 32, 33, 34, 35, 36, 37);
-INSERT INTO TEST_PARAMETERS VALUES (28, 'exomedepth', 'cnv.coverage_target', 29, 30, 31, 32, 33, 34, 35, 36, 37, 38);
-INSERT INTO TEST_PARAMETERS VALUES (29, 'exomedepth', 'cnv.coverage_target', 30, 31, 32, 33, 34, 35, 36, 37, 38, 39);
-INSERT INTO TEST_PARAMETERS VALUES (30, 'exomedepth', 'cnv.coverage_target', 31, 32, 33, 34, 35, 36, 37, 38, 39, 40);
+INSERT INTO TEST_PARAMETERS VALUES (1, 'codex', 'cnv.coverage_target', 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (2, 'codex', 'cnv.coverage_target', 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (3, 'codex', 'cnv.coverage_target', 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (4, 'codex', 'cnv.coverage_target', 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (5, 'codex', 'cnv.coverage_target', 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (6, 'codex', 'cnv.coverage_target', 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (7, 'codex', 'cnv.coverage_target', 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (8, 'codex', 'cnv.coverage_target', 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (9, 'codex', 'cnv.coverage_target', 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (10, 'xhmm', 'cnv.coverage_target', 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (11, 'xhmm', 'cnv.coverage_target', 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (12, 'xhmm', 'cnv.coverage_target', 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (13, 'xhmm', 'cnv.coverage_target', 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (14, 'xhmm', 'cnv.coverage_target', 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (15, 'xhmm', 'cnv.coverage_target', 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (16, 'xhmm', 'cnv.coverage_target', 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (17, 'xhmm', 'cnv.coverage_target', 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (18, 'xhmm', 'cnv.coverage_target', 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (19, 'xhmm', 'cnv.coverage_target', 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (20, 'xhmm', 'cnv.coverage_target', 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (21, 'exomedepth', 'cnv.coverage_target', 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (22, 'exomedepth', 'cnv.coverage_target', 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (23, 'exomedepth', 'cnv.coverage_target', 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (24, 'exomedepth', 'cnv.coverage_target', 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (25, 'exomedepth', 'cnv.coverage_target', 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (26, 'exomedepth', 'cnv.coverage_target', 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (27, 'exomedepth', 'cnv.coverage_target', 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (28, 'exomedepth', 'cnv.coverage_target', 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (29, 'exomedepth', 'cnv.coverage_target', 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 'exomedepth', 5);
+INSERT INTO TEST_PARAMETERS VALUES (30, 'exomedepth', 'cnv.coverage_target', 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 'exomedepth', 5);
 
 CREATE TABLE IF NOT EXISTS TEST_CALLS (
     id SERIAL PRIMARY KEY,