Skip to content
Snippets Groups Projects
Select Git revision
  • 604d22878b741001d201ec4971958e9fc22c1a32
  • master default protected
  • v0.2.0
  • v0.1.0
4 results

LR.R

Blame
  • LR.R 3.11 KiB
    
    #' Runs the LR function
    #'
    #' Returns the \code{|L||R|} value for each SNP location supplied to the function, where \code{|L|} and \code{|R|} are the number of SNPs to the left and right of the current locus within the given window \code{ws}.
    #' For more information about the \code{|L||R|} diversity statistic, please see Jacobs (2016).
    #'
    #' @param pos A numeric vector of SNP locations
    #' @param ws The window size which the \code{LR} statistic will be calculated over. This should be on the same scale as the \code{pos} vector.
    #' @param X Optional. Specify a region of the chromosome to calculate LR for in the format \code{c(startposition, endposition)}. The start position and the end position should be within the extremes of the positions given in the \code{pos} vector. If not supplied, the function will calculate LR for every SNP in the \code{pos} vector.
    #'
    #' @return A list containing the SNP positions and the \code{LR} values for those SNPs
    #' @references Jacobs, G.S., T.J. Sluckin, and T. Kivisild, \emph{Refining the Use of Linkage Disequilibrium as a Robust Signature of Selective Sweeps.} Genetics, 2016. \strong{203}(4): p. 1807
    #' @examples
    #' ## load the snps example dataset
    #' data(snps)
    #' ## run LR over all the SNPs with a window size of 3000 bp
    #' LR(snps$bp_positions,3000)
    #' ## only return results for SNPs between locations 600 and 1500 bp
    #' LR(snps$bp_positions,3000,X=c(600,1500))
    #'
    #' @export
    
    LR <- function(pos, ws, X = NULL) {
      #Check things are in the correct format
    
      #Check pos is a numeric vector
      if (is.numeric(pos) ==FALSE || is.vector(pos)==FALSE){
        stop("pos must be a numeric vector")
      }
      #Check windowsize is a number greater than 0
      if(is.numeric(ws) ==FALSE || ws <= 0){
        stop("ws must be a number greater than 0")
      }
      #If X is specified, check it is in the correct format
      if (is.null(X)==FALSE){
        if(is.numeric(X)==FALSE || is.vector(X)==FALSE){
          stop("X should be a numeric vector of length 2 e.g. c(100,200)")
        } else {
          if (length(X) != 2){
            stop("X should be a numeric vector of length 2 e.g. c(100,200)")
          } else {
            # X is in the correct format
            # Check that X will actually return a result (i.e. that the region specied by X overlaps with pos)
            if ((length(pos[pos>=X[1] & pos <= X[2]])>0) == FALSE){
              stop("The region specified by X is outside the region contained in the pos vector")
            }
          }
        }
      } else {
        # Set X equal to the extremes of pos
        X<-c(pos[1],pos[length(pos)])
      }
    
      # Set up output list
      outputLength<-length(pos[pos>=X[1] & pos <= X[2]])
      outputList<-list(position=pos[pos>=X[1] & pos <= X[2]],LR=rep(NA,outputLength))
    
      # Loop over each position in the output list and calculate LR
      for (i in 1:outputLength){
    
        # Current physical position in chromosome
        currentPos<-outputList$position[i]
    
        ## get L, R and LR
        noL <- length(pos[pos>=currentPos-ws/2 & pos < currentPos]) ## Number of SNPs to the left of the current SNP
        noR <- length(pos[pos<=currentPos+ws/2 & pos > currentPos]) ## Number of SNPs to the right of the current SNP
        outputList$LR[i]<-noL*noR
      }
      return(outputList)
    }