diff --git a/R/grn.R b/R/grn.R
index b4a30ee..cf3c891 100644
--- a/R/grn.R
+++ b/R/grn.R
@@ -201,7 +201,7 @@ fit_grn_models.GRNData <- function(
             group_name = aggregate_rna_col,
             verbose = FALSE
         )
-        gene_groups <- object@meta.data[[aggregate_rna_col]]
+        gene_groups <- object@data@meta.data[[aggregate_rna_col]]
     }
 
     if (is.null(aggregate_peaks_col)){
@@ -285,7 +285,7 @@ fit_grn_models.GRNData <- function(
             log_message('Warning: No peaks found near ', g, verbose=verbose==2)
             return()
         }
-
+        log_message('Selecting peaks correlating with target gene expression...', verbose=verbose)
         # Select peaks correlating with target gene expression
         g_x <- gene_data[gene_groups, g, drop=FALSE]
         peak_x <- peak_data[peak_groups, gene_peaks, drop=FALSE]
@@ -300,6 +300,7 @@ fit_grn_models.GRNData <- function(
         peak_motifs <- peaks2motif[gene_peaks, , drop=FALSE][peaks_use, , drop=FALSE]
 
         # Select TFs with motifs in peaks
+        log_message('Selecting TFs with motifs in peaks...', verbose=verbose)
         gene_peak_tfs <- map(rownames(peak_motifs), function(p){
             x <- as.logical(peak_motifs[p, ])
             peak_tfs <- colMaxs(motif2tf[x, , drop=FALSE])
@@ -310,6 +311,7 @@ fit_grn_models.GRNData <- function(
         names(gene_peak_tfs) <- rownames(peak_motifs)
 
         # Check correlation of peaks with target gene
+        log_message('Check correlation of peaks with target gene...', verbose=verbose)
         gene_tfs <- purrr::reduce(gene_peak_tfs, union)
         tf_x <- gene_data[gene_groups, gene_tfs, drop=FALSE]
         tf_g_cor <- as(sparse_cor(tf_x, g_x), 'generalMatrix')
@@ -334,7 +336,9 @@ fit_grn_models.GRNData <- function(
                 return()
             }
             peak_name <- str_replace_all(p, '-', '_')
+            peak_name <- paste0("`", peak_name, "`")
             tf_name <- str_replace_all(peak_tfs, '-', '_')
+            tf_name <- paste0("`", tf_name, "`")
             formula_str <- paste(
                 paste(peak_name, interaction_term, tf_name, sep=' '), collapse = ' + ')
             return(list(tfs=peak_tfs, frml=formula_str))
@@ -345,19 +349,23 @@ fit_grn_models.GRNData <- function(
             return()
         }
 
+        
         target <- str_replace_all(g, '-', '_')
+        target <- paste0("`", target, "`")
         model_frml <- as.formula(
-            paste0(target, ' ~ ', paste0(map(frml_string, function(x) x$frml),  collapse=' + '))
+            paste0(target, ' ~ ', paste0(map(frml_string, function(x) x$`frml`),  collapse=' + '))
         )
-
+        
         # Get expression data
-        nfeats <- sum(map_dbl(frml_string, function(x) length(x$tfs)))
-        gene_tfs <- purrr::reduce(map(frml_string, function(x) x$tfs), union)
+                                        
+        nfeats <- sum(map_dbl(frml_string, function(x) length(x$`tfs`)))               
+        gene_tfs <- purrr::reduce(map(frml_string, function(x) x$`tfs`), union)
         gene_x <- gene_data[gene_groups, union(g, gene_tfs), drop=FALSE]
         model_mat <- as.data.frame(cbind(gene_x, peak_x))
         if (scale) model_mat <- as.data.frame(scale(as.matrix(model_mat)))
+                                      
         colnames(model_mat) <- str_replace_all(colnames(model_mat), '-', '_')
-
+                                      
         log_message('Fitting model with ', nfeats, ' variables for ', g, verbose=verbose==2)
         result <- try(fit_model(
             model_frml,
@@ -381,13 +389,13 @@ fit_grn_models.GRNData <- function(
         log_message('Warning: Fitting model failed for all genes.', verbose=verbose)
     }
 
-    coefs <- map_dfr(model_fits, function(x) x$coefs, .id='target')
+    coefs <- map_dfr(model_fits, function(x) x$`coefs`, .id='target')
     coefs <- format_coefs(coefs, term=interaction_term, adjust_method=adjust_method)
-    corrs <- map_dfr(model_fits, function(x) x$corr, .id='target')
+    corrs <- map_dfr(model_fits, function(x) x$`corr`, .id='target')
     if (nrow(coefs)>0){
         coefs <- suppressMessages(left_join(coefs, corrs))
     }
-    gof <- map_dfr(model_fits, function(x) x$gof, .id='target')
+    gof <- map_dfr(model_fits, function(x) x$`gof`, .id='target')
 
     params <- list()
     params[['method']] <- method
@@ -445,15 +453,15 @@ format_coefs <- function(coefs, term=':', adjust_method='fdr'){
         ) %>%
         select(-region_, -tf_) %>%
         mutate(
-            region = str_replace_all(region, '_', '-'),
-            tf = str_replace_all(tf, '_', '-'),
-            target = str_replace_all(target, '_', '-')
+            term = str_replace_all(term, '_', '-') %>% str_replace_all('`', ''),
+            region = str_replace_all(region, '_', '-') %>% str_replace_all('`', ''),
+            tf = str_replace_all(tf, '_', '-') %>% str_replace_all('`', ''),
+            target = str_replace_all(target, '_', '-') %>% str_replace_all('`', '')
         ) %>%
         select(tf, target, region, term, everything())
     return(coefs_use)
 }
-
-
+                   
 
 #' Find TF modules in regulatory network
 #'