## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#", fig.path = "figures/", fig.height = 4.5, fig.width = 6 ) ## ----setup, eval = FALSE------------------------------------------------------ # # Attempt to load a package, if the package was not available, install and load # if(!require("diemr", character.only = TRUE)){ # install.packages("diemr", dependencies = TRUE) # library("diemr", character.only = TRUE) # } ## ----eval = FALSE------------------------------------------------------------- # filepaths <- system.file("extdata", "data7x3.txt", package = "diemr") # # Analysing six individuals # samples <- 1:6 # # Assuming diploid markers of all individuals # ploidies = rep(list(rep(2, nchar(readLines(filepaths[1])[1]) - 1)), length(filepaths)) # CheckDiemFormat(files = filepaths, # ChosenInds = samples, # ploidy = ploidies) # # File check passed: TRUE # # Ploidy check passed: TRUE ## ----eval = FALSE------------------------------------------------------------- # res <- diem(files = filepaths, # ploidy = ploidies, # markerPolarity = list(c(FALSE, FALSE, TRUE)), # ChosenInds = samples, # nCores = 1) ## ----eval = FALSE------------------------------------------------------------- # res$DI # # newPolarity DI Support Marker # # 1 FALSE -4.872256 15.930181 1 # # 2 TRUE -3.520647 18.633399 2 # # 3 TRUE -13.274822 6.130625 3 ## ----eval = FALSE------------------------------------------------------------- # genotypes <- importPolarized(file = filepaths, # changePolarity = res$markerPolarity, # ChosenInds = samples) # # h <- apply(X = res$I4, # MARGIN = 1, # FUN = pHetErrOnStateCount)[1, ] # # plotPolarized(genotypes = genotypes, # HI = h[samples]) ## ----echo = FALSE, fig.width = 6, fig.height = 4.5, fig.align = 'center'------ genotypes <- diemr::importPolarized(file = system.file("extdata", "data7x3.txt", package = "diemr"), changePolarity = c(FALSE, TRUE, TRUE), ChosenInds = 1:6) h <- c(.5,0,.33,.5,.83,1) diemr::plotPolarized(genotypes = genotypes, HI = h) ## ----eval = FALSE------------------------------------------------------------- # S0011222 # S1210001 # S02221U0 ## ----eval = FALSE------------------------------------------------------------- # filepaths2 <- c(system.file("extdata", "data7x3.txt", package = "diemr"), # system.file("extdata", "data7x10.txt", package = "diemr")) # # ploidies2 <- list(rep(2, 7), # c(2, 1, 2, 2, 2, 1, 2)) # # CheckDiemFormat(files = filepaths2, # ChosenInds = samples, # ploidy = ploidies2) # # File check passed: TRUE # # Ploidy check passed: TRUE # # # Set random seed for repeatibility of null polarities (optional) # set.seed(39583782) # # # Run diem with verbose = TRUE to store hybrid indices with ploidy-aware allele counts # res2 <- diem(files = filepaths2, # ploidy = ploidies2, # markerPolarity = FALSE, # ChosenInds = samples, # nCores = 1, # verbose = TRUE) ## ----eval = FALSE------------------------------------------------------------- # # Import polarized genotypes for all compartments # genotypes2 <- importPolarized(files = filepaths2, # changePolarity = res2$markerPolarity, # ChosenInds = samples) # # # Load individual hybrid indices from a stored file # h2 <- unlist(read.table("diagnostics/HIwithOptimalPolarities.txt")) # # # Plot the polarised genotypes # plotPolarized(genotypes = genotypes2, # HI = h2[samples]) ## ----echo = FALSE, warning = FALSE,fig.width = 6, fig.height = 4.5, fig.align = 'center'---- genotypes2 <- diemr::importPolarized(files = c(system.file("extdata", "data7x3.txt", package = "diemr"), system.file("extdata", "data7x10.txt", package = "diemr")), changePolarity = c(FALSE, TRUE, TRUE,T,T,T,T,F,T,T,F,T,T), ChosenInds = 1:6) h2 <- apply(genotypes2, 1, \(x) diemr::pHetErrOnStateCount(diemr::sStateCount(x)))[1, ] diemr::plotPolarized(genotypes = genotypes2, HI = h2, lwd = 2) ## ----eval = FALSE------------------------------------------------------------- # # Select a threshold for the top diagnostic markers. # threshold <- quantile(res2$DI$DI, prob = 0.6) # # Create a vector identifying the diagnostic markers at the given threshold # markers <- res2$DI$DI > threshold # # Import only the selected markers # genotypes3 <- importPolarized(files = filepaths2, # changePolarity = res2$markerPolarity, # ChosenInds = samples, # ChosenSites = markers) # # Calculate hybrid index from diagnostic markers # h3 <- apply(genotypes3, 1, FUN = function(x) pHetErrOnStateCount(sStateCount(x)))[1, ] # # Plot diagnostic markers # plotPolarized(genotypes3, h3) ## ----echo = FALSE, warning = FALSE,fig.width = 6, fig.height = 4.5, fig.align = 'center'---- diemr::plotPolarized(genotypes = genotypes2[, c(TRUE,TRUE,FALSE,TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE)], HI = c(.3,0,.5,.3,.4,1), lwd = 2) ## ----eval = FALSE------------------------------------------------------------- # install.packages(package = "diemr_1.4.tar.gz", # repos = NULL, type = "source") ## ----eval = FALSE------------------------------------------------------------- # apply(res$I4, MARGIN = 1, FUN = pHetErrOnStateCount) # # [,1] [,2] [,3] [,4] [,5] [,6] # # p 0.5000000 0 0.3333333 0.5000000 0.8333333 1.0000000 # # Het 0.3333333 0 0.6666667 0.3333333 0.3333333 0.0000000 # # Err 0.0000000 0 0.0000000 0.0000000 0.0000000 0.3333333 ## ----eval = FALSE------------------------------------------------------------- # samples2 <- c(1, 3:6) # CheckDiemFormat(files = filepaths, # ChosenInds = samples2, # ploidy = ploidies) # # File check passed: TRUE # # Ploidy check passed: TRUE # # res2 <- diem(files = filepaths, # ChosenInds = samples2, # ploidy = ploidies, # nCores = 1, # markerPolarity = list(c(FALSE, FALSE, TRUE))) # # # calculate hybrid indices from updated I4 # h.res2 <- apply(res2$I4, # MARGIN = 1, # FUN = pHetErrOnStateCount)[1, ] # # # set names for the hybrid index values # h.res2 <- setNames(h.res2, nm = samples2) # # 1 3 4 5 6 # # 0.50 0.33 0.50 0.83 1.00 # # # calculate the center of the maximum hybrid index change # diffs <- data.frame(rollmean = zoo::rollmean(sort(h.res2), k = 2), # diff = diff(sort(h.res2), lag = 1)) # h.res2.c <- diffs$rollmean[which.max(diffs$diff)] # # [1] 0.6666667 ## ----eval = FALSE------------------------------------------------------------- # # Select 40% of markers with the highest diagnostic index # threshold <- quantile(res2$DI$DI, prob = 0.6) # genotypes3 <- genotypes2[, res2$DI$DI > threshold] # # Recalculate I4 and hybrid indices # h3 <- apply(genotypes3, # MARGIN = 1, # FUN = \(x) pHetErrOnStateCount(sStateCount(x)))[1, ] # # Plot the polarised markers # plotPolarized(genotypes3, h3) ## ----echo = FALSE------------------------------------------------------------- genotypes3 <- genotypes2[, c(1,2,4,10,11)] h3 <- apply(genotypes3, MARGIN = 1, FUN = \(x) diemr::pHetErrOnStateCount(diemr::sStateCount(x)))[1, ] diemr::plotPolarized(genotypes3, h3)