diff --git a/NEWS.md b/NEWS.md
index 078fc677..e9e12070 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 ## QFeatures 1.7.3
 
-- Nothing yet.
+- fix: fixed filterFeatures when selection contains environment variables
 
 ## QFeatures 1.7.2
 
diff --git a/R/QFeatures-filter.R b/R/QFeatures-filter.R
index 20eb09e5..c405582b 100644
--- a/R/QFeatures-filter.R
+++ b/R/QFeatures-filter.R
@@ -309,6 +309,15 @@ filterFeaturesWithFormula <- function(object, filter, na.rm, keep, ...) {
 
 ## Internal function that 
 .checkFilterVariables <- function(object, vars) {
+    ## Ignore variables from the user environment. We search for 
+    ## variable in the 4th parent environment (may not always be 
+    ## .GlobalEnv). Here is a "traceback" counter:
+    ## 0 in .checkFilterVariables()
+    ## 1 in FilterFeaturesWithFormula()
+    ## 2 in .local()
+    ## 3 in filterFeatures()
+    ## 4 in environment the function was called
+    vars <- vars[!vars %in% ls(envir = parent.frame(4))]
     ## Get in which assays each variable comes from
     out <- sapply(rowDataNames(object), function(rdn) vars %in% rdn)
     if (!is.array(out)) out <- t(out)
diff --git a/R/QFeatures-join.R b/R/QFeatures-join.R
index 4a38e639..d0d817cc 100644
--- a/R/QFeatures-join.R
+++ b/R/QFeatures-join.R
@@ -25,22 +25,24 @@
     res
 }
 
-##' @importFrom methods as
-.merge_2_by_rows <- function(x, y) {
-    ## Save class to coerce at the end
-    cl <- class(x)
-    res <- merge(x, y,
-                 by = 0,
-                 all.x = TRUE, all.y = TRUE,
-                 sort = FALSE)
-    ## Set and remove row names
-    rownames(res) <- res[[1]]
-    res <- res[, -1, drop = FALSE]
-    as(res, cl[1])
-}
 
-.merge_by_rows <- function(x, y, ...) {
-    Reduce(.merge_2_by_rows, list(x, y, ...))
+.merge_assays_by_rows <- function(l) {
+    cn <- unlist(lapply(l, colnames))
+    rn <- unique(unlist(lapply(l, rownames)))
+    
+    ## Check for duplicate column (sample) names
+    if (any(duplicated(cn))) 
+        stop("Merging assays with columns in common is not allowed.")
+    
+    res <- matrix(NA, ncol = length(cn), nrow = length(rn), 
+                  dimnames = list(rn, cn))
+    for (i in seq_along(l)) {
+        x <- l[[i]]
+        res[rownames(x), colnames(x)] <- as.matrix(x) 
+        ## as.matrix in case x is an HDF5Array, note x (and res) are
+        ## realized in memory. 
+    }
+    res
 }
 
 
@@ -54,7 +56,7 @@ mergeSElist <- function(x) {
     if (length(x_classes) != 1)
         stop("Can't join assays from different classes.", call. = FALSE)
     joined_mcols <- Reduce(.merge_2_by_cols, lapply(x, rowData))
-    joined_assay <- Reduce(.merge_2_by_rows, lapply(x, assay))
+    joined_assay <- .merge_assays_by_rows(lapply(x, assay))
     joined_coldata <- Reduce(.merge_2_by_cols, lapply(x, colData))
     res <- SummarizedExperiment(joined_assay[rownames(joined_mcols), ],
                                 joined_mcols,
@@ -152,10 +154,9 @@ joinAssays <- function(x,
               "Need at least 2 assays to join" = length(i) >= 2)
     if (name %in% names(x))
         stop("Assay with name '", name, "' already exists.")
-    joined_se <- mergeSElist(as.list(experiments(x)[i]))
+    ## Join assays and add to x
+    joined_se <- mergeSElist(experiments(x)[i])
     x <- addAssay(x, joined_se, name = name)
     ## Add the multi-parent AssayLinks
-    if (is.numeric(i)) i <- names(x)[i]
-    al <- .create_assay_link(x, from = i, to = name)
-    .update_assay_links(x, al)
+    addAssayLink(x, from = i, to = name)
 }
diff --git a/R/QFeatures-missing-data.R b/R/QFeatures-missing-data.R
index eb176b01..be490f1c 100644
--- a/R/QFeatures-missing-data.R
+++ b/R/QFeatures-missing-data.R
@@ -174,9 +174,7 @@ setMethod("zeroIsNA", c("QFeatures", "integer"),
               el <- experiments(object)
               for (ii in i)
                   el[[ii]] <- zeroIsNA(el[[ii]])
-              BiocGenerics:::replaceSlots(object,
-                                          ExperimentList = el,
-                                          check = TRUE)
+              replaceAssay(object, el)
           })
 
 ##' @rdname QFeatures-missing-data
@@ -206,9 +204,7 @@ setMethod("infIsNA", c("QFeatures", "integer"),
               el <- experiments(object)
               for (ii in i)
                   el[[ii]] <- infIsNA(el[[ii]])
-              BiocGenerics:::replaceSlots(object,
-                                          ExperimentList = el,
-                                          check = TRUE)
+              replaceAssay(object, el)
           })
 
 ##' @rdname QFeatures-missing-data
diff --git a/tests/testthat/test_filterFeatures.R b/tests/testthat/test_filterFeatures.R
index 46926216..383e7eb6 100644
--- a/tests/testthat/test_filterFeatures.R
+++ b/tests/testthat/test_filterFeatures.R
@@ -20,6 +20,19 @@ test_that("filterFeatures", {
     expect_equal(filter1, filter7)
     expect_identical(lengths(filter1), c(6L, 2L, 1L))
     
+    ## Test filter stored in variable
+    target <- "Mitochondrion"
+    filter8 <- expect_message(filterFeatures(feat1, ~  location == target))
+    filter9 <- expect_message(filterFeatures(feat1, VariableFilter("location", target)))
+    expect_equal(filter1, filter8)
+    expect_equal(filter8, filter9)
+    ## Test filter stored in variable within function
+    runfilter <- function() {
+        target2 <- "Mitochondrion"
+        expect_message(filterFeatures(feat1, ~  location == target2))
+    }
+    expect_equal(filter8, runfilter())
+    
     ## Test numerical filters
     filter1 <- expect_message(filterFeatures(feat1, VariableFilter("pval", 0.03, "<=")))
     filter2 <- expect_message(filterFeatures(feat1, ~ pval <= 0.03))