R/read.bsmooth.R

In bsseq: Analyze, manage and store bisulfite sequencing data

read.umtab2 <- function(dirs, sampleNames = NULL, rmZeroCov = FALSE,
                        readCycle = FALSE, keepFilt = FALSE,
                        pattern = NULL, keepU, keepM, verbose = TRUE) {
    filesPerDir <- lapply(dirs, function(xx) sort(list.files(xx, pattern = pattern, full.names = TRUE)))
    if(!all(sapply(filesPerDir, function(xx) all(basename(xx) == basename(filesPerDir[[1]])))))
        warning("'dirs' does not contain the exact same file names")
    allChrs <- as.list(seq_along(filesPerDir[[1]]))
    for(ii in seq_along(filesPerDir[[1]])) {
        chrRead <- read.umtab2.chr(sapply(filesPerDir, function(xx) xx[ii]),
                                   keepM = keepM, keepU = keepU, sampleNames = sampleNames,
                                   readCycle = readCycle, keepFilt = keepFilt,
                                   verbose = verbose)
        if(rmZeroCov && length(chrRead) != 0){
            wh <- which(rowSums2(chrRead$U + chrRead$M) > 0)
            chrRead$M <- chrRead$M[wh,]
            chrRead$U <- chrRead$U[wh,]
            if(readCycle) {
                chrRead$Mcy <- chrRead$Mcy[wh,]
                chrRead$Ucy <- chrRead$Ucy[wh,]
            }
            chrRead$chr <- chrRead$chr[wh]
            chrRead$pos <- chrRead$pos[wh]
            if(keepFilt) {
                chrRead$filt_cycle <- chrRead$filt_cycle[wh,]
                chrRead$filt_allele <- chrRead$filt_allele[wh,]
                chrRead$filt_mapq <- chrRead$filt_mapq[wh,]
                chrRead$filt_baseq <- chrRead$filt_baseq[wh,]
            }
        }
        allChrs[[ii]] <- chrRead
        print(gc())
    }
    chr <- do.call(c, lapply(allChrs, function(xx) xx$chr))
    pos <- do.call(c, lapply(allChrs, function(xx) xx$pos))
    M <- do.call(rbind, lapply(allChrs, function(xx) xx$M))
    Cov <- do.call(rbind, lapply(allChrs, function(xx) xx$M)) +
        do.call(rbind, lapply(allChrs, function(xx) xx$U))
    if(readCycle) {
        Mcy <- do.call(rbind, lapply(allChrs, function(xx) xx$Mcy))
        Ucy <- do.call(rbind, lapply(allChrs, function(xx) xx$Ucy))
    } else {
        Mcy <- NULL
        Ucy <- NULL
    }
    if(keepFilt) {
        filt_cycle <- do.call(rbind, lapply(allChrs, function(xx) xx$filt_cycle))
        filt_allele <- do.call(rbind, lapply(allChrs, function(xx) xx$filt_allele))
        filt_mapq <- do.call(rbind, lapply(allChrs, function(xx) xx$filt_mapq))
        filt_baseq <- do.call(rbind, lapply(allChrs, function(xx) xx$filt_baseq))
    } else {
        filt_cycle <- filt_allele <- filt_mapq <- filt_baseq <- NULL
    }
    print(gc())
    csums <- lapply(allChrs, function(xx) xx$csums)
    names(csums) <- basename(filesPerDir[[1]])
    rm(allChrs)
    BSdata <- BSseq(chr = chr, pos = pos, M = M, Cov = Cov, rmZeroCov = rmZeroCov)
    print(gc())
    return(list(BSdata = BSdata, filt_cycle = filt_cycle, filt_allele = filt_allele,
                filt_mapq = filt_mapq, filt_baseq = filt_baseq,
                Mcy = Mcy, Ucy = Ucy, csums = csums))
}


read.umtab2.chr <- function(files, sampleNames = NULL,
                            keepM, keepU, readCycle = FALSE, keepFilt = FALSE,
                            verbose = TRUE) {
    columnHeaders <- strsplit(readLines(files[1], n = 1), "\t")[[1]]
    if("strand" %in% columnHeaders)
        stranded <- TRUE
    what0 <- replicate(length(columnHeaders), integer(0))
    names(what0) <- columnHeaders
    what0[["ref"]] <- character(0)
    if(stranded)
        what0[["strand"]] <- character(0)
    if(!readCycle) {
        what0 <- what0[! what0 %in% c("Mcy", "Ucy")]
    }
    if(!keepFilt) {
        what0 <- what0[!grepl("^filt", what0)]
    }
    if(verbose) cat("[read.umtab2.chr] reading", files[1], "\n")
    scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,
                     what = what0, na.strings = c("NA", "?"))
    intab <- do.call(scan, c(scanPars, file = files[1]))
    if(length(intab[[1]]) == 0)
        return(list())
    pos <- intab[["off"]]
    chr <- intab[["ref"]]
    nPos <- length(pos)
    nSamples <- length(files)
    M <- matrix(0L, nrow = nPos, ncol = nSamples)
    colnames(M)  <- sampleNames
    U <- M
    if(readCycle)
        Mcy <- Ucy <- M
    else
        Mcy <- Ucy <- NULL
    if(keepFilt)
        filt_cycle <- filt_allele <- filt_mapq <- filt_baseq <- M
    else
        filt_cycle <- filt_allele <- filt_mapq <- filt_baseq <- NULL
    allMnames <- grep("^M[[:digit:]][[:digit:]]*", names(intab), value = TRUE)
    allUnames <- grep("^U[[:digit:]][[:digit:]]*", names(intab), value = TRUE)
    csums <- matrix(0, nrow = length(c(allMnames, allUnames)), ncol = length(files))
    colnames(csums) <- sampleNames
    rownames(csums) <- c(allMnames, allUnames)
    if(missing(keepM)) keepM <- allMnames
    if(missing(keepU)) keepU <- allUnames
    stopifnot(all(c(keepM, keepU) %in% c(allMnames, allUnames)))
    M[,1] <- as.integer(Reduce("+", intab[keepM]))
    U[,1] <- as.integer(Reduce("+", intab[keepU]))
    csums[,1] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))
    if(readCycle) {
        Mcy[,1] <- intab[["Mcy"]]
        Ucy[,1] <- intab[["Ucy"]]
    }
    if(keepFilt) {
        filt_cycle[,1] <- intab[["filt_cycle"]]
        filt_allele[,1] <- intab[["filt_allele"]]
        filt_mapq[,1] <- intab[["filt_mapq"]]
        filt_baseq[,1] <- intab[["filt_baseq"]]
    }
    for(ii in seq_along(files[-1]) + 1) {
        if(verbose) cat("reading", files[ii], "\n")
        intab <- do.call(scan, c(scanPars, file = files[ii]))
        if(length(intab[[1]]) == 0)
            next
        stopifnot(all(pos == intab[["Off"]]) && all(chr == intab[["Chr"]]))
        M[,ii] <- as.integer(Reduce("+", intab[keepM]))
        U[,ii] <- as.integer(Reduce("+", intab[keepU]))
        csums[,ii] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))
        if(readCycle) {
            Mcy[,ii] <- intab[["Mcy"]]
            Ucy[,ii] <- intab[["Ucy"]]
        }
        if(keepFilt) {
            filt_cycle[,ii] <- intab[["filt_cycle"]]
            filt_allele[,ii] <- intab[["filt_allele"]]
            filt_mapq[,ii] <- intab[["filt_mapq"]]
            filt_baseq[,ii] <- intab[["filt_baseq"]]
        }
    }
    return(list(chr = chr, pos = pos, M = M, U = U,
                Mcy = Mcy, Ucy = Ucy,
                filt_cycle = filt_cycle, filt_allele = filt_allele,
                filt_mapq = filt_mapq, filt_baseq = filt_baseq,
                csums = csums))
}

read.umtab2.chr2 <- function(files, sampleNames = NULL,
                             keepM, keepU, verbose = TRUE) {
    columnHeaders <- strsplit(readLines(files[1], n = 1), "\t")[[1]]
    if("strand" %in% columnHeaders)
        stranded <- TRUE
    what0 <- replicate(length(columnHeaders), integer(0))
    names(what0) <- columnHeaders
    what0[["ref"]] <- character(0)
    if(stranded)
        what0[["strand"]] <- character(0)
    if(verbose) cat("[read.umtab2.chr] parsing all samples first\n")
    scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,
                     what = what0, na.strings = c("NA", "?"))
    ## First we get the locations of everything in all files
    getLocs <- function(file) {
        if(stranded) {
            intab <- scan(file, what = list("ref" = character(0),
                                "off" = integer(0),
                                "strand" = character(0)),
                          flush = TRUE, skip = 1, quote = "", sep = "")
            gr <- GRanges(seqnames = intab[["ref"]],
                          strand = ifelse(intab[["strand"]]== "W", 1L, -1L),
                          ranges = IRanges(start = intab[["off"]], width = 1))
        } else {
            intab <- scan(file, what = list("ref" = character(0),
                                "off" = integer(0)),
                          flush = TRUE, quote = "", sep = "")
            gr <- GRanges(seqnames = intab[["ref"]],
                          ranges = IRanges(start = intab[["pos"]], width = 1))
        }
        gr
    }
    grLoc <- as.list(files)
    grLoc[[1]] <- getLocs(files[1])
    grLoc <- Reduce(function(gr, file) {
        union(gr, getLocs(file))
    }, grLoc)
    nPos <- length(grLoc)
    nSamples <- length(files)

    M <- matrix(0L, nrow = nPos, ncol = nSamples)
    colnames(M)  <- sampleNames
    U <- Mcy <- Ucy <- M
    filt_cycle <- filt_allele <- filt_mapq <- filt_baseq <- M
    allMnames <- grep("^M[[:digit:]][[:digit:]]*", columnHeaders, value = TRUE)
    allUnames <- grep("^U[[:digit:]][[:digit:]]*", columnHeaders, value = TRUE)
    csums <- matrix(0, nrow = length(c(allMnames, allUnames)), ncol = length(files))
    colnames(csums) <- sampleNames
    rownames(csums) <- c(allMnames, allUnames)
    if(missing(keepM)) keepM <- allMnames
    if(missing(keepU)) keepU <- allUnames
    stopifnot(all(c(keepM, keepU) %in% c(allMnames, allUnames)))

    for(ii in seq_along(files)) {
        if(verbose) cat("[read.umtab2.chr] reading", files[ii], "\n")
        intab <- do.call(scan, c(scanPars, file = files[ii]))
        if(length(intab[[1]]) == 0)
            next
        if(stranded) {
            gr <- GRanges(seqnames = intab[["ref"]],
                          strand = ifelse(intab[["strand"]] == "W", 1L, -1L),
                          ranges = IRanges(start = intab[["off"]], width = 1))
        } else {
            gr <- GRanges(seqnames = intab[["ref"]],
                          ranges = IRanges(start = intab[["off"]], width = 1))
        }
        sh <- subjectHits(findOverlaps(gr, grLoc))
        stopifnot(length(sh) == length(intab$ref) & !anyDuplicated(sh))
        M[sh,ii] <- as.integer(Reduce("+", intab[keepM]))
        U[sh,ii] <- as.integer(Reduce("+", intab[keepU]))
        Mcy[sh,ii] <- intab[["Mcy"]]
        Ucy[sh,ii] <- intab[["Ucy"]]
        filt_cycle[sh,ii] <- intab[["filt_cycle"]]
        filt_allele[sh,ii] <- intab[["filt_allele"]]
        filt_mapq[sh,ii] <- intab[["filt_mapq"]]
        filt_baseq[sh,ii] <- intab[["filt_baseq"]]
        csums[,ii] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))
    }
    return(list(chr = as.character(seqnames(grLoc)), pos = start(grLoc), M = M, U = U,
                Mcy = Mcy, Ucy = Ucy,
                filt_cycle = filt_cycle, filt_allele = filt_allele,
                filt_mapq = filt_mapq, filt_baseq = filt_baseq,
                csums = csums))
}



read.umtab <- function(dirs, sampleNames = NULL, rmZeroCov = FALSE,
                       pattern = NULL,
                       keepU = c("U10", "U20", "U30", "U40"),
                       keepM = c("M10", "M20", "M30", "M40"), verbose = TRUE) {
    filesPerDir <- lapply(dirs, function(xx) sort(list.files(xx, pattern = pattern, full.names = TRUE)))
    if(!all(sapply(filesPerDir, function(xx) all(basename(xx) == basename(filesPerDir[[1]])))))
        warning("'dirs' does not contain the exact same file names")
    allChrs <- as.list(seq_along(filesPerDir[[1]]))
    for(ii in seq_along(filesPerDir[[1]])) {
        chrRead <- read.umtab.chr(sapply(filesPerDir, function(xx) xx[ii]),
                                  keepM = keepM, keepU = keepU, sampleNames = sampleNames,
                                  verbose = verbose)
        if(rmZeroCov && length(chrRead) != 0){
            wh <- which(rowSums2(chrRead$U + chrRead$M) > 0)
            chrRead$M <- chrRead$M[wh,]
            chrRead$U <- chrRead$U[wh,]
            chrRead$Mcy <- chrRead$Mcy[wh,]
            chrRead$Ucy <- chrRead$Ucy[wh,]
            chrRead$chr <- chrRead$chr[wh]
            chrRead$pos <- chrRead$pos[wh]
            chrRead$Map <- chrRead$Map[wh]
            chrRead$GC <- chrRead$GC[wh]
        }
        allChrs[[ii]] <- chrRead
    }
    BSdata <- BSseq(chr = do.call(c, lapply(allChrs, function(xx) xx$chr)),
                    pos = do.call(c, lapply(allChrs, function(xx) xx$pos)),
                    M = do.call(rbind, lapply(allChrs, function(xx) xx$M)),
                    Cov = do.call(rbind, lapply(allChrs, function(xx) xx$M)) +
                    do.call(rbind, lapply(allChrs, function(xx) xx$U)),
                    rmZeroCov = rmZeroCov)
    GC <- do.call(c, lapply(allChrs, function(xx) xx$GC))
    Map <- do.call(c, lapply(allChrs, function(xx) xx$Map))
    Mcy <- do.call(rbind, lapply(allChrs, function(xx) xx$Mcy))
    Ucy <- do.call(rbind, lapply(allChrs, function(xx) xx$Ucy))
    csums <- lapply(allChrs, function(xx) xx$csums)
    names(csums) <- basename(filesPerDir[[1]])
    return(list(BSdata = BSdata, GC = GC, Map = Map,
                Mcy = Mcy, Ucy = Ucy, csums = csums))
}


read.umtab.chr <- function(files, sampleNames = NULL,
                           keepM = c("M10", "M20", "M30", "M40"),
                           keepU = c("U10", "U20", "U30", "U40"), verbose = TRUE) {
    columnHeaders <- c("Chr", "Off", "M0", "M10", "M20", "M30", "M40", "Mqual", "Mcy",
                       "U0", "U10", "U20", "U30", "U40", "Uqual", "Ucy",
                       "MapaFW", "MapaRC", "CGContFW", "CGContRC")
    what0 <- c(list(character(0)), replicate(19, integer(0)))
    names(what0) <- columnHeaders
    line1 <- readLines(files[1], n = 1)
    if(length(line1) == 0)
        return(list())
    line1 <- strsplit(line1, "\t")[[1]]
    stopifnot(all(line1 == columnHeaders))
    if(verbose) cat("[read.umtab.chr] reading", files[1], "\n")
    scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,
                  what = what0, na.strings = c("NA", "?"))
    intab <- do.call(scan, c(scanPars, file = files[1]))
    if(length(intab[[1]]) == 0)
        return(list())
    GC <- as.integer(intab[["CGContFW"]] + intab[["CGContRC"]])
    Map <- as.integer(intab[["MapaFW"]] + intab[["MapaRC"]])
    pos <- intab[["Off"]]
    chr <- intab[["Chr"]]
    nPos <- length(pos)
    nSamples <- length(files)
    M <- matrix(0L, nrow = nPos, ncol = nSamples)
    U <- matrix(0L, nrow = nPos, ncol = nSamples)
    Mcy <- matrix(0L, nrow = nPos, ncol = nSamples)
    Ucy <- matrix(0L, nrow = nPos, ncol = nSamples)
    csums <- matrix(0, nrow = 10, ncol = length(files))
    colnames(M) <- colnames(U) <- colnames(csums) <- sampleNames
    allMnames <- c("M0", "M10", "M20", "M30", "M40")
    allUnames <- c("U0", "U10", "U20", "U30", "U40")
    stopifnot(all(c(keepM, keepU) %in% c(allMnames, allUnames)))
    rownames(csums) <- c(allMnames, allUnames)
    M[,1] <- as.integer(Reduce("+", intab[keepM]))
    U[,1] <- as.integer(Reduce("+", intab[keepU]))
    Mcy[,1] <- intab[["Mcy"]]
    Ucy[,1] <- intab[["Ucy"]]
    csums[,1] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))
    for(ii in seq_along(files[-1]) + 1) {
        if(verbose) cat("[read.umtab2] reading", files[ii], "\n")
        intab <- do.call(scan, c(scanPars, file = files[ii]))
        if(length(intab[[1]]) == 0)
            next
        stopifnot(all(pos == intab[["Off"]]) && all(chr == intab[["Chr"]]))
        M[,ii] <- as.integer(Reduce("+", intab[keepM]))
        U[,ii] <- as.integer(Reduce("+", intab[keepU]))
        Mcy[,ii] <- intab[["Mcy"]]
        Ucy[,ii] <- intab[["Ucy"]]
        csums[,ii] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))
    }
    return(list(chr = chr, pos = pos, Map = Map, GC = GC, M = M, U = U,
                Mcy = Mcy, Ucy = Ucy, csums = csums))
}

read.bsmoothDirRaw <- function(dir, seqnames = NULL, keepCycle = FALSE, keepFilt = FALSE,
                               header = TRUE, verbose = TRUE) {
    dir <- normalizePath(dir)
    inpattern <- "\\.cpg\\.tsv(|\\.gz)$"
    if(length(dir) != 1 || !file.info(dir)$isdir)
        stop("argument 'dir' needs to be a single directory")
    allChrFiles <- list.files(dir, pattern = inpattern, full.names = TRUE)
    if(!is.null(seqnames))
        allChrFiles <- allChrFiles[sub(inpattern, "", basename(allChrFiles)) %in% seqnames]
    if(length(allChrFiles) == 0) {
        warning(sprintf("dir '%s' is empty or has no output from 'seqnames'", dir))
        return(NULL)
    }
    if(header) {
        columnHeaders <- sapply(allChrFiles, function(thisfile) {
            if(grepl("\\.gz$", thisfile))
                con <- gzfile(thisfile)
            else
                con <- file(thisfile, open = "r")
            out <- readLines(con, n = 1)
            close(con)
            out
        })
        columnHeaders <- strsplit(columnHeaders, "\t")
        if(!all(sapply(columnHeaders, function(xx) all.equal(columnHeaders[[1]], xx))))
            stop(sprintf("input files in dir '%s' does not have the same headers", dir))
        columnHeaders <- columnHeaders[[1]]
    } else
        columnHeaders <- c("ref", "off", "strand", "Mstr", "Mcy", "Ustr", "Ucy",
                           "filt_cycle", "filt_readlen", "filt_allele", "filt_mapq", "filt_baseq")
    what0 <- replicate(length(columnHeaders), character(0))
    names(what0) <- columnHeaders
    int <- c(which(columnHeaders %in% c("off", "Mcy", "Ucy")), grep("^filt", columnHeaders))
    what0[int] <- replicate(length(int), integer(0))
    if(!keepCycle)
        what0[c("Mcy", "Ucy")] <- replicate(2, NULL)
    if(!keepFilt)
        what0[grep("^filt", names(what0))] <- replicate(length(grep("^filt", names(what0))), NULL)
    outList <- lapply(allChrFiles, function(thisfile) {
        if(verbose)
            cat(sprintf("[read.bsmoothDirRaw] Reading '%s'\n", thisfile))
        if(grepl("\\.gz$", thisfile))
            con <- gzfile(thisfile)
        else
            con <- file(thisfile)
        out <- scan(con, skip = header, what = what0, sep = "\t",
                    quote = "", na.strings = "NA", quiet = TRUE)
        close(con)
        out
    })
    listNames <- names(outList[[1]])
    names(listNames) <- listNames
    out <- lapply(listNames, function(name) {
        do.call(c, lapply(outList, function(xx) xx[[name]]))
    })
    rm(outList)
    gr <- GRanges(seqnames = paste0("chr", out[["ref"]]),
                  ranges = IRanges(start = out[["off"]], width = 1))
    out[["ref"]] <- out[["off"]] <- NULL
    names(out)[names(out) == "strand"] <- "bstrand"
    out <- out[!sapply(out, is.null)]
    df <- DataFrame(out)
    mcols(gr) <- df
    gr
}


sampleRawToBSseq <- function(gr, qualityCutoff = 20, sampleName = NULL, rmZeroCov = FALSE) {
    numberQualsGreaterThan <- function(cvec) {
        onestring <- paste(cvec, collapse = "")
        greater <- (as.integer(charToRaw(onestring)) - 33L >= qualityCutoff)
        out <- tapply(greater, rep(1:length(cvec), times = nchar(cvec)), sum)
        out
    }
    strToCov <- function(vec) {
        Cov <- rep(0, length(vec))
        wh <- which(! vec %in% c("", "0"))
        if(length(wh) > 0)
            Cov[wh] <- numberQualsGreaterThan(vec[wh])
        Cov
    }
    M <- matrix(strToCov(mcols(gr)[, "Mstr"]), ncol = 1)
    Cov <- M + strToCov(mcols(gr)[, "Ustr"])
    mcols(gr) <- NULL
    BSseq(gr = gr, M = M, Cov = Cov, sampleNames = sampleName, rmZeroCov = rmZeroCov)
}

read.bsmooth <- function(dirs, sampleNames = NULL, seqnames = NULL, returnRaw = FALSE,
                         qualityCutoff = 20, rmZeroCov = FALSE, verbose = TRUE) {
    dirs <- normalizePath(dirs, mustWork = TRUE)
    if(!(all(file.info(dirs)$isdir)))
        stop("argument 'dirs' has to be directories")
    if(anyDuplicated(dirs))
        stop("duplicate entries in 'dirs'")
    if(is.null(sampleNames) && !anyDuplicated(basename(dirs)))
        sampleNames <- basename(dirs)
    if(is.null(sampleNames))
        sampleNames <- dirs
    if(!is.null(sampleNames) && (length(sampleNames) != length(dirs) || anyDuplicated(sampleNames)))
        stop("argument 'sampleNames' (if not NULL) has to have the same length as argument 'dirs', without duplicate entries")
    idxes <- seq_along(dirs)
    names(idxes) <- sampleNames
    allOut <- lapply(idxes, function(ii) {
        if(verbose) cat(sprintf("[read.bsmooth] Reading dir '%s' ... ", dirs[ii]))
        ptime1 <- proc.time()
        if(returnRaw) {
            out <- read.bsmoothDirRaw(dir = dirs[ii], seqnames = seqnames, keepCycle = TRUE,
                                      keepFilt = TRUE, header = TRUE, verbose = FALSE)
        } else {
            raw <- read.bsmoothDirRaw(dir = dirs[ii], seqnames = seqnames, keepCycle = FALSE,
                                      keepFilt = FALSE, header = TRUE, verbose = FALSE)
            out <- sampleRawToBSseq(raw, qualityCutoff = qualityCutoff, rmZeroCov, sampleName = sampleNames[ii])
        }
        ptime2 <- proc.time()
        stime <- (ptime2 - ptime1)[3]
        if(verbose) cat(sprintf("done in %.1f secs\n", stime))
        out
    })
    if(!returnRaw) {
        if(verbose) cat(sprintf("[read.bsmooth] Joining samples ... "))
        ptime1 <- proc.time()
        allOut <- combineList(allOut)
        ptime2 <- proc.time()
        stime <- (ptime2 - ptime1)[3]
        if(verbose) cat(sprintf("done in %.1f secs\n", stime))
    }
    allOut
}

parsingPipeline <- function(dirs, qualityCutoff = 20, outDir, seqnames = NULL,
                            subdir = "ev_bt2_cpg_tab", timing = FALSE) {
    if(!all(file.exists(dirs)))
        stop("not all directories in 'dirs' exists.")
    cat("[parsingPipeline] Parsing all files.\n")
    oneDir <- function(dir) {
        cat("[parsingPipeline]  dir", basename(dir), ": ")
        base <- basename(dir)
        cat("parsing, ")
        ptime1 <- proc.time()
        raw <- read.bsmoothDirRaw(file.path(dir, subdir),
                                  keepCycle = TRUE, keepFilt = TRUE,
                                  verbose = FALSE)
        assign(paste0(base, ".raw"), raw)
        ptime2 <- proc.time()
        stime <- (ptime2 - ptime1)[3]
        if(timing) {
            cat(sprintf("\ndone in %.1f secs\n", stime))
            print(gc())
        }
        cat("saving, ")
        save(list = paste0(base, ".raw"),
             file = file.path(outDir, paste0(base, ".raw.rda")))
        cat("converting, ")
        ptime1 <- proc.time()
        bsseq <- sampleRawToBSseq(raw, qualityCutoff = qualityCutoff,
                                  sampleName = base)
        seqlevels(bsseq)[seqlevels(bsseq) == "chrgi|9626243|ref|NC_001416.1|"] <- "chrLambda"
        ptime2 <- proc.time()
        stime <- (ptime2 - ptime1)[3]
        if(timing) {
            cat(sprintf("\ndone in %.1f secs\n", stime))
            print(gc())
        }
        cat("ordering, ")
        ptime1 <- proc.time()
        bsseq <- chrSelectBSseq(bsseq, order = TRUE,
                                seqnames = seqnames)
        assign(paste0(base, ".bsseq"), bsseq)
        ptime2 <- proc.time()
        stime <- (ptime2 - ptime1)[3]
        if(timing) {
            cat(sprintf("\nIn %.1f secs\n", stime))
            print(gc())
        }
        cat("saving, ")
        save(list = paste0(base, ".bsseq"),
             file = file.path(outDir, paste0(base, ".bsseq.rda")))
        cat("done\n")
        NULL
    }
    for(dir in dirs)
        oneDir(dir)
    invisible(dirs)
}