R/labelManhattan.R

In gap: Genetic Analysis Package

#' Annotate Manhattan or Miami Plot
#'
#' @param chr A vector of chromosomes for the markers to be labelled.
#' @param pos A vector of positions on the chromosome for the markers to be labelled. These must correspond to markers in the GWAS dataset used to make the manhattan plot.
#' @param name A vector of labels to be added next to the points specified by \env{chr} and \env{pos}.
#' @param gwas The same GWAS dataset used to plot the existing Manhattan or Miami plot to be annotated.
#' @param gwasChrLab The name of the column in \env{gwas} containing chromosome number. Defaults to \samp{"chr"}.
#' @param gwasPosLab The name of the column in \env{gwas} containing position. Defaults to \samp{"pos"}.
#' @param gwasPLab The name of the column in \env{gwas} containing p-values. Defaults to \samp{"p"}.
#' @param gwasZLab The name of the column in \env{gwas} containing z-values. Defaults to \samp{"NULL"}.
#' @param chrmaxpos Data frame containing x coordinates for chromosome start positions, generated by [`labelManhattan`].
#' @param textPos An integer or vector dictating where the label should be plotted relative to each point. Good for avoiding overlapping labels. Provide an integer to plot all points in the same relative position or use a vector to specify position for each label. Passed to the \option{pos} option of [`graphics::text`]. Defaults to \samp{4}.
#' @param angle An integer or vector dictating the plot angle of the label for each point. rovide an integer to plot all points in the same relative position or use a vector to specify position for each label.Passed to the \option{srt} option of [`graphics::text`]. Defaults to \samp{0}.
#' @param miamiBottom If \samp{TRUE}, labels will be plotted on the lower region of a Miami plot. If \samp{FALSE}, labels will be plotted on the upper region. Defaults to \samp{FALSE}.
#' @export
#'
#' @details
#' Add labels beside specified points on a Manhattan or Miami plot. Ideal for adding locus names to peaks. Currently only designed to work with [`miamiplot2`].
#'
#' @return
#' Adds annotation to existing Manhattan or Miami plot
#'
#' @author Jonathan Marten
#' @note Extended to handle extreme P values.
#' @keywords GWAS Miami Manhattan annotation
#' @examples
#' \dontrun{
#' labelManhattan(c(4,5,11,19),c(9994215,16717922,45538760,51699256),
#'                c("GENE1","GENE2","GENE3","GENE4"),
#'                gwas1,chrmaxpos=chrmaxpos)
#' labelManhattan(geneLabels$chr,geneLabel$pos,geneLabel$geneName,gwas1,chrmaxpos=chrmaxpos)
#' }

labelManhattan <- function(chr,pos,name,gwas,gwasChrLab="chr",gwasPosLab="pos",gwasPLab="p",gwasZLab="NULL",
                           chrmaxpos,textPos=4,angle=0,miamiBottom=FALSE)
{
   # Identify which columns of input files correspond to each value
   posCol 	<- which(names(gwas)%in%gwasPosLab)
   chrCol 	<- which(names(gwas)%in%gwasChrLab)
   if (is.null(gwasZLab))
      pCol 	<- which(names(gwas)%in%gwasPLab)
   else
      zCol 	<- which(names(gwas)%in%gwasZLab)

   if(length(posCol)==0){
      stop(paste0("Could not find column ",gwasPosLab," in input"))
   }
   if(length(chrCol)==0){
      stop(paste0("Could not find column ",gwasChrLab," in input"))
   }
   if(is.null(gwasZLab)){
      if(length(pCol)==0) stop(paste("Could not find column",gwasPLab,"in input"))
   } else if(length(zCol)==0) stop(paste("Could not find column",gwasZLab,"in input"))

   if(!is.data.frame(chrmaxpos)){
      stop("chrmaxpos must be an output from the miamiplot2 function")
   }
   if(!identical(names(chrmaxpos),c("chr","maxpos","genomestartpos","labpos"))){
      stop("chrmaxpos must be an output from the miamiplot2 function")
   }

   if(length(chr) != length(pos) | length(chr) != length(name)){
      stop("Chr, Pos and Name vectors must be the same length")
   }

   gwas$chrpos <- paste0(gwas[,chrCol],":",gwas[,posCol])
   chrpos <- paste0(chr,":",pos)
   row 	<- which(gwas$chrpos %in% chrpos)
   row <- match(chrpos,gwas$chrpos)
   # Make sure label order matches position order. Match would do it faster but this is how I did it the first time around and if it ain't broke...
   #v1 <- order(chrpos)
   #v2 <- order(gwas$chrpos[row])
#
   #chr <- chr[v1]
   #pos <- pos[v1]
   #name <- name[v1]
   #textPos <- textPos[v1]
   #row <- row[v2]

   # Add position-modifier values for each chromosome to position of SNP to get x-axis position
   xpos <- chrmaxpos$genomestartpos[chr] + gwas[row,posCol]
   if (is.null(gwasZLab))
      ypos <- -log10(gwas[row,pCol])
   else
      ypos <- -log10p(gwas[row,zCol])
   # Flip y axis values to plot on bottom of plot
   if(miamiBottom){
      ypos <- -ypos
   }
   # Add labels
   text(x=xpos, y=ypos, pos=textPos, labels=name, cex=1, font=3, srt=angle)
}