diff --git a/R/read.modkit.R b/R/read.modkit.R
index d4f6af2..b79158d 100644
--- a/R/read.modkit.R
+++ b/R/read.modkit.R
@@ -9,13 +9,15 @@ read.modkit <- function(files,
         }
 
     for (i in seq_along(files)){
-        data <- read.table(files[i], header = FALSE, sep="\t",
+        data <- fread(files[i], header = FALSE, sep="\t",
                     stringsAsFactors=FALSE, quote="")
+        data$V6[data$V6 == "."] <- "*"
 
         if (length(unique(data$V4)) == 2){
             gr <- GRanges(seqnames = data[data$V4 == "m", ]$V1,
                     ranges = IRanges(start = data[data$V4 == "m", ]$V2+1,
-                                     end = data[data$V4 == "m", ]$V3))
+                                     end = data[data$V4 == "m", ]$V3),
+                    strand = data[data$V4 == "m", ]$V6)
 
             mcols(gr)$m <- data[data$V4 == "m", ]$V12
             mcols(gr)$h <- data[data$V4 != "m", ]$V12
@@ -23,7 +25,8 @@ read.modkit <- function(files,
             mcols(gr)$filter <- data[data$V4 == "m", ]$V16
         }else{
             gr <- GRanges(seqnames = data$V1,
-                          ranges = IRanges(start = data$V2+1, end = data$V3))
+                          ranges = IRanges(start = data$V2+1, end = data$V3),
+                          strand = data$V6)
 
             mcols(gr)$m <- data$V12
             mcols(gr)$u <- data$V13
@@ -58,19 +61,23 @@ read.modkit <- function(files,
                              Filtered = as.matrix(filter),
                              coef = NULL, se.coef = NULL,
                              pos = start(overlap_gr), trans = NULL,
-                             parameters = NULL, pData = colData, gr = NULL,
+                             parameters = NULL, pData = colData, gr = overlap_gr,
                              chr = as.vector(seqnames(overlap_gr)),
                              sampleNames = sampleNames, rmZeroCov = rmZeroCov)
-          if (strandCollapse) {strandCollapse(bsseq_obj)}
+          if (strandCollapse) {
+            bsseq_obj <- strandCollapse(bsseq_obj)
+          }
     }else{
       bsseq_obj <- BSseq(M = as.matrix(m), Cov = as.matrix(u + m),
                          Filtered = as.matrix(filter),
                          coef = NULL, se.coef = NULL,
                          pos = start(overlap_gr), trans = NULL,
-                         parameters = NULL, pData = colData, gr = NULL,
+                         parameters = NULL, pData = colData, gr = overlap_gr,
                          chr = as.vector(seqnames(overlap_gr)),
                          sampleNames = sampleNames, rmZeroCov = rmZeroCov)
-      if (strandCollapse) {strandCollapse(bsseq_obj)}
+      if (strandCollapse) {
+        bsseq_obj <- strandCollapse(bsseq_obj)
+      }
     }
 
     return(bsseq_obj)