diff --git a/DESCRIPTION b/DESCRIPTION
index f20b46f..d436b65 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: AlphaBeta
 Type: Package
 Title: Computational inference of epimutation rates and spectra from high-throughput DNA methylation data in plants
-Version: 1.3.2
+Version: 1.3.3
 Authors@R: c(
     person("Yadollah", "Shahryary Dizaji", email = "shahryary@gmail.com", role = c("cre","aut")),
     person("Frank", "Johannes", email = "frank@jlab.org", role = "aut"),
diff --git a/NEWS b/NEWS
index 8251d3d..9c1c228 100644
--- a/NEWS
+++ b/NEWS
@@ -23,3 +23,7 @@ Changes in version 1.2.1 (2020-04-20)
     + Added some plots.
 * Added SOMA functions.
 * Plotting estimates
+
+Changes in version 1.2.3 (2020-06-06)
+* Fixed bugs.
+* Compatible with new version of dplyr,data.tables
diff --git a/R/BOOTmodel.R b/R/BOOTmodel.R
index 5675ef7..0d2445d 100644
--- a/R/BOOTmodel.R
+++ b/R/BOOTmodel.R
@@ -36,15 +36,17 @@ BOOTmodel<-function(pedigree.data, Nboot, out.dir, out.name)
   {
 
       ## Reading the dataset for bootstrapping and extracting the parameter settings
-    	settings<-pedigree.data$settings
-    	est<-pedigree.data$estimates
-    	eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
-    	eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
-    	optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-    	p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-    	p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-    	p0um<-1-(p0mm + p0uu)
-    	pedigree<-pedigree.data$pedigree
+      settings<-pedigree.data$settings
+      est<-pedigree.data$estimates
+      eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
+      eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
+      optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
+      #p0mm<-round(as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2])),16)
+      p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
+      #p0um<-round(as.numeric(as.character(settings[which(settings[,1] == "p0um"),2])),15)
+      pedigree<-pedigree.data$pedigree
+      p0mm = 1-p0uu
+      p0um = 0
 
 
     	if(sum(c(p0mm, p0um, p0uu), na.rm =TRUE) != 1) {stop("The initial state probabilities don't sum to 1")}
@@ -258,10 +260,10 @@ BOOTmodel<-function(pedigree.data, Nboot, out.dir, out.name)
     eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
     eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
     optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-    p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-    p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-    p0um<-1-(p0mm + p0uu)
+    p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
     pedigree<-pedigree.data$pedigree
+    p0mm = 1-p0uu
+    p0um = 0
 
 
     if(sum(c(p0mm, p0um, p0uu), na.rm =TRUE) != 1) {stop("The initial state probabilities don't sum to 1")}
@@ -477,11 +479,10 @@ BOOTmodel<-function(pedigree.data, Nboot, out.dir, out.name)
     eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
     eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
     optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-    p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-    p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-    p0um<-1-(p0mm + p0uu)
-    sel<-as.numeric(as.character(settings$Setting[6]))
+    p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
     pedigree<-pedigree.data$pedigree
+    p0mm = 1-p0uu
+    p0um = 0
 
 
     if(sum(c(p0mm, p0um, p0uu), na.rm =TRUE) != 1) {stop("The initial state probabilities don't sum to 1")}
@@ -700,11 +701,10 @@ BOOTmodel<-function(pedigree.data, Nboot, out.dir, out.name)
     eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
     eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
     optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-    p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-    p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-    p0um<-1-(p0mm + p0uu)
-    sel<-as.numeric(as.character(settings$Setting[6]))
+    p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
     pedigree<-pedigree.data$pedigree
+    p0mm = 1-p0uu
+    p0um = 0
 
 
     if(sum(c(p0mm, p0um, p0uu), na.rm =TRUE) != 1) {stop("The initial state probabilities don't sum to 1")}
@@ -919,15 +919,15 @@ BOOTmodel<-function(pedigree.data, Nboot, out.dir, out.name)
   if (model == "ABneutral.R")
   {
   ## Reading the dataset for bootstrapping and extracting the parameter settings
-  settings<-pedigree.data$settings
-  est<-pedigree.data$estimates
-  eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
-  eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
-  optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-  p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-  p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-  p0um<-1-(p0mm + p0uu)
-  pedigree<-pedigree.data$pedigree
+    settings<-pedigree.data$settings
+    est<-pedigree.data$estimates
+    eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
+    eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
+    optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
+    p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
+    pedigree<-pedigree.data$pedigree
+    p0mm = 1-p0uu
+    p0um = 0
 
 
   if(sum(c(p0mm, p0um, p0uu), na.rm =TRUE) != 1) {stop("The initial state probabilities don't sum to 1")}
@@ -1102,15 +1102,15 @@ allow.neg.intercept="yes"
 
 
   ##### Reading the dataset for bootstrapping and extracting the parameter settings
-  settings<-pedigree.data$settings
-  est<-pedigree.data$estimates
-  eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
-  eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
-  optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
-  p0mm<-as.numeric(as.character(settings[which(settings[,1] == "p0mm"),2]))
-  p0uu<-as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2]))
-  p0um<-1-(p0mm + p0uu)
-  pedigree<-pedigree.data$pedigree
+    settings<-pedigree.data$settings
+    est<-pedigree.data$estimates
+    eqp<-as.numeric(as.character(settings[which(settings[,1] == "eqp"),2]))
+    eqp.weight<-as.numeric(as.character(settings[which(settings[,1] == "eqp.weight"),2]))
+    optim.method<-as.character(settings[which(settings[,1] == "optim.method"),2])
+    p0uu<-round(as.numeric(as.character(settings[which(settings[,1] == "p0uu"),2])),16)
+    pedigree<-pedigree.data$pedigree
+    p0mm = 1-p0uu
+    p0um = 0
   #props<-pedigree.data$prop.data.used
 
 
diff --git a/R/dataRead.R b/R/dataRead.R
index e3f2d02..a0be1df 100644
--- a/R/dataRead.R
+++ b/R/dataRead.R
@@ -46,7 +46,7 @@ DM.dataRead <- function(dFile, cytosine, posteriorMaxFilter) {
 
 
     drops <- c("context", "posteriorMax")
-    tmp <- tmp %>% select(!(drops))
+    tmp <- tmp %>% select(!all_of(drops))
     #tmp <- tmp[ !(names(tmp) %in% drops)]
     # filter out based on context & posteriorMax
 
@@ -83,7 +83,7 @@ RC.dataRead <- function(dFile, cytosine, posteriorMaxFilter) {
     # filter out based on context & posteriorMax
     drops <- c("context", "posteriorMax")
     #tmp <- tmp[,!(names(tmp) %in% drops)]
-    tmp <- tmp %>% select(!(drops))
+    tmp <- tmp %>% select(!all_of(drops))
     return(tmp)
 }
 
diff --git a/R/rcMethlvl.R b/R/rcMethlvl.R
index ac3ba2e..448d3ce 100644
--- a/R/rcMethlvl.R
+++ b/R/rcMethlvl.R
@@ -52,17 +52,18 @@ RCsaveResult<-function(list.rc,cytosine,posteriorMaxFilter){
       tmp_dmr <- data.frame(matrix(ncol = 3, nrow =1 ))
       x <- c("Sample_name", "context", "rc.meth.lvls")
       colnames(tmp_dmr) <- x
-      mainRC<-NULL
+      mainRC <- tmp_dmr
       for (nlist in seq_len(length(list.rc))){
-        tmp_dmr$Sample_name<-list.rc[[nlist]][[1]]
+        tmp_dmr$Sample_name<-list.rc[[nlist]][[1]]$node
         tmp_dmr$context<-cytosine
         tmp_dmr$rc.meth.lvls <-list.rc[[nlist]][[2]]
         mainRC<-rbind(mainRC,tmp_dmr)
+        tmp_dmr<-NULL
       }
-      mainRC<-mainRC[mixedorder(mainRC$Sample_name),]
-      #saveFile <- paste0(getwd(), "/", "AB-methProp-", cytosine, "-", posteriorMaxFilter, ".csv")
-      #fwrite(mainRC,file = saveFile ,quote = FALSE,sep = '\t',row.names = FALSE,col.names = TRUE)
-      #cat(paste0("Methylation proportions results saved in: ",saveFile,"\n"))
+      mainRC <- mainRC[!(rowSums(x = is.na(x = mainRC)) == ncol(x = mainRC)),]
+      rownames(mainRC) <- NULL
+      mainRC %>% mutate_if(is.factor, as.character) -> mainRC
+      mainRC <- mainRC[mixedorder(mainRC$Sample_name), ]
       return(mainRC)
 
 }