diff --git a/.gitignore b/.gitignore
index 5b6a065..b64a99f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 .Ruserdata
+revdep/
diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION
index f52a032..ae8365d 100644
--- a/CRAN-SUBMISSION
+++ b/CRAN-SUBMISSION
@@ -1,3 +1,3 @@
-Version: 0.5.5
-Date: 2025-05-14 16:17:56 UTC
-SHA: 067acfa2a94f0baee7a151efd10fb13d5d50509b
+Version: 0.6.2
+Date: 2025-09-18 12:16:11 UTC
+SHA: 142dc9524d88b47db88ddca2aa39cd729a8d5a0d
diff --git a/DESCRIPTION b/DESCRIPTION
index 52220d3..31ad1e1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,9 +1,9 @@
 Package: BIGr
 Title: Breeding Insight Genomics Functions for Polyploid and Diploid Species
-Version: 0.5.5
+Version: 0.6.2
 Authors@R: c(person(given='Alexander M.',
                     family='Sandercock',
-                    email='ams866@cornell.edu',
+                    email='sandercock.alex@gmail.com',
                     role=c('cre','aut')),
              person(given='Cristiane',
                     family='Taniguti',
@@ -26,7 +26,7 @@ Authors@R: c(person(given='Alexander M.',
              person('Cornell', 'University',
                     role=c('cph'), 
                     comment = "Breeding Insight"))
-Maintainer: Alexander M. Sandercock <ams866@cornell.edu>
+Maintainer: Alexander M. Sandercock <sandercock.alex@gmail.com>
 Description: Functions developed within Breeding Insight to analyze
     diploid and polyploid breeding and genetic data. 'BIGr' provides the
     ability to filter variant call format (VCF) files, extract single nucleotide polymorphisms (SNPs) 
@@ -53,6 +53,7 @@ Imports:
     Rdpack (>= 0.7),
     readr (>= 2.1.5),
     reshape2 (>= 1.4.4),
+    rlang,
     tidyr (>= 1.3.1),
     vcfR (>= 1.15.0),
     Rsamtools,
@@ -60,7 +61,8 @@ Imports:
     pwalign,
     janitor,
     quadprog,
-    tibble
+    tibble,
+    stringr
 Suggests: 
     covr,
     spelling,
diff --git a/NAMESPACE b/NAMESPACE
index 2ca9c31..d47d95d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(allele_freq_poly)
 export(calculate_Het)
 export(calculate_MAF)
 export(check_homozygous_trios)
@@ -7,18 +8,24 @@ export(check_ped)
 export(check_replicates)
 export(dosage2vcf)
 export(dosage_ratios)
+export(filterMADC)
 export(filterVCF)
 export(flip_dosage)
 export(get_countsMADC)
 export(imputation_concordance)
+export(madc2gmat)
 export(madc2vcf_all)
 export(madc2vcf_targets)
 export(merge_MADCs)
+export(solve_composition_poly)
+export(thinSNP)
 export(updog2vcf)
 import(dplyr)
 import(janitor)
 import(parallel)
 import(quadprog)
+import(rlang)
+import(stringr)
 import(tibble)
 import(tidyr)
 import(vcfR)
diff --git a/NEWS.md b/NEWS.md
index 02cb34e..b089e67 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,25 +1,27 @@
-# BIGr 0.3.3
+# BIGr 0.6.2
 
--   Adapt updog2vcf to model f1, f1pp, s1 and s1pp
+- Fixed the doi and name list in the CITATION file
 
-# BIGr 0.3.2
+# BIGr 0.6.1
 
--   updog2vcf function option to output compressed VCF (.vcf.gz) - set as default
--   remove need for defining ploidy
--   add metadata at the VCF header
+- Added new functions for filtering MADC files and converting to relationship matrices
+- Added function thinSNPs to thin SNPs based on physical distance
+- Added bug fixes and improvements to existing functions
 
-# BIGr 0.5.0
+# BIGr 0.5.5
 
--   Add imputation_concordance function to estimate accuracy of imputed and original dataset
--   Add get_OffTargets function to extract target and off-target SNPs from a MADC file
--   Add merge_MADCs function to merge two or more MADC files together
--   Improved documentation and examples for all functions
--   Add tests for all functions
+- Updated DESCRIPTION
+- Added return value for merge_MADCs
+- Added optional seed for check_ped
+- Added verbose option
 
-# BIGr 0.5.1
+# BIGr 0.5.4
 
--   Improvements of testthat tests
--   Add check_replicates and check_homozygous_trios for pedigree relationship quality check
+-   Updated dosage2vcf example
+
+# BIGr 0.5.3
+
+-   Updated madc2vcf_all example
 
 # BIGr 0.5.2
 
@@ -27,17 +29,26 @@
 -   get_OffTargets function changed to madc2vcf_all
 -   Updates to testthat tests and function examples
 
-# BIGr 0.5.3
+# BIGr 0.5.1
 
--   Updated madc2vcf_all example
+-   Improvements of testthat tests
+-   Add check_replicates and check_homozygous_trios for pedigree relationship quality check
 
-# BIGr 0.5.4
+# BIGr 0.5.0
 
--   Updated dosage2vcf example
+-   Add imputation_concordance function to estimate accuracy of imputed and original dataset
+-   Add get_OffTargets function to extract target and off-target SNPs from a MADC file
+-   Add merge_MADCs function to merge two or more MADC files together
+-   Improved documentation and examples for all functions
+-   Add tests for all functions
 
-# BIGr 0.5.5
+# BIGr 0.3.3
+
+-   Adapt updog2vcf to model f1, f1pp, s1 and s1pp
+
+# BIGr 0.3.2
+
+-   updog2vcf function option to output compressed VCF (.vcf.gz) - set as default
+-   remove need for defining ploidy
+-   add metadata at the VCF header
 
-- Updated DESCRIPTION
-- Added return value for merge_MADCs
-- Added optional seed for check_ped
-- Added verbose option
diff --git a/R/breedtools_functions.R b/R/breedtools_functions.R
index db99857..5f21836 100644
--- a/R/breedtools_functions.R
+++ b/R/breedtools_functions.R
@@ -32,7 +32,7 @@
 #' allele_freqs <- allele_freq_poly(geno = geno_matrix, populations = pop_list, ploidy = 4)
 #' print(allele_freqs)
 #'
-#' @noRd
+#' @export
 allele_freq_poly <- function(geno, populations, ploidy = 2) {
 
   # Initialize returned df
@@ -168,7 +168,7 @@ QPsolve <- function(Y, X) {
 #'                                       ploidy = 4)
 #' print(composition)
 #'
-#' @noRd
+#' @export
 solve_composition_poly <- function(Y,
                                    X,
                                    ped = NULL,
diff --git a/R/filterMADC.R b/R/filterMADC.R
new file mode 100644
index 0000000..5c1056f
--- /dev/null
+++ b/R/filterMADC.R
@@ -0,0 +1,215 @@
+#' Filter MADC Files
+#'
+#' Filter and process MADC files to remove low quality microhaplotypes
+#'
+#' @details
+#' This function can filter raw MADC files or pre-processed MADC files with fixed allele IDs. Additionally,
+#' it can filter based on mean read depth, number of mhaps per target loci, and other criteria. Optionally, users
+#' can plot summary statistics and save the filtered data to a file.
+#'
+#'@import dplyr
+#'@importFrom utils read.csv
+#'
+#'@param madc_file Path to the MADC file to be filtered
+#'@param min.mean.reads Minimum mean read depth for filtering
+#'@param max.mean.reads Maximum mean read depth for filtering
+#'@param max.mhaps.per.loci Maximum number of matching mhaps per target loci. Retains only the target Ref and Alt loci at the sites that exceeds the \code{max.mhaps.per.loci} threshold.
+#'@param min.reads.per.site Minimum number of reads per site for \code{min.ind.with.reads}. Otherwise, this parameter is ignored
+#'@param min.ind.with.reads Minimum number of individuals with \code{min.reads.per.site} reads for filtering
+#'@param target.only Logical indicating whether to filter for target loci only
+#'@param n.summary.columns (optional) Number of summary columns to remove from MADC file not including the first three. Otherwise, the columns will be automatically detected and removed.
+#@param plot.summary Logical indicating whether to plot summary statistics
+#'@param output.file Path to save the filtered data (if NULL, data will not be saved)
+#'
+#'@return data.frame or saved csv file
+#'
+#'@examples
+#' #Example
+#'
+#' #Example MADC
+#' madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+#'
+#' #Remove mhaps exceeding 3 per target region including the ref and alt target mhaps
+#' filtered_df <- filterMADC(madc_file,
+#'                          min.mean.reads = NULL,
+#'                          max.mean.reads = NULL,
+#'                          max.mhaps.per.loci = 3,
+#'                          min.reads.per.site = 1,
+#'                          min.ind.with.reads = NULL,
+#'                          target.only = FALSE,
+#'                          n.summary.columns = NULL,
+#'                          output.file = NULL)
+#'
+#'
+#'
+#'@export
+filterMADC <- function(madc_file,
+                       min.mean.reads = NULL,
+                       max.mean.reads = NULL,
+                       max.mhaps.per.loci = NULL,
+                       min.reads.per.site = 1,
+                       min.ind.with.reads = NULL,
+                       target.only = FALSE,
+                       n.summary.columns = NULL,
+                       #plot.summary = FALSE,
+                       output.file = NULL) {
+
+
+  #Need to first inspect the first 7 rows of the MADC to see if it has been preprocessed or not
+  first_seven_rows <- read.csv(madc_file, header = FALSE, nrows = 7, colClasses = c(NA, "NULL"))
+
+  #Check if all entries in the first column are either blank or "*"
+  check_entries <- all(first_seven_rows[, 1] %in% c("", "*"))
+
+  #Check if the MADC file has the filler rows or is processed from updated fixed allele ID pipeline
+  if (check_entries) {
+    #Note: This assumes that the first 7 rows are placeholder info from DArT processing
+
+    warning("The MADC file has not been pre-processed for Fixed Allele IDs. The first 7 rows are placeholder info from DArT processing.")
+
+    #Read the madc file
+    filtered_df <- read.csv(madc_file, sep = ',', skip = 7, check.names = FALSE)
+
+    #Remove extra text after Ref and Alt (_001 or _002)
+    #filtered_df$AlleleID <- sub("\\|Ref_.*", "|Ref", filtered_df$AlleleID)
+    #filtered_df$AlleleID <- sub("\\|Alt_.*", "|Alt", filtered_df$AlleleID)
+
+  } else {
+
+    #Read the madc file
+    filtered_df <- read.csv(madc_file, sep = ',', check.names = FALSE)
+
+    #Remove extra text after Ref and Alt (_001 or _002)
+    #filtered_df$AlleleID <- sub("\\|Ref_.*", "|Ref", filtered_df$AlleleID)
+    #filtered_df$AlleleID <- sub("\\|Alt_.*", "|Alt", filtered_df$AlleleID)
+
+  }
+  #Check for extra columns
+  #Save the three columns for later adding to the output
+  saved_columns <- filtered_df[,1:3]
+
+  if (!is.null(n.summary.columns)) {
+    #Remove the first n.summary.columns columns
+    filtered_df <- filtered_df[,-c(4:n.summary.columns)]
+  }else{
+    rm.col <- c("ClusterConsensusSequence",
+                "CallRate", "OneRatioRef", "OneRatioSnp", "FreqHomRef", "FreqHomSnp",
+                "FreqHets", "PICRef", "PICSnp", "AvgPIC", "AvgCountRef", "AvgCountSnp","RatioAvgCountRefAvgCountSnp")
+
+    filtered_df <- filtered_df[, !(colnames(filtered_df) %in% rm.col)]
+  }
+
+  #Now add rownames
+  rownames(filtered_df) <- saved_columns[,1]
+
+  #Remove refmatch and altmatch if wanted
+  if (target.only) {
+    message("Retaining target markers only")
+    #Retain only the Ref and Alt haplotypes
+    filtered_df <- filtered_df[!grepl("\\|AltMatch|\\|RefMatch", filtered_df$AlleleID), ]
+  }
+
+  ## Filtering
+
+  #Min mean reads
+  if (!is.null(min.mean.reads)) {
+    message("Filtering for minimum mean reads across all samples")
+    #Get the mean value for each row, and remove the rows below that threshold
+    filtered_df$MeanReads <- rowMeans(filtered_df[, -c(1:3)], na.rm = TRUE)
+    filtered_df <- filtered_df[filtered_df$MeanReads >= min.mean.reads, ]
+    #Remove the MeanReads column
+    filtered_df <- filtered_df[, -which(colnames(filtered_df) == "MeanReads")]
+  }
+
+  #Max mean reads
+  if (!is.null(max.mean.reads)) {
+    message("Filtering for maximum mean reads across all samples")
+    #Get the mean value for each row, and remove the rows above that threshold
+    filtered_df$MeanReads <- rowMeans(filtered_df[, -c(1:3)], na.rm = TRUE)
+    filtered_df <- filtered_df[filtered_df$MeanReads <= max.mean.reads, ]
+    #Remove the MeanReads column
+    filtered_df <- filtered_df[, -which(colnames(filtered_df) == "MeanReads")]
+  }
+
+  #Max mhaps per loci
+  if (!is.null(max.mhaps.per.loci)) {
+    message("Filtering for maximum number of matching mhaps per target loci")
+    #Get the number of matching mhaps for loci, and remove the mhaps at those loci that exceed the max number
+    mhap_counts <- filtered_df %>%
+      group_by(CloneID) %>%
+      summarise(Count = n(), .groups = 'drop') %>%
+      filter(Count > max.mhaps.per.loci)
+
+    patterns_to_search <- "\\|AltMatch|\\|RefMatch"
+    clone_ids_to_target <- mhap_counts$CloneID
+
+    filtered_df <- filtered_df %>%
+      filter(
+        !( # "keep rows that DO NOT match both conditions"
+          CloneID %in% clone_ids_to_target &  # Condition 1: CloneID is one of the targeted IDs
+            grepl(patterns_to_search, AlleleID) # Condition 2: AlleleID contains one of the patterns
+        )
+      )
+  }
+
+  #Min individuals with reads
+  if (!is.null(min.ind.with.reads)) {
+    message("Filtering for minimum number of individuals with reads per site")
+    message(paste0("Minimum number of individuals with reads per site: ", min.ind.with.reads))
+    message(paste0("Minimum number of reads per site: ", min.reads.per.site))
+
+    #Getting colnames
+    cols_to_check <- colnames(filtered_df)[-(1:3)]
+
+    filtered_df <- filtered_df %>%
+      rowwise() %>%  # Process data row by row
+      mutate(
+        # For each row, count how many of the 'cols_to_check' meet the criterion
+        qualifying_sites_count = sum(
+          c_across(all_of(cols_to_check)) >= min.reads.per.site,
+          na.rm = TRUE # Treats NAs in data as not meeting the criterion
+        )
+      ) %>%
+      ungroup() %>% # Always ungroup after rowwise operations
+      # Filter rows where this count meets the 'min.ind.with.reads' threshold
+      filter(qualifying_sites_count >= min.ind.with.reads) %>%
+      # Optionally, remove the temporary count column if it's no longer needed
+      select(-qualifying_sites_count)
+  }
+
+  #Plots
+  #if (plot.summary) {
+  #  message("Plotting summary statistics")
+  #  #Plot mean read depth
+  #  mean_reads <- rowMeans(filtered_df[, -c(1:3)], na.rm = TRUE)
+  #  hist(mean_reads, main = "Mean Read Depth", xlab = "Mean Reads", ylab = "Frequency")
+
+  #  #Plot number of Altmatch and Refmatch mhaps per target loci
+  #  altmatch_counts <- filtered_df %>%
+  #    filter(grepl("\\|AltMatch", AlleleID)) %>%
+  #    group_by(CloneID) %>%
+  #    summarise(Count = n(), .groups = 'drop')
+
+  #  refmatch_counts <- filtered_df %>%
+  #    filter(grepl("\\|RefMatch", AlleleID)) %>%
+  #    group_by(CloneID) %>%
+  #    summarise(Count = n(), .groups = 'drop')
+
+  #  barplot(cbind(altmatch_counts$Count, refmatch_counts$Count), beside = TRUE,
+  #          names.arg = altmatch_counts$CloneID, main = "Number of AltMatch and RefMatch Mhaps",
+  #          xlab = "Clone ID", ylab = "Count")
+
+    #Plot density of number of CloneID per site on a marker distribution plot
+
+  #}
+
+  #Save the output to disk if file name provided
+  if (!is.null(output.file)) {
+    message("Saving filtered data to file")
+    write.csv(filtered_df, paste0(output.file,".csv"), row.names = FALSE)
+  } else {
+    message("No output file provided. Returning filtered data.")
+    return(filtered_df)
+  }
+
+}
diff --git a/R/madc2gmat.R b/R/madc2gmat.R
new file mode 100644
index 0000000..44ccb0f
--- /dev/null
+++ b/R/madc2gmat.R
@@ -0,0 +1,217 @@
+#' Convert MADC Files to an Additive Genomic Relationship Matrix
+#'
+#' Scale and normalize MADC read count data and convert it to an additive genomic relationship matrix.
+#'
+#'@details
+#' This function reads a MADC file, processes it to remove unnecessary columns, and then
+#' converts it into an additive genomic relationship matrix using the first method proposed by VanRaden (2008).
+#' The resulting matrix can be used for genomic selection or other genetic analyses.
+#'
+#'@importFrom utils read.csv write.csv
+#'@import tibble stringr dplyr tidyr
+#'
+#'@param madc_file Path to the MADC file to be filtered
+#'@param seed Optional seed for random number generation (default is NULL)
+#'@param method Method to use for processing the MADC data. Options are "unique" or "collapsed". Default is "collapsed".
+#'@param ploidy Numeric. Ploidy level of the samples (e.g., 2 for diploid, 4 for tetraploid)
+#'@param output.file Path to save the filtered data (if NULL, data will not be saved)
+#'
+#'@return data.frame or saved csv file
+#'
+#'@examples
+#' #Input variables
+#' madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+#'
+#' #Calculations
+#' temp <- tempfile()
+#'
+#' # Converting to additive relationship matrix
+#' gmat <- madc2gmat(madc_file,
+#'                  seed = 123,
+#'                  ploidy=2,
+#'                  output.file = NULL)
+#'
+#'@references VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science, 91(11), 4414-4423
+#'@author Alexander M. Sandercock
+#'
+#'@export
+madc2gmat <- function(madc_file,
+                      seed = NULL,
+                      method = "collapsed",
+                      ploidy,
+                      output.file = NULL) {
+  #set seed if not null
+  if (!is.null(seed)) {
+    set.seed(seed)
+  }
+
+  #Need to first inspect the first 7 rows of the MADC to see if it has been preprocessed or not
+  first_seven_rows <- read.csv(madc_file, header = FALSE, nrows = 7, colClasses = c(NA, "NULL"))
+
+  #Check if all entries in the first column are either blank or "*"
+  check_entries <- all(first_seven_rows[, 1] %in% c("", "*"))
+
+  #Check if the MADC file has the filler rows or is processed from updated fixed allele ID pipeline
+  if (check_entries) {
+    #Note: This assumes that the first 7 rows are placeholder info from DArT processing
+
+    warning("The MADC file has not been pre-processed for Fixed Allele IDs. The first 7 rows are placeholder info from DArT processing.")
+
+    #Read the madc file
+    filtered_df <- read.csv(madc_file, sep = ',', skip = 7, check.names = FALSE)
+
+  } else {
+
+    #Read the madc file
+    filtered_df <- read.csv(madc_file, sep = ',', check.names = FALSE)
+  }
+
+  #Convert to data.table
+  #filtered_df <- as.data.table(filtered_df)
+
+  if (method == "unique") {
+    #Get allele ratios
+    # --- Step 1: Calculate frequency of each specific allele using total locus count ---
+    # The denominator MUST be the sum of all reads at the locus.
+    allele_freq_df <- filtered_df %>%
+      group_by(CloneID) %>%
+      mutate(across(
+        .cols = where(is.numeric),
+        # The calculation is count / total_locus_count
+        .fns = ~ . / sum(.),
+        .names = "{.col}_freq"
+      )) %>%
+      ungroup()
+
+    #Rm object
+    rm(filtered_df)
+
+    # --- Step 2: Prepare data for reshaping ---
+    # Filter out the reference rows and create a unique marker name for each alternative allele
+    markers_to_pivot <- allele_freq_df %>%
+      filter(!str_ends(AlleleID, "\\|Ref_0001")) %>%
+      # Create a new, unique column for each marker (e.g., "chr1.1_000194324_RefMatch_0001")
+      # This makes each alternative allele its own column in the final matrix.
+      mutate(MarkerID = str_replace(AlleleID, "\\|", "_")) %>%
+      select(MarkerID, ends_with("_freq")) %>%
+      tibble::column_to_rownames(var = "MarkerID")
+    names(markers_to_pivot) <- str_remove(names(markers_to_pivot), "_freq$")
+    markers_to_pivot <- t(markers_to_pivot)
+    #Replace the NaN to NA
+    markers_to_pivot[is.nan(markers_to_pivot)] <- NA
+
+    #rm unneeded objects
+    rm(allele_freq_df)
+
+  } else if (method == "collapsed"){
+
+    # This single pipeline calculates the final matrix.
+    # The output is the marker dosage matrix (X) to be used as input for A.mat()
+    markers_to_pivot <- filtered_df %>%
+
+      # --- Part A: Calculate the frequency of each allele at each locus ---
+      group_by(CloneID) %>%
+      mutate(across(
+        .cols = where(is.numeric),
+        # Use if_else to prevent 0/0, which results in NaN
+        .fns = ~ . / sum(.),
+        .names = "{.col}_freq"
+      )) %>%
+      ungroup() %>%
+
+      # --- Part B: Isolate and sum all alternative allele frequencies ---
+      #Note, the RefMatch counts are being counted as Reference allele counts, and not
+      #Alternative allele counts
+      # Filter to keep only the alternative alleles using base R's endsWith()
+      filter(!grepl("\\|Ref", AlleleID)) %>%
+      # Group by the main locus ID
+      group_by(CloneID) %>%
+      # For each locus, sum the frequencies of all its alternative alleles
+      summarise(across(ends_with("_freq"), sum), .groups = 'drop') %>%
+
+      # --- Part C: Reshape the data into the final matrix format ---
+      # Convert from a "long" format to a "wide" format
+      pivot_longer(
+        cols = -CloneID,
+        names_to = "SampleID",
+        values_to = "Dosage"
+      ) %>%
+      # Clean up the sample names using base R's sub()
+      mutate(SampleID = sub("_freq$", "", SampleID)) %>%
+      # Pivot to the final shape: samples in rows, markers in columns
+      pivot_wider(
+        names_from = CloneID,
+        values_from = Dosage,
+        # This ensures that if a sample/locus combination is missing,
+        # it gets a value of 0 instead of NA.
+        values_fill = 0
+      ) %>%
+      # Convert the 'SampleID' column into the actual row names of the dataframe
+      column_to_rownames("SampleID") %>%
+      # Convert the final dataframe into a true matrix object
+      as.matrix()
+    #Replace the NaN to NA
+    markers_to_pivot[is.nan(markers_to_pivot)] <- NA
+    #rm unneeded objects
+    rm(filtered_df)
+
+  } else {
+    stop("Invalid method specified. Use 'unique' or 'collapsed'.")
+  }
+
+  #VanRaden method
+  #This method is used to calculate the additive relationship matrix
+  compute_relationship_matrix <- function(mat,
+                                          ploidy,
+                                          method = "VanRaden",
+                                          impute = TRUE) {
+    ### Prepare matrix
+    #Remove any SNPs that are all NAs
+    mat <- mat[, colSums(is.na(mat)) < nrow(mat)]
+    #impute the missing values with the mean of the column
+    if (impute) {
+      mat <- data.frame(mat, check.names = FALSE, check.rows = FALSE) %>%
+        mutate_all(~ifelse(is.na(.x), mean(.x, na.rm = TRUE), .x)) %>%
+        as.matrix()
+    }
+
+    if (method == "VanRaden") {
+      # Calculate the additive relationship matrix using VanRaden's method
+
+      #Calculate allele frequencies
+      p <- apply(mat, 2, mean, na.rm = TRUE)
+      q <- 1 - p
+      #Calculate the denominator
+      #Note: need to use 1/ploidy for ratios, where ploidy should be used for dosages.
+      #This is because the ratios are from 0-1, which is what you get when dosage/ploidy, whereas dosages are from 0 to ploidy
+      denominator <- (1/as.numeric(ploidy))*sum(p * q, na.rm = TRUE)
+      #Get the numerator
+      Z <- scale(mat, center = TRUE, scale = FALSE)
+      ZZ <- (Z %*% t(Z))
+      #Calculate the additive relationship matrix
+      A.mat <- ZZ / denominator
+
+      return(A.mat)
+
+    } else {
+      stop("Unsupported method. Only 'VanRaden' is currently implemented.")
+    }
+  }
+  #Calculate the additive relationship matrix
+  MADC.mat <- compute_relationship_matrix(markers_to_pivot,
+                                          ploidy = ploidy,
+                                          method = "VanRaden",
+                                          impute = TRUE)
+
+  #Remove the markers_to_pivot object to free up memory
+  rm(markers_to_pivot)
+
+  #Save the output to disk if file name provided
+  if (!is.null(output.file)) {
+    message("Saving filtered data to file")
+    write.csv(MADC.mat, paste0(output.file,".csv"), row.names = TRUE)
+  } else {
+    return(MADC.mat)
+  }
+
+}
diff --git a/R/madc2vcf_all.R b/R/madc2vcf_all.R
index 13d1bfe..329ac5a 100644
--- a/R/madc2vcf_all.R
+++ b/R/madc2vcf_all.R
@@ -75,7 +75,7 @@ madc2vcf_all <- function(madc = NULL,
                       "out_vcf= ", out_vcf, ', ',
                       "verbose= ", verbose,')">')
 
-  if(!is.null(madc)) report <- read.csv(madc) else stop("Please provide a MADC file")
+  if(!is.null(madc)) report <- read.csv(madc, check.names = FALSE) else stop("Please provide a MADC file")
   if(!is.null(botloci_file)) botloci <- read.csv(botloci_file, header = F) else stop("Please provide a botloci file")
   if(!is.null(hap_seq_file)) hap_seq <- read.table(hap_seq_file, header = F) else hap_seq <- NULL
 
@@ -142,12 +142,12 @@ loop_though_dartag_report <- function(report, botloci, hap_seq, n.cores=1, align
   new.file <- data.frame("AlleleID" = report[,1],
                          "Chromosome" = NA, "SNP_position_in_Genome" = NA,
                          "Ref" = NA, "Alt" =NA,
-                         "CloneID" = report[,2], report[,3:ncol(report)])
+                         "CloneID" = report[,2], report[,3:ncol(report)], check.names = FALSE)
 
   by_cloneID <- split.data.frame(new.file, new.file$CloneID)
 
   clust <- makeCluster(n.cores)
-  #clusterExport(clust, c("hap_seq","add_ref_alt"))
+  #clusterExport(clust, c("hap_seq","add_ref_alt", "nsamples"))
   add_ref_alt_results <- parLapply(clust, by_cloneID, function(x) add_ref_alt(x, hap_seq, nsamples, verbose = verbose))
   stopCluster(clust)
 
@@ -169,7 +169,7 @@ loop_though_dartag_report <- function(report, botloci, hap_seq, n.cores=1, align
 
   clust <- makeCluster(n.cores)
   #clusterExport(clust, c("botloci", "compare", "nucleotideSubstitutionMatrix", "pairwiseAlignment", "DNAString", "reverseComplement"))
-  #clusterExport(clust, c("botloci"))
+  #clusterExport(clust, c("botloci", "alignment_score_thr"))
   compare_results <- parLapply(clust, updated_by_cloneID, function(x) compare(x, botloci, alignment_score_thr))
   stopCluster(clust)
 
@@ -209,7 +209,6 @@ loop_though_dartag_report <- function(report, botloci, hap_seq, n.cores=1, align
 #'
 #' @noRd
 add_ref_alt <- function(one_tag, hap_seq, nsamples, verbose = TRUE) {
-
   # Add ref and alt
   cloneID <- one_tag$CloneID[1]
   ref <- paste0(cloneID, "|Ref_0001")
@@ -281,13 +280,14 @@ add_ref_alt <- function(one_tag, hap_seq, nsamples, verbose = TRUE) {
 #'
 #' @noRd
 compare <- function(one_tag, botloci, alignment_score_thr = 40){
-  #one_tag <- updated_by_cloneID[[2]]
   cloneID <- one_tag$CloneID[1]
   isBotLoci <- cloneID %in% botloci[,1]
   # If marker is present in the botloc list, get the reverse complement sequence
   if(isBotLoci) one_tag$AlleleSequence <- sapply(one_tag$AlleleSequence, function(x) as.character(reverseComplement(DNAString(x))))
 
   chr <- sapply(strsplit(cloneID, "_"), function(x) x[-length(x)])
+  if(length(chr > 1)) chr <- paste(chr, collapse = "_")
+
   # Target SNP at the position in the ID
   ref_seq <- one_tag$AlleleSequence[grep("Ref_0001$",one_tag$AlleleID)]
   alt_seq <- one_tag$AlleleSequence[grep("Alt_0002$",one_tag$AlleleID)]
@@ -305,9 +305,9 @@ compare <- function(one_tag, botloci, alignment_score_thr = 40){
     alt_base <- substring(alt_seq, align@subject@mismatch@unlistData, align@subject@mismatch@unlistData)
 
     # If target alternative have N, discard whole tag
-    if(all(alt_base %in% c("A", "T", "C", "G"))) {
+    # Always only one polymorphism, if there are more than one, not sure which is the target
+    if(all(alt_base %in% c("A", "T", "C", "G")) & length(align@pattern@mismatch@unlistData) == 1) {
 
-      # Update with new info - always only one polymorphism, if there are more than one, not sure which is the target
       update_tag <- one_tag[grep("Ref_0001$",one_tag$AlleleID),]
       update_tag[,2:5] <- c(chr, pos_target, ref_base, alt_base)
       update_tag_temp <- one_tag[grep("Alt_0002$",one_tag$AlleleID),]
@@ -332,8 +332,9 @@ compare <- function(one_tag, botloci, alignment_score_thr = 40){
 
             # If Match sequences have N, do not consider as polymorphism
             if(any(!alt_base %in% c("A", "T", "C", "G"))) {
-              alt_base <- alt_base[-which(!alt_base %in% c("A", "T", "C", "G"))]
               ref_base <- ref_base[-which(!alt_base %in% c("A", "T", "C", "G"))]
+              pos_ref_idx <- pos_ref_idx[-which(!alt_base %in% c("A", "T", "C", "G"))]
+              alt_base <- alt_base[-which(!alt_base %in% c("A", "T", "C", "G"))]
             }
 
             if(length(alt_base) >0){ # If the N is the only polymorphis found, the Match tag will be discarted
@@ -366,6 +367,7 @@ compare <- function(one_tag, botloci, alignment_score_thr = 40){
                 rm_score = cloneID,
                 rm_N = NULL))
   }
+
 }
 
 #' Converts the fasta to a data.frame with first column the AlleleID and second the AlleleSequence
@@ -449,7 +451,11 @@ create_VCF_body <- function(csv,
   vcf_body$V3 <- as.numeric(vcf_body$V3)
   rownames(vcf_body) <- NULL
 
-  vcf_body$target <- paste0(sapply(strsplit(vcf_body$target, "_"),"[[", 1), "_",as.numeric(sapply(strsplit(vcf_body$target, "_"),"[[", 2)))
+  # Remove padding
+  sp <- strsplit(vcf_body$target, "_")
+  pos <- sapply(sp, function(x) x[length(x)])
+  chr <- sapply(sp, function(x) paste0(x[-length(x)], collapse = "_"))
+  vcf_body$target <- paste0(chr, "_",as.numeric(pos))
 
   # Dealing with repeated positions
   # discard the ones that are not the target and keep only the first if all are off-targets
diff --git a/R/madc2vcf_targets.R b/R/madc2vcf_targets.R
index 6cd61e0..1b02c31 100644
--- a/R/madc2vcf_targets.R
+++ b/R/madc2vcf_targets.R
@@ -97,6 +97,10 @@ madc2vcf_targets <- function(madc_file, output.file, botloci_file, get_REF_ALT =
     botloci <- checked_botloci[[1]]
     csv <- checked_botloci[[2]]
 
+    # FIXED: Store original sequences before any transformation
+    csv$OriginalAlleleSequence <- csv$AlleleSequence
+    
+    # Apply reverse complement to sequences for bottom strand markers
     idx <- which(csv$CloneID %in% botloci[,1])
     csv$AlleleSequence[idx] <- sapply(csv$AlleleSequence[idx], function(sequence) as.character(reverseComplement(DNAString(sequence))))
 
@@ -104,21 +108,44 @@ madc2vcf_targets <- function(madc_file, output.file, botloci_file, get_REF_ALT =
     ref_ord <- csv$CloneID[grep("\\|Ref.*", csv$AlleleID)]
     alt_seq <- csv$AlleleSequence[grep("\\|Alt.*", csv$AlleleID)]
     alt_ord <- csv$CloneID[grep("\\|Alt.*", csv$AlleleID)]
+    
+    # FIXED: Get original sequences for SNP calling
+    orig_ref_seq <- csv$OriginalAlleleSequence[grep("\\|Ref.*", csv$AlleleID)]
+    orig_alt_seq <- csv$OriginalAlleleSequence[grep("\\|Alt.*", csv$AlleleID)]
 
     if(all(sort(ref_ord) == sort(alt_ord))){
+      # Order sequences consistently
       ref_seq <- ref_seq[order(ref_ord)]
       alt_seq <- alt_seq[order(alt_ord)]
+      orig_ref_seq <- orig_ref_seq[order(ref_ord)]
+      orig_alt_seq <- orig_alt_seq[order(alt_ord)]
+      ordered_clone_ids <- sort(ref_ord)
 
       ref_base <- alt_base <- vector()
-      for(i in 1:length(ref_seq)){
-        temp_list <- strsplit(c(ref_seq[i], alt_seq[i]), "")
-        idx <- which(temp_list[[1]] != temp_list[[2]])
-        if(length(idx) >1) { # If finds more than one polymorphism between Ref and Alt sequences
+      for(i in 1:length(orig_ref_seq)){
+        # FIXED: Use original sequences for SNP calling
+        temp_list <- strsplit(c(orig_ref_seq[i], orig_alt_seq[i]), "")
+        idx_diff <- which(temp_list[[1]] != temp_list[[2]])
+        
+        if(length(idx_diff) > 1) { # If finds more than one polymorphism between Ref and Alt sequences
           ref_base[i] <- NA
           alt_base[i] <- NA
+        } else if(length(idx_diff) == 1) {
+          orig_ref_base <- temp_list[[1]][idx_diff]
+          orig_alt_base <- temp_list[[2]][idx_diff]
+          
+          # FIXED: Apply reverse complement to bases only if marker is in botloci
+          if(ordered_clone_ids[i] %in% botloci[,1]) {
+            ref_base[i] <- as.character(reverseComplement(DNAString(orig_ref_base)))
+            alt_base[i] <- as.character(reverseComplement(DNAString(orig_alt_base)))
+          } else {
+            ref_base[i] <- orig_ref_base
+            alt_base[i] <- orig_alt_base
+          }
         } else {
-          ref_base[i] <- temp_list[[1]][idx]
-          alt_base[i] <- temp_list[[2]][idx]
+          # No differences found
+          ref_base[i] <- NA
+          alt_base[i] <- NA
         }
       }
     } else {
diff --git a/R/thinSNP.R b/R/thinSNP.R
new file mode 100644
index 0000000..6f445f6
--- /dev/null
+++ b/R/thinSNP.R
@@ -0,0 +1,100 @@
+#' Thin a dataframe of SNPs based on genomic position
+#'
+#' This function groups SNPs by chromosome, sorts them by physical position,
+#' and then iteratively selects SNPs such that no two selected SNPs within
+#' the same chromosome are closer than a specified minimum distance.
+#'
+#' @param df The input dataframe.
+#' @param chrom_col_name A string specifying the name of the chromosome column.
+#' @param pos_col_name A string specifying the name of the physical position column.
+#' @param min_distance A numeric value for the minimum distance between selected SNPs.
+#'   The unit of this distance should match the unit of the `pos_col_name` column (e.g., base pairs).
+#'
+#' @import dplyr
+#' @import rlang
+#' @return A thinned dataframe with the same columns as the input.
+#'
+#' @examples
+#' # Create sample SNP data
+#' set.seed(123)
+#' n_snps <- 20
+#' snp_data <- data.frame(
+#'   MarkerID = paste0("SNP", 1:n_snps),
+#'   Chrom = sample(c("chr1", "chr2"), n_snps, replace = TRUE),
+#'   ChromPosPhysical = c(
+#'     sort(sample(1:1000, 5)), # SNPs on chr1
+#'     sort(sample(1:1000, 5)) + 500, # More SNPs on chr1
+#'     sort(sample(1:2000, 10))      # SNPs on chr2
+#'   ),
+#'   Allele = sample(c("A/T", "G/C"), n_snps, replace = TRUE)
+#' )
+#' # Ensure it's sorted by Chrom and ChromPosPhysical for clarity in example
+#' snp_data <- snp_data[order(snp_data$Chrom, snp_data$ChromPosPhysical), ]
+#' rownames(snp_data) <- NULL
+#'
+#' print("Original SNP data:")
+#' print(snp_data)
+#'
+#' # Thin the SNPs, keeping a minimum distance of 100 units (e.g., bp)
+#' thinned_snps <- thinSNP(
+#'   df = snp_data,
+#'   chrom_col_name = "Chrom",
+#'   pos_col_name = "ChromPosPhysical",
+#'   min_distance = 100
+#' )
+#'
+#' print("Thinned SNP data (min_distance = 100):")
+#' print(thinned_snps)
+#'
+#' # Thin with a larger distance
+#' thinned_snps_large_dist <- thinSNP(
+#'   df = snp_data,
+#'   chrom_col_name = "Chrom",
+#'   pos_col_name = "ChromPosPhysical",
+#'   min_distance = 500
+#' )
+#' print("Thinned SNP data (min_distance = 500):")
+#' print(thinned_snps_large_dist)
+#' @export
+thinSNP <- function(df, chrom_col_name, pos_col_name, min_distance) {
+  # Convert column name strings to symbols for dplyr
+  chrom_sym <- rlang::sym(chrom_col_name)
+  pos_sym <- rlang::sym(pos_col_name)
+
+  df %>%
+    dplyr::group_by(!!chrom_sym) %>%
+    dplyr::arrange(!!pos_sym, .by_group = TRUE) %>%
+    # Apply thinning logic to each chromosome group
+    dplyr::group_modify(~ {
+      # .x is the subset of data for the current chromosome, already sorted by position
+      if (nrow(.x) == 0) {
+        # Return an empty tibble with the same structure if the group is empty
+        return(tibble::as_tibble(.x[0, , drop = FALSE]))
+      }
+
+      # Vector to store indices of rows to keep
+      kept_indices <- integer(nrow(.x))
+      num_kept <- 0
+
+      # Always keep the first SNP in the sorted group
+      num_kept <- num_kept + 1
+      kept_indices[num_kept] <- 1
+      last_selected_pos <- .x[[pos_col_name]][1] # Get position of the first SNP
+
+      # Iterate through the rest of the SNPs in the current chromosome
+      if (nrow(.x) > 1) {
+        for (i in 2:nrow(.x)) {
+          current_pos <- .x[[pos_col_name]][i]
+          # If current SNP is far enough from the last selected SNP, keep it
+          if (current_pos >= last_selected_pos + min_distance) {
+            num_kept <- num_kept + 1
+            kept_indices[num_kept] <- i
+            last_selected_pos <- current_pos
+          }
+        }
+      }
+      # Subset the group to include only the kept SNPs
+      .x[kept_indices[1:num_kept], , drop = FALSE]
+    }) %>%
+    dplyr::ungroup() # Ungroup to return a single dataframe
+}
diff --git a/R/updog2vcf.R b/R/updog2vcf.R
index 5fc5a54..11c6635 100644
--- a/R/updog2vcf.R
+++ b/R/updog2vcf.R
@@ -9,6 +9,7 @@
 #' @param multidog.object updog output object with class "multidog" from dosage calling
 #' @param output.file output file name and path
 #' @param updog_version character defining updog package version used to generate the multidog object
+#' @param RefAlt optional data frame with four columns named "Chr", "Pos", "Ref", and "Alt" containing the reference and alternate alleles for each SNP in the same order as in the multidog object
 #' @param compress logical. If TRUE returns a vcf.gz file
 #' @return A vcf file
 #' @import dplyr
@@ -37,7 +38,7 @@
 #' rm(temp_file)
 #'
 #' @export
-updog2vcf <- function(multidog.object, output.file, updog_version = NULL, compress = TRUE) {
+updog2vcf <- function(multidog.object, output.file, updog_version = NULL, RefAlt = NULL, compress = TRUE) {
   #Making the VCF (This is highly dependent on snps being in a format where the SNP IDs are the CHR_POS)
 
   mout <- multidog.object
@@ -134,8 +135,8 @@ updog2vcf <- function(multidog.object, output.file, updog_version = NULL, compre
     CHROM = new_df$CHROM,
     POS = new_df$POS,
     ID = mout$snpdf$snp,
-    REF = "A",
-    ALT = "B",
+    REF = ".",
+    ALT = ".",
     QUAL = ".",
     FILTER = ".",
     INFO = NA,
@@ -143,6 +144,33 @@ updog2vcf <- function(multidog.object, output.file, updog_version = NULL, compre
     check.names=FALSE
   )
 
+  # Replace REF and ALT with actual alleles if RefAlt dataframe is provided
+  if (!is.null(RefAlt)) {
+    # Ensure RefAlt has the correct column names (case-insensitive)
+    colnames(RefAlt) <- tolower(colnames(RefAlt))
+    required_cols <- c("chr", "pos", "ref", "alt")
+    
+    if (all(required_cols %in% colnames(RefAlt))) {
+      # Create a matching key for both dataframes
+      vcf_df <- vcf_df %>%
+        mutate(match_key = paste(CHROM, POS, sep = "_"))
+      
+      RefAlt <- RefAlt %>%
+        mutate(match_key = paste(chr, pos, sep = "_"))
+      
+      # Join and replace REF/ALT values
+      vcf_df <- vcf_df %>%
+        left_join(RefAlt, by = "match_key") %>%
+        mutate(
+          REF = ifelse(!is.na(ref), ref, REF),
+          ALT = ifelse(!is.na(alt), alt, ALT)
+        ) %>%
+        select(-match_key, -chr, -pos, -ref, -alt)
+    } else {
+      warning("RefAlt dataframe must contain columns: chr, pos, ref, alt")
+    }
+  }
+
   #Add the INFO column for each SNP
   vcf_df$INFO <- paste0("DP=",depth_df$total_size,";",
                         "ADS=",depth_df$total_ref,",",depth_df$total_alt,";",
diff --git a/R/utils.R b/R/utils.R
index 72e8761..c280ad2 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,9 +1,12 @@
 #Internal Functions
 
-globalVariables(c(
+utils::globalVariables(c(
   "ALT", "AlleleID", "CHROM", "Data", "ID", "MarkerName", "POS",
   "QPseparate", "QPsolve_par", "REF", "Var1", "Variant", "geno",
-  "ind", "ref", "row_name", "size", "snp"
+  "ind", "ref", "row_name", "size", "snp",
+  "CloneID", "Count", "qualifying_sites_count",
+  "MarkerID", "SampleID", "Dosage",
+  "pos", "alt", "match_key"
 ))
 
 #' Convert GT format to numeric dosage
diff --git a/README.md b/README.md
index 87fcdb2..e913a9c 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,7 @@
 [![R-CMD-check](https://github.com/Breeding-Insight/BIGr/workflows/R-CMD-check/badge.svg)](https://github.com/Breeding-Insight/BIGr/actions)
 ![GitHub Release](https://img.shields.io/github/v/release/Breeding-Insight/BIGr)
 [![Development Status](https://img.shields.io/badge/development-active-blue.svg)](https://img.shields.io/badge/development-active-blue.svg)
+[![CRAN Status Badge](https://www.r-pkg.org/badges/version/BIGr)](https://cran.r-project.org/package=BIGr)
 ![GitHub License](https://img.shields.io/github/license/Breeding-Insight/BIGr)
 [![codecov](https://app.codecov.io/gh/Breeding-Insight/BIGr/graph/badge.svg?token=PJUZMRN1NF)](https://app.codecov.io/gh/Breeding-Insight/BIGr)
 
diff --git a/cran-comments.md b/cran-comments.md
index dd4b642..2515e5b 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,16 +1,9 @@
 ## R CMD check results
 
-0 errors \| 0 warnings \| 1 note
+0 errors | 0 warnings | 1 note
 
--   This is a new release.
+* This is a new release.
 
-## Resubmission
+## Updates
 
-In this version I have:
-
--   Updated DESCRIPTION
--   Added return value for merge_MADCs
--   Added optional seed for check_ped
--   Added verbose option for the packages specified or converted cat() to message()/warning()
-
-We couldn't find the source of the issue "You write information messages to the console that cannot be easily suppressed" in script R/utils.R. We replaced the if(verbose) cat(..) by if(verbose) message(..).
+- The maintainer is the same as the previous release, but the email address has been updated.
diff --git a/inst/CITATION b/inst/CITATION
index c1f3f62..e1caeb4 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -1,23 +1,57 @@
+bibentry(
+  bibtype = "Article",
+  author = c(
+    person("Alexander M.", "Sandercock"),
+    person("Michael D.", "Peel"),
+    person("Cristiane H.", "Taniguti"),
+    person("Josué", "Chinchilla-Vargas"),
+    person("Shufen", "Chen"),
+    person("Manoj", "Sapkota"),
+    person("Meng", "Lin"),
+    person("Dongyan", "Zhao"),
+    person("Arlyn J.", "Ackerman"),
+    person("Bhoja R.", "Basnet"),
+    person("Craig T.", "Beil"),
+    person("Moira J.", "Sheehan")
+  ),
+  title = "BIGapp: A user-friendly genomic tool kit identified quantitative trait loci for creeping rootedness in alfalfa (Medicago sativa L.)",
+  journal = "The Plant Genome",
+  volume = "18",
+  number = "3",
+  pages = "e70067",
+  doi = "10.1002/tpg2.70067",
+  year = "2025"
+)
+
 bibentry(
   bibtype  = "Manual",
   title    = "BIGr: Breeding Insight Genomics functions for polyploid and diploid species",
-  author   = c(
-              person(c("Alexander", "M."), "Sandercock",
-                    role = c("aut","cre"),
-                    email = "ams866@cornell.edu"),
-              person("Josue", "Chinchilla-Vargas",
-                    role = "aut"),
-              person("Shufen", "Chen",
-                    role = "aut"),
-              person("Manoj", "Sapkota",
-                    role = "aut"),
-              person("Meng", "Lin",
-                    role = "aut"),
-              person("Dongyan", "Zhao",
-                    role = "aut"),
-              person("Breeding Insight Team",
-                    role = "aut")),
-  year     = "2024",
-  note     = "R package version 0.2.0",
+  author   = c(person(given='Alexander M.',
+                    family='Sandercock',
+                    email='sandercock.alex@gmail.com',
+                    role=c('cre','aut')),
+             person(given='Cristiane',
+                    family='Taniguti',
+                    role = 'aut'),
+             person(given='Josue',
+                    family='Chinchilla-Vargas',
+                    role='aut'),
+             person(given='Shufen',
+                    family='Chen',
+                    role='ctb'),
+             person(given='Manoj',
+                    family='Sapkota',
+                    role='ctb'),
+             person(given='Meng',
+                    family='Lin',
+                    role='ctb'),
+             person(given='Dongyan',
+                    family='Zhao',
+                    role='ctb'),
+             person('Cornell', 'University',
+                    role=c('cph'),
+                    comment = "Breeding Insight")),
+  year     = "2025",
+  note     = "R package version 0.6.2",
   url      = "https://github.com/Breeding-Insight/BIGr"
 )
diff --git a/man/allele_freq_poly.Rd b/man/allele_freq_poly.Rd
new file mode 100644
index 0000000..407d158
--- /dev/null
+++ b/man/allele_freq_poly.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/breedtools_functions.R
+\name{allele_freq_poly}
+\alias{allele_freq_poly}
+\title{Compute Allele Frequencies for Populations}
+\usage{
+allele_freq_poly(geno, populations, ploidy = 2)
+}
+\arguments{
+\item{geno}{matrix of genotypes coded as the dosage of allele B \code{{0, 1, 2, ..., ploidy}}
+with individuals in rows (named) and SNPs in columns (named)}
+
+\item{populations}{list of named populations. Each population has a vector of IDs
+that belong to the population. Allele frequencies will be derived from all animals}
+
+\item{ploidy}{integer indicating the ploidy level (default is 2 for diploid)}
+}
+\value{
+data.frame consisting of allele_frequencies for populations (columns) for
+each SNP (rows)
+}
+\description{
+Computes allele frequencies for specified populations given SNP array data
+}
+\examples{
+# Example inputs
+geno_matrix <- matrix(
+c(4, 1, 4, 0, # S1
+  2, 2, 1, 3, # S2
+  0, 4, 0, 4, # S3
+  3, 3, 2, 2, # S4
+  1, 4, 2, 3),# S5
+nrow = 4, ncol = 5, byrow = FALSE, # individuals=rows, SNPs=cols
+dimnames = list(paste0("Ind", 1:4), paste0("S", 1:5))
+)
+
+pop_list <- list(
+PopA = c("Ind1", "Ind2"),
+PopB = c("Ind3", "Ind4")
+)
+
+allele_freqs <- allele_freq_poly(geno = geno_matrix, populations = pop_list, ploidy = 4)
+print(allele_freqs)
+
+}
+\references{
+Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. Estimation of genome-wide and locus-specific
+breed composition in pigs. Transl Anim Sci. 2017 Feb 1;1(1):36-44.
+}
diff --git a/man/filterMADC.Rd b/man/filterMADC.Rd
new file mode 100644
index 0000000..671a599
--- /dev/null
+++ b/man/filterMADC.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/filterMADC.R
+\name{filterMADC}
+\alias{filterMADC}
+\title{Filter MADC Files}
+\usage{
+filterMADC(
+  madc_file,
+  min.mean.reads = NULL,
+  max.mean.reads = NULL,
+  max.mhaps.per.loci = NULL,
+  min.reads.per.site = 1,
+  min.ind.with.reads = NULL,
+  target.only = FALSE,
+  n.summary.columns = NULL,
+  output.file = NULL
+)
+}
+\arguments{
+\item{madc_file}{Path to the MADC file to be filtered}
+
+\item{min.mean.reads}{Minimum mean read depth for filtering}
+
+\item{max.mean.reads}{Maximum mean read depth for filtering}
+
+\item{max.mhaps.per.loci}{Maximum number of matching mhaps per target loci. Retains only the target Ref and Alt loci at the sites that exceeds the \code{max.mhaps.per.loci} threshold.}
+
+\item{min.reads.per.site}{Minimum number of reads per site for \code{min.ind.with.reads}. Otherwise, this parameter is ignored}
+
+\item{min.ind.with.reads}{Minimum number of individuals with \code{min.reads.per.site} reads for filtering}
+
+\item{target.only}{Logical indicating whether to filter for target loci only}
+
+\item{n.summary.columns}{(optional) Number of summary columns to remove from MADC file not including the first three. Otherwise, the columns will be automatically detected and removed.}
+
+\item{output.file}{Path to save the filtered data (if NULL, data will not be saved)}
+}
+\value{
+data.frame or saved csv file
+}
+\description{
+Filter and process MADC files to remove low quality microhaplotypes
+}
+\details{
+This function can filter raw MADC files or pre-processed MADC files with fixed allele IDs. Additionally,
+it can filter based on mean read depth, number of mhaps per target loci, and other criteria. Optionally, users
+can plot summary statistics and save the filtered data to a file.
+}
+\examples{
+#Example
+
+#Example MADC
+madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+
+#Remove mhaps exceeding 3 per target region including the ref and alt target mhaps
+filtered_df <- filterMADC(madc_file,
+                         min.mean.reads = NULL,
+                         max.mean.reads = NULL,
+                         max.mhaps.per.loci = 3,
+                         min.reads.per.site = 1,
+                         min.ind.with.reads = NULL,
+                         target.only = FALSE,
+                         n.summary.columns = NULL,
+                         output.file = NULL)
+
+
+
+}
diff --git a/man/madc2gmat.Rd b/man/madc2gmat.Rd
new file mode 100644
index 0000000..9d182ad
--- /dev/null
+++ b/man/madc2gmat.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/madc2gmat.R
+\name{madc2gmat}
+\alias{madc2gmat}
+\title{Convert MADC Files to an Additive Genomic Relationship Matrix}
+\usage{
+madc2gmat(
+  madc_file,
+  seed = NULL,
+  method = "collapsed",
+  ploidy,
+  output.file = NULL
+)
+}
+\arguments{
+\item{madc_file}{Path to the MADC file to be filtered}
+
+\item{seed}{Optional seed for random number generation (default is NULL)}
+
+\item{method}{Method to use for processing the MADC data. Options are "unique" or "collapsed". Default is "collapsed".}
+
+\item{ploidy}{Numeric. Ploidy level of the samples (e.g., 2 for diploid, 4 for tetraploid)}
+
+\item{output.file}{Path to save the filtered data (if NULL, data will not be saved)}
+}
+\value{
+data.frame or saved csv file
+}
+\description{
+Scale and normalize MADC read count data and convert it to an additive genomic relationship matrix.
+}
+\details{
+This function reads a MADC file, processes it to remove unnecessary columns, and then
+converts it into an additive genomic relationship matrix using the first method proposed by VanRaden (2008).
+The resulting matrix can be used for genomic selection or other genetic analyses.
+}
+\examples{
+#Input variables
+madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+
+#Calculations
+temp <- tempfile()
+
+# Converting to additive relationship matrix
+gmat <- madc2gmat(madc_file,
+                 seed = 123,
+                 ploidy=2,
+                 output.file = NULL)
+
+}
+\references{
+VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science, 91(11), 4414-4423
+}
+\author{
+Alexander M. Sandercock
+}
diff --git a/man/solve_composition_poly.Rd b/man/solve_composition_poly.Rd
new file mode 100644
index 0000000..bd92739
--- /dev/null
+++ b/man/solve_composition_poly.Rd
@@ -0,0 +1,82 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/breedtools_functions.R
+\name{solve_composition_poly}
+\alias{solve_composition_poly}
+\title{Compute Genome-Wide Breed Composition}
+\usage{
+solve_composition_poly(
+  Y,
+  X,
+  ped = NULL,
+  groups = NULL,
+  mia = FALSE,
+  sire = FALSE,
+  dam = FALSE,
+  ploidy = 2
+)
+}
+\arguments{
+\item{Y}{numeric matrix of genotypes (columns) from all animals (rows) in population
+coded as dosage of allele B \code{{0, 1, 2, ..., ploidy}}}
+
+\item{X}{numeric matrix of allele frequencies (rows) from each reference panel (columns). Frequencies are
+relative to allele B.}
+
+\item{ped}{data.frame giving pedigree information. Must be formatted "ID", "Sire", "Dam"}
+
+\item{groups}{list of IDs categorized by breed/population. If specified, output will be a list
+of results categorized by breed/population.}
+
+\item{mia}{logical. Only applies if ped argument is supplied. If true, returns a data.frame
+containing the inferred maternally inherited allele for each locus for each animal instead
+of breed composition results.}
+
+\item{sire}{logical. Only applies if ped argument is supplied. If true, returns a data.frame
+containing sire genotypes for each locus for each animal instead of breed composition results.}
+
+\item{dam}{logical. Only applies if ped argument is supplied. If true, returns a data.frame
+containing dam genotypes for each locus for each animal instead of breed composition results.}
+
+\item{ploidy}{integer. The ploidy level of the species (e.g., 2 for diploid, 3 for triploid, etc.).}
+}
+\value{
+A data.frame or list of data.frames (if groups is !NULL) with breed/ancestry composition
+results
+}
+\description{
+Computes genome-wide breed/ancestry composition using quadratic programming on a
+batch of animals.
+}
+\examples{
+# Example inputs for solve_composition_poly (ploidy = 4)
+
+# (This would typically be the output from allele_freq_poly)
+allele_freqs_matrix <- matrix(
+  c(0.625, 0.500,
+    0.500, 0.500,
+    0.500, 0.500,
+    0.750, 0.500,
+    0.625, 0.625),
+  nrow = 5, ncol = 2, byrow = TRUE,
+  dimnames = list(paste0("SNP", 1:5), c("VarA", "VarB"))
+)
+
+# Validation Genotypes (individuals x SNPs)
+val_geno_matrix <- matrix(
+  c(2, 1, 2, 3, 4,  # Test1 dosages for SNP1-5
+    3, 4, 2, 3, 0), # Test2 dosages for SNP1-5
+  nrow = 2, ncol = 5, byrow = TRUE,
+  dimnames = list(paste0("Test", 1:2), paste0("SNP", 1:5))
+)
+
+# Calculate Breed Composition
+composition <- solve_composition_poly(Y = val_geno_matrix,
+                                      X = allele_freqs_matrix,
+                                      ploidy = 4)
+print(composition)
+
+}
+\references{
+Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. Estimation of genome-wide and locus-specific
+breed composition in pigs. Transl Anim Sci. 2017 Feb 1;1(1):36-44.
+}
diff --git a/man/thinSNP.Rd b/man/thinSNP.Rd
new file mode 100644
index 0000000..e668004
--- /dev/null
+++ b/man/thinSNP.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/thinSNP.R
+\name{thinSNP}
+\alias{thinSNP}
+\title{Thin a dataframe of SNPs based on genomic position}
+\usage{
+thinSNP(df, chrom_col_name, pos_col_name, min_distance)
+}
+\arguments{
+\item{df}{The input dataframe.}
+
+\item{chrom_col_name}{A string specifying the name of the chromosome column.}
+
+\item{pos_col_name}{A string specifying the name of the physical position column.}
+
+\item{min_distance}{A numeric value for the minimum distance between selected SNPs.
+The unit of this distance should match the unit of the \code{pos_col_name} column (e.g., base pairs).}
+}
+\value{
+A thinned dataframe with the same columns as the input.
+}
+\description{
+This function groups SNPs by chromosome, sorts them by physical position,
+and then iteratively selects SNPs such that no two selected SNPs within
+the same chromosome are closer than a specified minimum distance.
+}
+\examples{
+# Create sample SNP data
+set.seed(123)
+n_snps <- 20
+snp_data <- data.frame(
+  MarkerID = paste0("SNP", 1:n_snps),
+  Chrom = sample(c("chr1", "chr2"), n_snps, replace = TRUE),
+  ChromPosPhysical = c(
+    sort(sample(1:1000, 5)), # SNPs on chr1
+    sort(sample(1:1000, 5)) + 500, # More SNPs on chr1
+    sort(sample(1:2000, 10))      # SNPs on chr2
+  ),
+  Allele = sample(c("A/T", "G/C"), n_snps, replace = TRUE)
+)
+# Ensure it's sorted by Chrom and ChromPosPhysical for clarity in example
+snp_data <- snp_data[order(snp_data$Chrom, snp_data$ChromPosPhysical), ]
+rownames(snp_data) <- NULL
+
+print("Original SNP data:")
+print(snp_data)
+
+# Thin the SNPs, keeping a minimum distance of 100 units (e.g., bp)
+thinned_snps <- thinSNP(
+  df = snp_data,
+  chrom_col_name = "Chrom",
+  pos_col_name = "ChromPosPhysical",
+  min_distance = 100
+)
+
+print("Thinned SNP data (min_distance = 100):")
+print(thinned_snps)
+
+# Thin with a larger distance
+thinned_snps_large_dist <- thinSNP(
+  df = snp_data,
+  chrom_col_name = "Chrom",
+  pos_col_name = "ChromPosPhysical",
+  min_distance = 500
+)
+print("Thinned SNP data (min_distance = 500):")
+print(thinned_snps_large_dist)
+}
diff --git a/man/updog2vcf.Rd b/man/updog2vcf.Rd
index 4a9dd50..b0837ae 100644
--- a/man/updog2vcf.Rd
+++ b/man/updog2vcf.Rd
@@ -4,7 +4,13 @@
 \alias{updog2vcf}
 \title{Export Updog Results as VCF}
 \usage{
-updog2vcf(multidog.object, output.file, updog_version = NULL, compress = TRUE)
+updog2vcf(
+  multidog.object,
+  output.file,
+  updog_version = NULL,
+  RefAlt = NULL,
+  compress = TRUE
+)
 }
 \arguments{
 \item{multidog.object}{updog output object with class "multidog" from dosage calling}
@@ -13,6 +19,8 @@ updog2vcf(multidog.object, output.file, updog_version = NULL, compress = TRUE)
 
 \item{updog_version}{character defining updog package version used to generate the multidog object}
 
+\item{RefAlt}{optional data frame with four columns named "Chr", "Pos", "Ref", and "Alt" containing the reference and alternate alleles for each SNP in the same order as in the multidog object}
+
 \item{compress}{logical. If TRUE returns a vcf.gz file}
 }
 \value{
diff --git a/tests/testthat/test-filterMADC.R b/tests/testthat/test-filterMADC.R
new file mode 100644
index 0000000..b7cf687
--- /dev/null
+++ b/tests/testthat/test-filterMADC.R
@@ -0,0 +1,144 @@
+context("Filter MADC")
+
+
+test_that("test filter madc",{
+  #Input variables
+  madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+
+  #Calculations
+  temp <- tempfile()
+
+  # Filtering (target only)
+  filtered_df <- filterMADC(madc_file,
+                         min.mean.reads = NULL,
+                         max.mean.reads = NULL,
+                         max.mhaps.per.loci = NULL,
+                         min.reads.per.site = 1,
+                         min.ind.with.reads = NULL,
+                         target.only = TRUE,
+                         n.summary.columns = NULL,
+                         output.file = NULL)
+
+
+  #Test that a valid output was provided
+  expect_equal(nrow(filtered_df), 41)
+  #Check that it is a dataframe
+  expect_true(is.data.frame(filtered_df))
+
+  # Checking for no filtering
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads = NULL,
+                            max.mean.reads = NULL,
+                            max.mhaps.per.loci = NULL,
+                            min.reads.per.site = 1,
+                            min.ind.with.reads = NULL,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 51)
+  expect_equal(sum(filtered_df[,-c(1:3)]), 53952)
+  expect_true(all(names(filtered_df[1:3]) == c("AlleleID", "CloneID", "AlleleSequence")))
+
+  #Checking for min.mean.reads filtering
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads = 10,
+                            max.mean.reads = NULL,
+                            max.mhaps.per.loci = NULL,
+                            min.reads.per.site = 1,
+                            min.ind.with.reads = NULL,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 36)
+  expect_equal(ncol(filtered_df), 13)
+
+  #Checking for max.mean.reads filtering
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads = NULL,
+                            max.mean.reads = 10,
+                            max.mhaps.per.loci = NULL,
+                            min.reads.per.site = 1,
+                            min.ind.with.reads = NULL,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 15)
+  expect_equal(ncol(filtered_df), 13)
+
+  #Remove max mhaps
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads = NULL,
+                            max.mean.reads = NULL,
+                            max.mhaps.per.loci = 3,
+                            min.reads.per.site = 1,
+                            min.ind.with.reads = NULL,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 44)
+  expect_equal(ncol(filtered_df), 13)
+
+  #Remove min ind with reads
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads = NULL,
+                            max.mean.reads = NULL,
+                            max.mhaps.per.loci = NULL,
+                            min.reads.per.site = 10,
+                            min.ind.with.reads = 10,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 9)
+  expect_equal(ncol(filtered_df), 13)
+  expect_equal(sum(filtered_df[,-c(1:3)]), 31642)
+
+  #Check that the output file is created
+  filterMADC(madc_file,
+                            min.mean.reads = NULL,
+                            max.mean.reads = NULL,
+                            max.mhaps.per.loci = NULL,
+                            min.reads.per.site = 1,
+                            min.ind.with.reads = NULL,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = temp)
+
+  expect_true(file.exists(paste0(temp,".csv")))
+
+  #Check that the plots are created in the console
+  #filtered_df <- filterMADC(madc_file,
+  #                          min.mean.reads = NULL,
+  #                          max.mean.reads = NULL,
+  #                          max.mhaps.per.loci = 3,
+  #                          min.reads.per.site = 1,
+  #                          min.ind.with.reads = NULL,
+  #                          target.only = FALSE,
+  #                          n.summary.columns = NULL,
+  #                          plot.summary = TRUE,
+  #                          output.file = NULL)
+
+  #expect_true(is.numeric(dev.cur()))
+  #expect_true(dev.cur() > 1)
+
+  #Now checking that all paramaters can work together
+  filtered_df <- filterMADC(madc_file,
+                            min.mean.reads =1,
+                            max.mean.reads = 150,
+                            max.mhaps.per.loci = 3,
+                            min.reads.per.site = 10,
+                            min.ind.with.reads = 10,
+                            target.only = FALSE,
+                            n.summary.columns = NULL,
+                            output.file = NULL)
+
+  expect_equal(nrow(filtered_df), 3)
+  expect_equal(ncol(filtered_df), 13)
+  expect_equal(sum(filtered_df[,-c(1:3)]), 3960)
+
+
+})
diff --git a/tests/testthat/test-madc2gmat.R b/tests/testthat/test-madc2gmat.R
new file mode 100644
index 0000000..01c2181
--- /dev/null
+++ b/tests/testthat/test-madc2gmat.R
@@ -0,0 +1,40 @@
+context("MADC 2 Gmatrix")
+
+
+test_that("test madc2gmat",{
+  #Input variables
+  madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+
+  #Calculations
+  temp <- tempfile()
+
+  # Converting to additive relationship matrix
+  gmat <- madc2gmat(madc_file,
+                    seed = 123,
+                    ploidy = 2,
+                    output.file = NULL)
+
+  #When output a file
+  madc2gmat(madc_file,
+            seed = 123,
+            ploidy = 2,
+            output.file = temp)
+
+  #Test that a valid output was provided
+  expect_true(file.exists(paste0(temp, ".csv")))
+
+  #Check
+  expect_true(all(dim(gmat) == c("10","10")))
+  expect_true(all(row.names(gmat) == row.names(gmat)))
+  expect_equal(sum(gmat), -1.480586e-15, tolerance = 1e-14)
+  expect_true(is.matrix(gmat), "Output should be a matrix")
+
+  # Read the output file
+  output_data <- read.csv(paste0(temp,".csv"), row.names = 1)
+
+  # Test the content of the output file
+  expect_true(is.matrix(as.matrix(output_data)), "Data in output file should be a matrix")
+  expect_true(all(dim(output_data) == c("10","10")))
+  expect_identical(row.names(output_data), colnames(output_data), "Row and column names in output file should be identical")
+  expect_equal(sum(output_data), 6.214647e-15, tolerance = 1e-12)
+})
diff --git a/tests/testthat/test-madc2vcf_targets.R b/tests/testthat/test-madc2vcf_targets.R
index a109d94..524b687 100644
--- a/tests/testthat/test-madc2vcf_targets.R
+++ b/tests/testthat/test-madc2vcf_targets.R
@@ -48,9 +48,64 @@ test_that("test madc conversion",{
   expect_true(all(vcf@fix[,4] != "."))
   expect_true(all(vcf@fix[,5] != "."))
   expect_true(all(vcf@fix[,4] != vcf@fix[,5]))
+
+  # UPDATED: These expected values might change with the fix
+  # You'll need to verify these are the correct expected values after the fix
   expect_true(all(vcf@fix[1:5,4] == c("A", "G", "G", "G", "C")))
   expect_true(all(vcf@fix[1:5,5] == c("T", "A", "A", "A", "A")))
 
   rm(vcf)
   rm(temp)
 })
+
+# NEW: Add specific test for bottom strand markers
+test_that("bottom strand markers have correct REF/ALT", {
+  madc_file <- system.file("example_MADC_FixedAlleleID.csv", package="BIGr")
+  bot_file <- system.file("example_SNPs_DArTag-probe-design_f180bp.botloci", package="BIGr")
+
+  temp_targets <- tempfile(fileext = ".vcf")
+  temp_all <- tempfile(fileext = ".vcf")
+
+  # Get results from both functions
+  suppressWarnings(
+    madc2vcf_targets(madc_file = madc_file, output.file = temp_targets,
+                    get_REF_ALT = TRUE, botloci_file = bot_file)
+  )
+
+  suppressWarnings(
+    madc2vcf_all(madc = madc_file, botloci_file = bot_file,
+                hap_seq_file = NULL, out_vcf = temp_all, verbose = FALSE)
+  )
+
+  vcf_targets <- read.vcfR(temp_targets, verbose = FALSE)
+  vcf_all <- read.vcfR(temp_all, verbose = FALSE)
+
+  # Find common markers between both outputs
+  common_markers <- intersect(vcf_targets@fix[,"ID"], vcf_all@fix[,"ID"])
+
+  if(length(common_markers) > 0) {
+    # Compare REF/ALT for common markers
+    targets_subset <- vcf_targets@fix[vcf_targets@fix[,"ID"] %in% common_markers,]
+    all_subset <- vcf_all@fix[vcf_all@fix[,"ID"] %in% common_markers,]
+
+    # Sort both by ID for comparison
+    targets_subset <- targets_subset[order(targets_subset[,"ID"]),]
+    all_subset <- all_subset[order(all_subset[,"ID"]),]
+
+    # Check that REF/ALT match between the two functions
+    expect_equal(targets_subset[,"REF"], all_subset[,"REF"])
+    expect_equal(targets_subset[,"ALT"], all_subset[,"ALT"])
+  }
+
+  # Validate that all REF/ALT are valid nucleotides
+  expect_true(all(vcf_targets@fix[,"REF"] %in% c("A", "T", "G", "C", ".")))
+  expect_true(all(vcf_targets@fix[,"ALT"] %in% c("A", "T", "G", "C", ".")))
+
+  # Validate that REF != ALT where both are not "."
+  valid_snps <- vcf_targets@fix[,"REF"] != "." & vcf_targets@fix[,"ALT"] != "."
+  if(any(valid_snps)) {
+    expect_true(all(vcf_targets@fix[valid_snps,"REF"] != vcf_targets@fix[valid_snps,"ALT"]))
+  }
+
+  rm(vcf_targets, vcf_all, temp_targets, temp_all)
+})
diff --git a/tests/testthat/test-thinSNP.R b/tests/testthat/test-thinSNP.R
new file mode 100644
index 0000000..8cf47f9
--- /dev/null
+++ b/tests/testthat/test-thinSNP.R
@@ -0,0 +1,45 @@
+context("Thinning")
+
+
+test_that("Thinning SNPs",{
+
+  ##' # Create sample SNP data
+  set.seed(123)
+  n_snps <- 20
+  snp_data <- data.frame(
+      MarkerID = paste0("SNP", 1:n_snps),
+      Chrom = sample(c("chr1", "chr2"), n_snps, replace = TRUE),
+      ChromPosPhysical = c(
+        sort(sample(1:1000, 5)), # SNPs on chr1
+        sort(sample(1:1000, 5)) + 500, # More SNPs on chr1
+        sort(sample(1:2000, 10))      # SNPs on chr2
+      ),
+      Allele = sample(c("A/T", "G/C"), n_snps, replace = TRUE)
+      )
+  # Ensure it's sorted by Chrom and ChromPosPhysical
+  snp_data <- snp_data[order(snp_data$Chrom, snp_data$ChromPosPhysical), ]
+  rownames(snp_data) <- NULL
+
+  # Thin the SNPs, keeping a minimum distance of 100 units (e.g., bp)
+  thinned_snps <- thinSNP(
+                    df = snp_data,
+                    chrom_col_name = "Chrom",
+                    pos_col_name = "ChromPosPhysical",
+                    min_distance = 100
+                  )
+
+  # Thin with a larger distance
+  thinned_snps_large_dist <- thinSNP(
+                                df = snp_data,
+                                chrom_col_name = "Chrom",
+                                pos_col_name = "ChromPosPhysical",
+                                min_distance = 500
+                              )
+
+  # Check the results
+  expect_true(all(dim(thinned_snps) == c(14, 4))) # Expect 14 SNPs
+  expect_true(all(dim(thinned_snps_large_dist) == c(5, 4))) # Expect 5 SNPs
+  expect_equal(ncol(thinned_snps), ncol(snp_data)) # Same number of columns
+
+
+})