R/utils.read.ped.r

In dartR.base: Analysing 'SNP' and 'Silicodart' Data - Basic Functions

#' @name utils.read.ped
#' @title An internal script [Custodian to provide a title]
#' @family utilities
#' 
#' @description 
#' WARNING: UTILITY SCRIPTS ARE FOR INTERNAL USE ONLY AND SHOULD NOT BE USED BY END USERS AS THEIR USE OUT OF CONTEXT COULD LEAD TO UNPREDICTABLE OUTCOMES.
#' @param file Custodian to provide
#' @param snps Custodian to provide
#' @param which Custodian to provide
#' @param split Custodian to provide
#' @param sep Custodian to provide
#' @param na.strings Custodian to provide
#' @param lex.order Custodian to provide
#' @param show_warnings Custodian to provide
#'  
#' #[Custodian to provide other parameter descriptions]

#' @details
#' #[Custodian to provide details]
#' 
#' @author Custodian: Luis Mijangos (Post to
#' \url{https://groups.google.com/d/forum/dartr})

# @export
#' @return The resultant genlight object
#' @importFrom snpStats switch.alleles

utils.read.ped <- function (file, 
                            # n, 
                            snps, 
                            which, 
                            split = "\t| +",
                            sep = ".", 
                            na.strings = "0", 
                            lex.order = FALSE,
                            show_warnings = TRUE){
  
# check if packages are installed
  pkg <- "snpStats"
  if (!(requireNamespace(pkg, quietly = TRUE))) {
    cat(error(
      "Package",
      pkg,
      " needed for this function to work. Please install it.\n"
    ))
    return(-1)
  }



  r0 <- as.raw(0)
  r1 <- as.raw(1)
  r2 <- as.raw(2)
  r3 <- as.raw(3)
  con <- gzfile(file)
  open(con)
  # if (missing(n)) {
  n <- 0
  repeat {
    line <- readLines(con, n = 1)
    if (length(line) == 0) 
      break
    n <- n + 1
  }
  if (n == 0){ 
    stop("Nothing read")
  }
  seek(con, 0)
  # }
  gen <- missing(snps)
  map <- NULL
  if (!gen) {
    m <- length(snps)
    if (m == 1) {
      map <- read.table(snps, comment.char = "")
      m <- nrow(map)
      if (missing(which)) {
        which <- 1
        repeat {
          snps <- map[, which]
          if (!any(duplicated(snps))) 
            break
          if (which == ncol(map)) 
            stop("No unambiguous snp names found on file")
          which <- which + 1
        }
      } else {
        snps <- map[, which]
      }
    }
  } else {
    line <- readLines(con, n = 1)
    fields <- strsplit(line, split)[[1]]
    nf <- length(fields)
    if (nf%%2 != 0) 
      stop("Odd number of fields")
    m <- (nf - 6)/2
    seek(con, 0)
  }
  nf <- 6 + 2 * m
  result <- matrix(raw(n * m), nrow = n)
  ped <- character(n)
  mem <- character(n)
  pa <- character(n)
  ma <- character(n)
  sex <- numeric(n)
  aff <- numeric(n)
  rownms <- character(n)
  a1 <- a2 <- rep(NA, m)
  a1m <- a2m <- rep(TRUE, m)
  mallelic <- rep(FALSE, m)
  for (i in 1:n) {
    line <- readLines(con, n = 1)
    fields <- strsplit(line, "\t| +")[[1]]
    to.na <- fields %in% na.strings
    fields[to.na] <- NA
    ped[i] <- fields[1]
    mem[i] <- fields[2]
    pa[i] <- fields[3]
    ma[i] <- fields[4]
    sex[i] <- suppressWarnings(as.numeric(fields[5]))
    aff[i] <- as.numeric(fields[6])
    alleles <- matrix(fields[7:nf], byrow = TRUE, ncol = 2)
    one <- two <- rep(FALSE, m)
    for (k in 1:2) {
      ak <- alleles[, k]
      akm <- is.na(ak)
      br1 <- !akm & a1m
      a1[br1] <- ak[br1]
      a1m[br1] <- FALSE
      br2 <- !akm & (a1 == ak)
      one[br2] <- TRUE
      br3 <- !akm & !a1m & (a1 != ak)
      br4 <- br3 & a2m
      a2[br4] <- ak[br4]
      a2m[br4] <- FALSE
      br5 <- br3 & (a2 == ak)
      two[br5] <- TRUE
      mallelic <- mallelic | !(akm | one | two)
    }
    gt <- rep(r0, m)
    gt[one & !two] <- r1
    gt[one & two] <- r2
    gt[two & !one] <- r3
    result[i, ] <- gt
  }
  close(con)
  if (any(a1m & show_warnings==T)){ 
    warning("no data for ", sum(a1m), " loci")
  }
  mono <- (a2m & !a1m)
  if (any(mono & show_warnings==T)){ 
    warning(sum(mono), " loci were monomorphic")
  }
  if (any(mallelic & show_warnings==T)) {
    result[, mallelic] <- r0
    warning(sum(mallelic), " loci were multi-allelic --- set to NA")
  }
  if (gen){ 
    snps <- paste("locus", 1:m, sep = sep)
  }
  if (any(duplicated(ped))) {
    if (any(duplicated(mem))) {
      rnames <- paste(ped, mem, sep = sep)
      if (any(duplicated(rnames))){ 
        stop("could not create unique subject identifiers")
      }
    } else{
      rnames <- mem
    }
  }  else {
    rnames <- ped
  }
  dimnames(result) <- list(rnames, snps)
  result <- new("SnpMatrix", result)
  if (lex.order) {
    swa <- (!(is.na(a1) | is.na(a2)) & (a1 > a2))
    snpStats::switch.alleles(result, swa)
    a1n <- a1
    a1n[swa] <- a2[swa]
    a2[swa] <- a1[swa]
    a1 <- a1n
  }
  fam <- data.frame(row.names = rnames, pedigree = ped, member = mem, 
                    father = pa, mother = ma, sex = sex, affected = aff)
  if (is.null(map)){ 
    map <- data.frame(row.names = snps, snp.name = snps, 
                      allele.1 = a1, allele.2 = a2)
  } else {
    map$allele.1 <- a1
    map$allele.2 <- a2
    names(map)[which] <- "snp.names"
  }
  list(genotypes = result, fam = fam, map = map)
}

rotate.matrix <- function (x, 
                           angle = 10, 
                           method = "bilinear"){
  
  img <- x
  angle.rad <- angle * pi/180
  co.x <- matrix(rep(-(ncol(img)/2 - 0.5):(ncol(img)/2 - 
                                             0.5), nrow(img)), nrow = nrow(img), byrow = T)
  co.y <- matrix(rep(-(nrow(img)/2 - 0.5):(nrow(img)/2 - 
                                             0.5), ncol(img)), ncol = ncol(img))
  co.xn <- round(co.x * cos(angle.rad) - co.y * sin(angle.rad))
  co.yn <- round(co.x * sin(angle.rad) + co.y * cos(angle.rad))
  co.xn2 <- co.xn + max(co.xn) + 1
  co.yn2 <- co.yn + max(co.yn) + 1
  img.rot <- numeric(max(co.yn2) * max(co.xn2))
  img.rot[(co.xn2 - 1) * max(co.yn2) + co.yn2] <- img
  dim(img.rot) <- c(max(co.yn2), max(co.xn2))
  attr(img.rot, "bits.per.sample") <- attr(img, "bits.per.sample")
  attr(img.rot, "samples.per.pixel") <- attr(img, "samples.per.pixel")
  return(img.rot)
  
}