R/scan1_summary.R

In jendelman/diaQTL: QTL Analysis in Diallel Populations

#' Summary of scan1 result
#' 
#' Summary of scan1 result
#' 
#' Plots the "score" (-deltaDIC) for each marker vs. genome position. The \code{thresh} argument should be a positive number.
#' 
#' @param scan1_data output from scan1
#' @param thresh optional -deltaDIC threshold for plotting
#' @param chrom optional, subset of chromosomes to plot 
#' @param position Either "cM" (default) or "bp"
#' 
#' @return List containing 
#' \describe{
#' \item{peaks}{data frame of markers with the highest score on each chromosome}
#' \item{plot}{ggplot object}
#' }
#' @examples
#' \dontrun{
#'   scan1_summary( scan1_example )
#'   scan1_summary( scan1_example, chrom = "10" )
#'   scan1_summary( scan1_example, chrom = c( "10", "12" ) ) 
#'   scan1_summary( scan1_example, chrom = "10", thresh = 20)
#'   }
#' @export
#' @import ggplot2
#' @importFrom rlang .data

scan1_summary <- function(scan1_data,
                          thresh=NULL,
                          chrom = NULL,
                          position = "cM") {

  stopifnot(position %in% colnames(scan1_data))
  stopifnot(thresh > 0)
  
  if(!is.null(chrom)) {
    scan1_data <- scan1_data[scan1_data$chrom %in% chrom,]
  }
  if(position=="bp") {
    x <- scan1_data[,position]/1e6
    x.label <- "Position (Mb)"
  } else {
    x <- scan1_data[,position]
    x.label <- "Position (cM)"
  }
  
  allchr <- unique(scan1_data$chrom)
  nchr <- length(allchr)
  #if (statistic=="DIC") {
    y.col <- "deltaDIC"
    y.label <- "-\U0394 DIC"
    y.sign <- -1
  #} else {
  #  y.col <- "LL"
  #  y.label <- "LL"
  #  y.sign <- 1
  #}
  
  k <- integer(nchr)
  for (i in 1:nchr) {
    y <- y.sign * scan1_data[,y.col]
    y[scan1_data$chrom!=allchr[i]] <- NA
    k[i] <- which.max(y)
  }

  p <- NULL
  if (nchr==1) {
    plotme <- data.frame(x=x,
                         y=y.sign*scan1_data[,y.col],
                         chrom=scan1_data$chrom) 
    p <- ggplot(data=plotme,aes(x=.data$x,y=.data$y)) +
        geom_line(color="#440154") +
        ylab(y.label) + 
        theme_bw() +
        theme(text = element_text(size=13),panel.grid = element_blank()) +
        xlab(x.label)
  } else {
    col <- ifelse(as.integer(factor(scan1_data$chrom))%%2==1,"1","0")
    x <- get_x(map=data.frame(chrom=scan1_data$chrom,position=x,stringsAsFactors = F))
    plotme <- data.frame(x=x,y=y.sign*scan1_data[,y.col],col=col)
    breaks <- (tapply(x,scan1_data$chrom,max) + tapply(x,scan1_data$chrom,min))/2
    p <- ggplot(data=plotme,aes(x=.data$x,y=.data$y,colour=.data$col)) +
        ylab(y.label) +
        theme_bw() +
        scale_x_continuous(name="Chromosome",breaks=breaks,labels=allchr) +
        scale_colour_manual(values=c("#21908c","#440154")) + 
        theme(text = element_text(size=13),panel.grid = element_blank(),legend.position = "none")
    
    for (i in allchr){
      ix <- which(scan1_data$chrom==i)
      p <- p + geom_line(data=plotme[ix,])
    }
  }
  
  
  if (!is.null(thresh)) {
    p <- p + geom_hline(yintercept = thresh,linetype="dashed",colour="#B89600")
  }
  
  return(list(peaks=scan1_data[k,],plot=p))
}