Development

.Rbuildignore

@@ -1,4 +1,15 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
-.github
-BIGf90_0.3.0.pdf
+^\.github$
+^\.git$
+^\.gitignore$
+^\.Rhistory$
+^\.RData$
+^\.DS_Store$
+^.*\.DS_Store$
+^BIGf90.*\.pdf$
+^inst/software
+^inst/exec_mac
+^inst/exec_windows
+^inst/data
+^inst/test_data_salmon

.gitignore

@@ -1,9 +1,21 @@
-.DS_store
-.Rproj.user
-.Rbuildignore
-BIGf90.Rproj
-inst/exec_mac
-inst/exec_windows/
-inst/data
-tests/
-.Rhistory
+# ---- OS / editor cruft ----
+.DS_Store
+**/.DS_Store
+.Rhistory
+.RData
+.Rproj.user
+
+# ---- RStudio project file ----
+BIGf90.Rproj
+
+# ---- Bundled BLUPf90 executables: keep locally, never push (platform-specific, ~165 MB) ----
+inst/software/
+inst/exec_mac/
+inst/exec_windows/
+inst/data/
+
+# ---- CVA/validation run outputs & local test data: keep locally, never push (~427 MB) ----
+inst/test_data_salmon/
+
+# NOTE: .Rbuildignore and tests/ are intentionally NOT ignored so they stay tracked
+#       (.Rbuildignore controls the built package; tests/ will hold testthat unit tests).

DESCRIPTION

@@ -1,22 +1,25 @@
 Package: BIGf90
-Title: R face front for running K-fold crossvalidation, estimating ebvs and variance component estimation with Blupf90 modules
-Version: 0.3.2
+Title: R Face Front for Running K-Fold Crossvalidation, Estimating EBVs and Variance Component Estimation with Blupf90 Modules
+Version: 0.4.0
 Authors@R: c(person("Josue", "Chinchilla-Vargas", role = c("aut", "cre"), email = "jc3635@cornell.edu"),
              person("Cristiane", "Taniguti", role = "aut"),
              person("Alexander", "Sandercock", role = "aut"),
              person("Breeding Insight Team", role = "ctb"))
 Maintainer: Josue Chinchilla-Vargas <jc3635@cornell.edu>
-Description: This package provides R functions to run several BLUPf90 modules. Along with an R function to run K-fold crossvalidation for univariate analyses through blupf90+. 
-    The output table from the K-fold crossvalidation function calculates accuracy as cor(y,ebv^) and cor(y*,ebv^) along with bias of ebvs calculated as reg(y,ebv.)
-    You will need to create the .par file to feed into Renumf90 manually. Once this is done, there are functions to run Renumf90, Blupf90+, Predictf90, Gibbsf90+ and Postgibbsf90.
+Description: Provides R functions to run several BLUPf90 modules. Along with an R function to run K-fold crossvalidation for univariate analyses through blupf90+. 
+    The output table from the K-fold crossvalidation function reports predictive ability as cor(y*, ebv^) and bias as the regression of y* on ebv^, with accuracy as cor(y*, ebv^)/sqrt(h2), where y* is the phenotype corrected for all fixed effects.
+    You can create the .par file to feed into Renumf90 manually or with write_par(). Once this is done, there are functions to run Renumf90, Blupf90+, Predictf90, Gibbsf90+ and Postgibbsf90.
     If you need to learn how to use the blupf90 suite of programs refer to nce.ads.uga.edu/wiki/doku.php?id=start. Please remember to cite Blupf90 appropriately along with this package  when used for publications.
     To run the functions in this package you will need to have a directory with all the blupf90 executables so that you can indicate the path for R to find the execs.
     As of version 0.3.0, all functions have been tested in PC and Unix environments, we have noticed that PC environments tend to have issues about permissions, so it is advised to run RStudio as administrator when planning to use this package.
 License: Apache License (== 2.0)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
-Imports: 
-    base (>= 4.3.1),
+Imports:
     dplyr (>= 1.1.4),
+    stats,
     utils (>= 4.3.1)
+Suggests:
+    testthat (>= 3.0.0)
+Config/testthat/edition: 3
+Config/roxygen2/version: 8.0.0

NAMESPACE

@@ -7,4 +7,5 @@ export(run_gibbs)
 export(run_postgibbs)
 export(run_predict)
 export(run_renum)
+export(write_par)
 import(dplyr)

NEWS.md

@@ -1,3 +1,14 @@
+# BIGf90 0.4.0
+
+* New function write_par() to build the raw RENUMF90 parameter (.par) file from R instead of writing it by hand
+* Effects are given as a named list keyed by a label; each effect gives its data-file column (col, or pos as an alias), type (cross/cov) and class (numer/alpha); an optional comment overrides the EFFECT #... label
+* random marks the random effect by its column (or label); the random effect is written last automatically so RANDOM follows it, with the pedigree FILE / FILE_POS block after
+* Supports genomic models (snp_file), permanent-environment / repeatability terms (optional = "pe") and (CO)VARIANCES starting values
+* Multi-trait models: give a trait vector with matrix RESIDUAL_VARIANCE / (CO)VARIANCES, and effect columns repeat per trait
+* Reproduces the sealice test .par files (pedigree and genomic)
+* Added a testthat suite covering write_par() and the create_folds() helper
+
+
 # BIGf90 0.3.1
 
 * Users now can have executable files, input files, and output files in different locations

R/bf90_cv.R

@@ -4,16 +4,14 @@
 #'
 #' This function sets up and runs a K-fold cross-validation analysis (CVA) using blupf90+ and predictf90.
 #' The function run_renumf90 needs to be used beforehand to process a .par file created by the user.
-#' This function calculates 2 accuracy estimates: correlation between raw phenotypes and ebvs divided by the square-root of narrow sense heritabilty and correlation between corrected phenotypes and ebvs along with bias estimations calculated as the regression of the phenotypes on the ebvs.
+#' Using the phenotype corrected for all fixed effects (y*), the function reports predictive ability as the correlation between y* and the ebvs, accuracy as that correlation divided by the square-root of the narrow-sense heritability, and bias as the regression of y* on the ebvs.
 #'
 #' @param missing_value_code code used in the .par file after OPTION MISSING to indicate missing phenotype, if this option is no use, this value must be 0.
 #' @param random_effect_col Column where random effect is located, found under RANDOM_GROUP in the renf90.par file.
 #' @param h2 estimate of narrow-sense heritabilty.This value is use to calculate accuracy of ebvs
 #' @param num_runs Number of independent cross-validation runs to be performed.
 #' @param num_folds Number of folds to be generated within each independent run.
 #' @param output_table_name Name of the final tab-separated out-up file. This field should be in quotes "".
-#' @param renf90_ped_name Name of pedigree file generated by renumf90. This field should be in quotes "".
-#' @param snp_file_name the name of the genotype file to be used. If use, this field should be in quotes"". Default for the function is no genotype file.
 #' @param path_2_execs path to a folder that holds all blupf90 executables that ill be used (blupf90+,predictf90). This field should be in quotes "".
 #' @param input_files_dir directory containing files renf90.par renf90.fields renf90.inb renf90.tables renf90.dat
 #' @param output_files_dir path to directory to store the output files
@@ -33,8 +31,8 @@
 #' #      num_runs = 5,
 #' #      num_folds = 10,
 #' #      output_table_name = "example_run",
-#' #     renf90_ped_name = "renadd03.ped",
-#' #     snp_file_name = "my_genos.geno" )
+#' #      input_files_dir = ".",
+#' #      output_files_dir = ".")
 #'
 #'
 #' @export
@@ -47,21 +45,19 @@ bf90_cv <- function(missing_value_code = NULL,
                     input_files_dir = ".",
                     output_files_dir = ".",
                     output_table_name = NULL,
-                    renf90_ped_name= NULL,
-                    snp_file_name = NULL,
                     seed = 101919,
                     verbose = TRUE) {
 
+  input_files_dir <- normalizePath(input_files_dir)
 
-  #check id pedigree and genotype file are present
-  if(is.null(renf90_ped_name)) stop("Specify renf90_ped_name")
-
-  if (!file.exists(renf90_ped_name)) {
-    stop("Parameter file not found at: ", renf90_ped_name)
+  # Checks
+  output_files_dir <- normalizePath(output_files_dir)
+  if(file.exists(output_files_dir)){
+    check_files <- list.files(output_files_dir)
+    if(length(check_files) > 0) warning(paste("Directory", output_files_dir, "is not empty. Some files may be replaced."))
+  } else {
+    stop(paste("Directory", output_files_dir, "does not exist. Create it before running the function."))
   }
-  renf90_ped_name <- normalizePath(renf90_ped_name)
-
-  input_files_dir <- normalizePath(input_files_dir)
 
   if(is.null(missing_value_code)) stop("Specify missing_value_code")
   if(is.null(random_effect_col)) stop("Specify random_effect_col")
@@ -70,36 +66,42 @@ bf90_cv <- function(missing_value_code = NULL,
   if(is.null(num_folds)) stop("Specify num_folds")
 
   if (!file.exists(file.path(input_files_dir,"renf90.par"))) stop("File 'renf90.par' not found at: ", input_files_dir)
-  if (!file.exists(file.path(input_files_dir,"renf90.inb"))) stop("File 'renf90.inb' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.fields"))) stop("File 'renf90.fields' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.tables"))) stop("File 'renf90.tables' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.dat"))) stop("File 'renf90.dat' not found at: ", input_files_dir)
 
-  if (!is.null(snp_file_name)) {
-    if (!file.exists(snp_file_name)) {
-      stop("Parameter file not found at:", snp_file_name)
-      snp_file_name <- normalizePath(snp_file_name)
-    }
+  # Check files in the parameter files
+  parfile <- readLines(file.path(input_files_dir,"renf90.par"))
+
+  snpfile <- grep("SNP_file", parfile)
+  if(length(snpfile) != 0) {
+    snp_file_name <- sapply(strsplit(parfile[snpfile], " "), function(x) x[length(x)])
+    if(!file.exists(file.path(input_files_dir, snp_file_name))) stop(paste("File", file.path(input_files_dir, snp_file_name),
+                                                                           "especified in parameter file",file.path(input_files_dir,"renf90.par"),
+                                                                           "line", snpfile, "does not exist."))
+    snp_file_name <- normalizePath(file.path(input_files_dir, snp_file_name))
+  } else snp_file_name <- NULL
+
+  # cleaned-SNP and saved G-inverse filenames (built once on the full data, reused per fold)
+  if(!is.null(snp_file_name)){
+    snp_base  <- base::basename(snp_file_name)      # e.g. sealice_match.geno
+    clean_snp <- base::paste0(snp_base, "_clean")   # produced by saveCleanSNPs
   }
 
-  if(is.null(output_table_name)) stop("Define output table name.")
+  pedfile <- grep(" FILE", parfile) + 1
+  if(length(pedfile) != 0) {
+    renf90_ped_name <- gsub(" ", "", parfile[pedfile])
+    if(!file.exists(file.path(input_files_dir, renf90_ped_name))) stop(paste("File", file.path(input_files_dir, renf90_ped_name),
+                                                                             "especified in parameter file",file.path(input_files_dir,"renf90.par"),
+                                                                             "line", pedfile, "does not exist."))
+    renf90_ped_name <- normalizePath(file.path(input_files_dir, renf90_ped_name))
+  } else renf90_ped_name <- NULL
 
-  if (!file.exists(renf90_ped_name)) {
-    stop("Parameter file not found at: ", renf90_ped_name)
-  }
+  if(is.null(output_table_name)) stop("Define output table name.")
 
   # set seed for reproducibility
   base::set.seed(seed)
 
-  # Checks
-  output_files_dir <- normalizePath(output_files_dir)
-  if(file.exists(output_files_dir)){
-    check_files <- list.files(output_files_dir)
-    if(length(check_files) > 0) warning(paste("Directory", output_files_dir, "is not empty. Some files may be replaced."))
-  } else {
-    stop(paste("Directory", output_files_dir, "does not exist. Create it before running the function."))
-  }
-
   path_2_execs <- normalizePath(path_2_execs)
 
   # Inform parameters and directories set
@@ -130,19 +132,52 @@ bf90_cv <- function(missing_value_code = NULL,
 
   # Run BF90 programs for the whole dataset
   setwd(input_files_dir)
-  output <- execute_command(command = paste0(file.path(path_2_execs, blup)," ", "renf90.par"), logfile = "run_blup.log")
-  command_predict <- paste0(file.path(path_2_execs, predict), " ", paste0("renf90.par"))
+  # For genomic analyses, make sure the G-inverse (Gi) and cleaned SNPs are saved on
+  # this full-data run so every fold can reuse them (built once here).
+  blup_parfile <- "renf90.par"
+  if(!is.null(snp_file_name)){
+    blup_par <- base::readLines("renf90.par")
+    if(!any(base::grepl("saveGInverse",  blup_par))) blup_par <- c(blup_par, "OPTION saveGInverse")
+    if(!any(base::grepl("saveCleanSNPs", blup_par))) blup_par <- c(blup_par, "OPTION saveCleanSNPs")
+    base::writeLines(blup_par, "renf90_blup.par")
+    blup_parfile <- "renf90_blup.par"
+  }
+  output <- execute_command(command = paste0(file.path(path_2_execs, blup)," ", blup_parfile), logfile = "run_blup.log")
+
+  # predictf90 on a copy of renf90.par that adjusts the phenotype for all effects
+  # except the random (animal) effect -> yhat = corrected phenotype (y*)
+  predict_par <- c(base::readLines("renf90.par"), base::paste("OPTION include_effects", random_effect_col))
+  base::writeLines(predict_par, "renf90_predict.par")
+  command_predict <- paste0(file.path(path_2_execs, predict), " ", "renf90_predict.par")
   output <- execute_command(command = command_predict, logfile = "run_predict.log")
 
-  files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions")
-  for(i in 1:length(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
+  if(is.null(snp_file_name)){
+    files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions")
+  } else {
+    files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions",
+                   "freqdata.count", "freqdata.count.after.clean", "Gen_call_rate", "Gen_conflicts",
+                   "run_predict.log", "sum2pq",
+                   "Gi", clean_snp, base::paste0(clean_snp, "_XrefID"))   # reused by every fold
+  }
+
+  for(i in 1:length(files_res)) if(file.exists(files_res[i])) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
 
   # Prepare files for each BLUP run
   renf90 <- base::readLines(paste0("renf90.par"))
   #ped_file <- dirname(renf90_ped_name)
   fields_file <- base::readLines(paste0("renf90.fields"))
-  inb_file <- base::readLines(paste0("renf90.inb"))
   tables_file <- base::readLines(paste0("renf90.tables"))
+
+  # Genomic: point each fold at the G-inverse + cleaned SNPs built on the full data,
+  # so folds skip genotype QC and G construction (identical results, much faster).
+  if(!is.null(snp_file_name)){
+    renf90 <- base::gsub("saveGInverse",  "readGInverse",       renf90, fixed = TRUE)
+    renf90 <- base::gsub("saveCleanSNPs", "no_quality_control", renf90, fixed = TRUE)
+    if(!any(base::grepl("readGInverse",       renf90))) renf90 <- c(renf90, "OPTION readGInverse")
+    if(!any(base::grepl("no_quality_control", renf90))) renf90 <- c(renf90, "OPTION no_quality_control")
+    renf90 <- base::gsub(snp_base, clean_snp, renf90, fixed = TRUE)   # SNP_file -> cleaned set
+    reuse_files <- base::file.path(output_files_dir, c("Gi", clean_snp, base::paste0(clean_snp, "_XrefID")))
+  } else reuse_files <- NULL
 
   # Read and preprocess the phenotype data
   bf90_phenos <- utils::read.table(paste0("renf90.dat"), sep = " ", header = FALSE) %>%
@@ -153,15 +188,20 @@ bf90_cv <- function(missing_value_code = NULL,
                                   dplyr::select((random_effect_col + 1)))
 
   folds <- base::lapply(data_shuffled, function(x) create_folds(x, num_folds))
-  mutated_data <- base::lapply(folds, function(f) mutate_folds(bf90_phenos, f, num_folds,missing_value_code))
+  mutated_data <- base::lapply(folds, function(f) mutate_folds(bf90_phenos, f, num_folds, missing_value_code, random_effect_col + 1))  # id col = random_effect_col+1
 
   setwd(output_files_dir)
   for (run in 1:num_runs) {
     for (fold in 1:num_folds) {
       dir_path <- base::sprintf("run%d/fold%d", run, fold)
-      create_cv_datasets(run, fold, mutated_data[[run]][[fold]],
-                         dir_path, renf90, renf90_ped_name,
-                         input_files_dir, snp_file_name)
+      create_cv_datasets(run, 
+                         fold, 
+                         data_frame =  mutated_data[[run]][[fold]],
+                         dir_path, 
+                         renf90, 
+                         renf90_ped_name,
+                         input_files_dir,
+                         reuse_files)
     }
   }
 
@@ -178,7 +218,7 @@ bf90_cv <- function(missing_value_code = NULL,
       data_file <- base::sprintf("renf90_run%d_fold%d.dat", run, fold)
       masked_ids <- utils::read.table(data_file) %>%
         dplyr::filter(V1 == missing_value_code) %>%
-        dplyr::select(4) %>%
+        dplyr::select(random_effect_col + 1) %>%   # id column
         base::unlist()
 
       ebvs_for_cv <- utils::read.table("solutions", header = FALSE, skip = 1) %>%
@@ -195,7 +235,6 @@ bf90_cv <- function(missing_value_code = NULL,
 
   # Calculate correlations and bias
   corrected_phenos <- utils::read.table(paste0(output_files_dir, "/yhat_residual"), header = FALSE) %>% dplyr::select(1, 2)
-  raw_phenos <- utils::read.table(paste0(input_files_dir, "/renf90.dat")) %>% dplyr::select(4, 1)
 
   ystar_correlations <- base::numeric(num_runs)
   bias_list <- base::numeric(num_runs)
@@ -205,9 +244,7 @@ bf90_cv <- function(missing_value_code = NULL,
     ystar <- dplyr::inner_join(ebvs_for_cv_runs[[base::sprintf("ebvs_for_cv_run%d", i)]], corrected_phenos, by = c("V3" = "V1"))
     ystar_correlations[i] <- base::round(stats::cor(ystar$V4, ystar$V2), 3)
 
-    yraw <- dplyr::inner_join(ebvs_for_cv_runs[[base::sprintf("ebvs_for_cv_run%d", i)]], raw_phenos, by = c("V3" = "V4"))
-
-    model <- stats::lm(V1 ~ V4, data = yraw)
+    model <- stats::lm(V2 ~ V4, data = ystar)   # bias (LR b1): corrected phenotype (V2) regressed on EBV (V4)
     bias_list[i] <- base::round(stats::coefficients(model)["V4"], 3)
 
     accuracy_list[i] <- base::round(ystar_correlations[i] / base::sqrt(h2), 3)

R/run_blup.R

@@ -33,15 +33,15 @@ run_blup <- function(path_2_execs) {
   }
 
   # Check if executable exists
-  if (!file.exists(paste0(path_2_execs, blup))) {
-    stop("Executable not found at: ", paste0(path_2_execs, blup))
+  if (!file.exists(file.path(path_2_execs, blup))) {
+    stop("Executable not found at: ", file.path(path_2_execs, blup))
   }
   if (!file.exists("renf90.par")) {
     stop("Parameter file not found: renf90.par")
   }
 
   # Run the command
-  output <- execute_command(command = paste0(path_2_execs, blup, " renf90.par"), logfile = "run_blup.log")
+  output <- execute_command(command = paste0(file.path(path_2_execs, blup), " renf90.par"), logfile = "run_blup.log")
 
   # Capture and print the log file content
   if (file.exists("run_blup.log")) {

R/run_postgibbs.R

@@ -77,14 +77,20 @@ run_postgibbs <- function(path_2_execs,
   # Remove the temporary file
   unlink(temp_input_file)
 
+  # Move result files to output directory if different from input
+  if(output_files_dir != input_files_dir){
+    files_res <- c(base::list.files(pattern = "^post"), "run_postgibbs.log")
+    for(i in seq_along(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir, files_res[i]))
+  }
+
   # Debugging: Print the output
   if(verbose) cat("Command output:", output, "\n")
 
   # Capture and print the log file content
   if(verbose){
-    if (file.exists("run_postgibbs.log")) {
+    if (file.exists(file.path(output_files_dir,"run_postgibbs.log"))) {
       cat("Log file content:\n")
-      cat(readLines("run_postgibbs.log"), sep = "\n")
+      cat(readLines(file.path(output_files_dir,"run_postgibbs.log")), sep = "\n")
     } else {
       cat("Log file not created.\n")
     }