R/preprocessSingleCell.R

In rhondabacher/methylscaper: Visualization of Methylation Data

#' Process single-cell data
#'
#' This function subsets the data and prepares it for visualizing by
#' generating representation methylation-state matrices from
#' single-cell methylation data (for example, sc-MNT data).
#' We assume the data has already been preprocess using the subsetSC
#' function in methylscaper. See the vignette for a more thorough
#' explanation of each parameter. The output should be passed directly to the
#' plotting function.
#'
#' @param dataIn A list object containing two elements labelled gch and hcg (already pre-processed.)
#' @param startPos The index of the first position to include
#'  in the visualization. If using this within the R console it is
#'  recomended to specify the start and end directly.
#'  In the Shiny app, a slider will let the user refine these positions.
#' @param endPos The index of the final position to include in the visualization.
#' @param updateProgress A function for generating progress bars in the Shiny app.
#'   Should be left NULL otherwise.
#' @return The output is a list containing the elements 'gch' and 'hcg.
#'      Each is a dataframe with reads/cells on the rows and each column
#'      is a base-pair. The matrix is coded as follows:
#'          -2: unmethylated GCH or HCG site
#'          -1: base pairs between two unmethylated GCH or HCG sites
#'          0: base pairs between mismatching methylation states of
#'              two GCH or HCG sites
#'          1: base pairs between two methylated GCH or HCG sites
#'          2: methylated GCH or HCG site
#' @importFrom utils tail
#' @import data.table
#' @importFrom methods is
#' @export
#' @examples
#'
#' data(singlecell_subset)
#' prepsc.out <- prepSC(singlecell_subset,
#'     startPos = 105636488, endPos = 105636993
#' )
prepSC <- function(
        dataIn, startPos = NULL, endPos = NULL,
        updateProgress = NULL) {
    if (is(dataIn, "SummarizedExperiment") | is(dataIn, "SingleCellExperiment")) {
        dataIn <- reformatSCE(dataIn)
    }

    if (is.null(updateProgress) & (endPos - startPos) >= 50000) {
        message("Selected range is longer than 50k bp,
                    plot may take a few seconds to render.")
    }

    if (is.null(updateProgress) & (endPos - startPos) >= 100000) {
        message("Selected range is longer than 100k bp,
                    this may not be optimal.
                    Plot will render but may take a few minutes.")
    }

    gc_seq_data <- dataIn$gch
    cg_seq_data <- dataIn$hcg
    if (is.function(updateProgress)) {
        updateProgress(message = "Filtering HCG (endogenous methylation)
                                    site data", value = 0.1)
    }

    cg_seq_sub <- lapply(cg_seq_data, function(x) {
        tempSub <- x[order(x$pos), ]
        tempSub <- subset(tempSub, tempSub$pos >= startPos & tempSub$pos <= endPos)
        return(tempSub)
    })

    if (is.function(updateProgress)) {
        updateProgress(message = "Filtering GCH (accessibility/occupancy) site data", value = 0.5)
    }


    gc_seq_sub <- lapply(gc_seq_data, function(x) {
        tempSub <- x[order(x$pos), ]
        tempSub <- subset(tempSub, tempSub$pos >= startPos & tempSub$pos <= endPos)
        return(tempSub)
    })

    all_cg_sites <- unique(do.call(c, cg_seq_sub))
    all_gc_sites <- unique(do.call(c, gc_seq_sub))
    useseq <- intersect(
        which(vapply(cg_seq_sub, function(x) nrow(x), numeric(1)) > 1),
        which(vapply(gc_seq_sub, function(x) nrow(x), numeric(1)) > 1)
    )

    if (length(useseq) == 0) {
        return(NULL)
    }

    cg_seq_sub <- cg_seq_sub[useseq]
    gc_seq_sub <- gc_seq_sub[useseq]

    if (is.function(updateProgress)) {
        updateProgress(message = "Mapping HCG (endogenous methylation) site data", value = 0.75)
    }

    cg_outseq <- lapply(cg_seq_sub, function(x) {
        if (nrow(x) == 0) {
            NULL
        } else {
            mapSC(x, startPos, endPos)
        }
    })

    cg_outseq <- cg_outseq[!vapply(cg_outseq, is.null, logical(1))]
    hcg <- data.matrix(do.call(rbind, cg_outseq))
    rownames(hcg) <- as.character(seq(1, nrow(hcg)))

    if (is.function(updateProgress)) {
        updateProgress(message = "Mapping GCH (accessibility/occupancy) site data", value = 0.9)
    }

    gc_outseq <- lapply(gc_seq_sub, function(x) {
        if (nrow(x) == 0) {
            NULL
        } else {
            mapSC(x, startPos, endPos)
        }
    })

    gc_outseq <- gc_outseq[!vapply(gc_outseq, is.null, logical(1))]
    gch <- data.matrix(do.call(rbind, gc_outseq))
    rownames(gch) <- as.character(seq(1, nrow(gch)))

    return(list(gch = gch, hcg = hcg))
}

mapSC <- function(inSeq, startPos, endPos) {
    inSeq$basePos <- as.numeric(inSeq$pos) - startPos + 1
    fill_1 <- seq(startPos, endPos) - startPos + 1
    someMethyl <- which(inSeq$rate > 0)
    noMethyl <- which(inSeq$rate <= 0)
    fill_1[fill_1 %in% inSeq[someMethyl, ]$basePos] <- 2
    fill_1[fill_1 %in% inSeq[noMethyl, ]$basePos] <- -2
    fill_1[abs(fill_1) != 2] <- "."

    sites <- inSeq$basePos
    sites <- sites[sites > 0]
    editseq <- fill_1
    sites_temp <- c(0, sites, max(sites) + 1)

    for (j in seq(1, (length(sites_temp) - 1))) {
        if (sites_temp[j + 1] == 1) { # skip
        } else {
            tofill <- seq(sites_temp[j] + 1, (sites_temp[j + 1] - 1))
            s1 <- editseq[pmax(1, sites_temp[j])]
            s2 <- editseq[pmin(length(editseq), sites_temp[j + 1])]

            if (s1 == "2" & s2 == "2") {
                fillvec <- 1
            } else if (s1 == "2" & s2 == "-2") {
                fillvec <- 0
            } else if (s1 == "-2" & s2 == "2") {
                fillvec <- 0
            } else if (s1 == "-2" & s2 == "-2") {
                fillvec <- -1
            } else {
                fillvec <- 0
            }
            fillvec <- rep(fillvec, length(tofill))
            editseq[tofill] <- fillvec
        }
    }
    return(editseq)
}



#' Load in methylation data
#'
#' This function loads the single-cell files. It takes a path to the data files
#' and a chromosome number as arguments and returns the desired subset of the
#' data. Processing by chromosome is important for speed and memory efficiency.
#' The input files should be tab separated with three columns.
#' The first column is the chromosome, the second is the position (basepair), and the third
#' is the methylation indicator/rate. The folder should contain two subfolders titled
#' met and acc, with the endogenous methylation and accessibility methylation files,
#' respectively.
#'
#' @param path Path to the folder containing the single-cell files.
#' @param chromosome The chromosome to subset the files to.
#' @param startPos The index of the first position to include
#'  in the subsetting. This is optional as further narrowing of the
#'  position can be done in the visualization step/tab.
#'  In the Shiny app, a slider will let the user refine the positions.
#' @param endPos The index of the final position to include in subset.
#' @param updateProgress A function for generating progress bars in the Shiny app.
#'   Should be left NULL otherwise.
#' @return The output is RDS files that can be loaded into the visualization
#'  tab on the Shiny app
#' @import data.table
#' @export
#'
#' @examples
#' # example not run since needs directory input from user
#' # subsc.out <- subsetSC("filepath", chromosome=19)
subsetSC <- function(
        path, chromosome, startPos = NULL,
        endPos = NULL, updateProgress = NULL) {
    if (!is.list(path)) {
        cgfiles <- paste0(path, "/", "met/", sort(grep("met", list.files(paste0(path, "/met")), value = TRUE)))
        gcfiles <- paste0(path, "/", "acc/", sort(grep("acc", list.files(paste0(path, "/acc")), value = TRUE)))
    } else {
        cgfiles <- path[[1]]
        gcfiles <- path[[2]]

        getord <- order(path[[3]])
        cgfiles <- cgfiles[getord]

        getord <- order(path[[4]])
        gcfiles <- gcfiles[getord]
    }


    if (length(cgfiles) != length(gcfiles)) {
        stop("Must have the same number of methylation and
                accessibility files!")
    }

    useChr <- chromosome

    cg_seq <- list()
    for (i in seq(1, length(cgfiles))) {
        in_cg_seq <- fread(cgfiles[i],
            header = FALSE, stringsAsFactors = FALSE
        )
        # A couple possible formats from bismark:
        if (ncol(in_cg_seq) %in% c(4, 6)) {
            in_cg_seq <- in_cg_seq[, c(1, 2, 4)]
        }

        if (ncol(in_cg_seq) != 3) {
            stop("Data is not formatted correctly. See vignette for details.")
        }

        if (in_cg_seq[1, 1] == "chr") {
            in_cg_seq <- in_cg_seq[-1, ]
        }
        colnames(in_cg_seq) <- c("chr", "pos", "rate")
        in_cg_seq$pos <- as.numeric(in_cg_seq$pos)
        in_cg_seq$rate <- as.numeric(in_cg_seq$rate)

        in_cg_seq <- subset(in_cg_seq, in_cg_seq$chr == useChr)
        if (!is.null(startPos) & !is.null(endPos)) {
            in_cg_seq <- subset(in_cg_seq, in_cg_seq$pos >= startPos & in_cg_seq$pos <= endPos)
        }

        cg_seq[[i]] <- in_cg_seq
        if (is.function(updateProgress)) {
            updateProgress(message = "Reading HCG (endogenous methylation) site files", value = i / length(cgfiles))
        }
    }

    gc_seq <- list()
    for (i in seq(1, length(gcfiles))) {
        in_gc_seq <- fread(gcfiles[i], header = FALSE, stringsAsFactors = FALSE)
        if (ncol(in_gc_seq) %in% c(4, 6)) {
            in_gc_seq <- in_gc_seq[, c(1, 2, 4)]
        }
        if (ncol(in_gc_seq) != 3) {
            stop("Data is not formatted correctly. See vignette for details.")
        }

        colnames(in_gc_seq) <- c("chr", "pos", "rate")
        if (in_gc_seq[1, 1] == "chr") {
            in_gc_seq <- in_gc_seq[-1, ]
        }
        colnames(in_gc_seq) <- c("chr", "pos", "rate")
        in_gc_seq$pos <- as.numeric(in_gc_seq$pos)
        in_gc_seq$rate <- as.numeric(in_gc_seq$rate)

        in_gc_seq <- subset(in_gc_seq, in_gc_seq$chr == useChr)
        if (!is.null(startPos) & !is.null(endPos)) {
            in_gc_seq <- subset(in_gc_seq, in_gc_seq$pos >= startPos & in_gc_seq$pos <= endPos)
        }

        gc_seq[[i]] <- in_gc_seq
        if (is.function(updateProgress)) {
            updateProgress(message = "Reading GCH (accessibility/occupancy) site files", value = i / length(gcfiles))
        }
    }

    list(hcg = cg_seq, gch = gc_seq)
}