coralesce Tutorial

Author

Colin Shea ([email protected])

What does coralesce do?

coralesce is an R package that uses SNP data to assign coral colonies to genets, calculate population-level mean kinship and gene diversity, calculate individual-level mean kinship, and generate a customized list of the n least-related colonies of a given species based on individual-level mean kinship.

There are two ways to work with the package:

  • Folder-based wrappers (runGenets, runKinship): point them at a folder of CSV files and they read, analyze, and write results to disk. Best for routine batch processing.
  • In-memory functions (computeGenets, computeKinship, collapseToGenets): take a data frame directly and return results, with no folders involved. Best when your genotypes are already in R, or when you want a clone-corrected (one-colony-per-genet) analysis (see Section 9).

Basic steps to use the folder-based workflow for a single species:

  1. Install coralesce if you haven’t done so already
  2. Create a folder for the target species
  3. In the target species’ folder, create a Data and Results folder
  4. Put the target species’ SNP data (a .csv) in the Data folder
  5. Open R and set the working directory to the target species’ folder
  6. Load coralesce and run runGenets and/or runKinship
  7. Check the target species’ Results folder for the results

Getting started

For the wrapper functions to work properly, you’ll need to:

  1. Choose a location for storing the genetic data for each species.
  2. Create folders in that location for each species.
  3. In each species’ folder, there must be two additional folders, one called Data and another called Results. The Results folder starts empty; the Data folder holds a .csv of that species’ SNP data. The folder-name matching is case-insensitive, so Data/DATA/data all work.
  4. You don’t need to import the data files yourself: runGenets and runKinship read them automatically provided the working directory is set to the focal species’ folder.
  5. If there is more than one .csv in Data, both wrappers loop through the files, analyze each separately, and write uniquely named results to Results.

If you don’t have it installed already, you can install coralesce from its R-universe repository:

install.packages("coralesce", repos = "https://colinpshea.r-universe.dev")

Once installed, load and attach it like any other package. Note that dplyr and tidyr are among the imported dependencies.

library(coralesce)

There are many functions available in coralesce, and there is also a table called IUPAC used to convert base-pair combinations such as G:G to single-letter codes such as G.

ls("package:coralesce")
 [1] "calcPercentNotNull"              "checkforAllowableData"          
 [3] "classifyAllelePairs"             "classifyAllelePairsOthers"      
 [5] "collapseToGenets"                "computeGenets"                  
 [7] "computeKinship"                  "convertBasePairstoCodes"        
 [9] "determineAllAlleleMatches"       "determineAllAlleleMatchesOthers"
[11] "find_dups"                       "groupByGenets"                  
[13] "handleError_allZeros"            "handleError_ColumnContract"     
[15] "handleError_ProhibitedData"      "isolateAllNAColonies"           
[17] "IUPAC"                           "kinshipCalcs"                   
[19] "kinshipCalcsNoInvar"             "omitInvariantLoci"              
[21] "readGeneticData"                 "returnGenetIdentity"            
[23] "runGenets"                       "runKinship"

Below is the IUPAC table, which becomes available when the package is loaded.

 

The input data format

Both wrappers read .csv files from the Data folder automatically, but let’s look at an example to see how the data must be formatted.

The input contract (required): the first column must be named Coral_ID and the second column must be named MatchMaker_Index and contain whole numbers (an integer bookkeeping key that maps back to your database). Every remaining column is a locus, holding paired-allele SNP calls such as A:G, or ? / blank for missing data. If the first two columns are not named exactly Coral_ID and MatchMaker_Index, or MatchMaker_Index is not an integer, the functions stop with an informative error rather than guessing (see Section 11). Any extra, non-SNP columns (e.g., site names) are ignored and dropped automatically.

Here, rows are individual colonies and columns are Coral_ID, MatchMaker_Index, and the SNP data at n loci.

 

Assigning colonies to genets

runGenets assigns each coral colony to a genet, either with other colonies or, at minimum, as the sole member of its own genet. Two colonies are grouped into the same genet if (a) the percentage of loci whose alleles match is ≥ a user-defined threshold (PctMatchThreshold) and (b) the percentage of loci with valid SNP data across the compared colonies is ≥ a second user-defined threshold (PctNotNullThreshold).

To use runGenets properly:

  1. Make sure the working directory is the species folder, its data are in the Data folder, and a Results folder exists.
  2. Supply PctMatchThreshold (the percent allele match at which two colonies are considered clones) and PctNotNullThreshold (the minimum percent of loci with valid data for a comparison to be trusted). Below we use PctMatchThreshold = 99 and PctNotNullThreshold = 78; good values may differ by species.
  3. Optionally set getPairwiseAlleleMatches = TRUE to also write the full pairwise comparison table (a large file) to Results. The default is FALSE.
genets <- runGenets(PctMatchThreshold = 99, PctNotNullThreshold = 78, getPairwiseAlleleMatches = TRUE)
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID

runGenets returns a named list; genetAssignments is a list with one entry per input file. With a single data file, the assignments are the first (and only) element.

Take a look at the resulting genetAssignments table. There are five columns: Coral_ID; MatchMaker_Index; genet (a 4-letter species code followed by a 5-digit number); pctNull (the percentage of loci with NA SNPs for that colony); and AdequateData (Yes = the colony had adequate data in a pairwise comparison with at least one other colony; No = the colony had inadequate data in all comparisons and could not be assigned to a genet).

 

Now the corresponding pairwiseAlleleMatches table (this list is absent if getPairwiseAlleleMatches = FALSE). It has five columns: coral1 and coral2 (the two colonies in the comparison), pctMatch (the percentage of loci whose SNPs matched), pctNotNull (the percentage of loci with valid data for that comparison), and PartOfGenet (Yes/No: are the two colonies in the same genet?).

 

Calculating individual- and population-level mean kinship and gene diversity

To use runKinship properly, make sure the working directory, Data, and Results folders are set up as above, then specify the arguments. By default subset = FALSE and targetN = NULL, so a bare runKinship() returns population-level mean kinship and gene diversity plus individual-level mean kinship for every colony.

As a reminder: population-level mean kinship is the average pairwise kinship across all colonies; individual-level mean kinship is the average pairwise kinship for each colony; and population-level mean gene diversity is 1 - population-level mean kinship.

The three calls below are equivalent under the defaults:

kinships_nosubset <- runKinship(subset = FALSE, targetN = NULL)
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID
kinships_nosubset <- runKinship(subset = FALSE)
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID
kinships_nosubset <- runKinship()
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID

runKinship returns a named list; each element (PopAvgMKGD, kinship_init, and — when subset = TRUE — kinship_targetN) is itself a per-file list. Population-level mean kinship and gene diversity live in PopAvgMKGD, and are written to Results with popAvgMKGD in the file name.

 

Individual-level mean kinship is in kinship_init (one row per colony), written to Results with kinship_Init in the file name.

 

Under the defaults there is no subset result, so kinship_targetN is absent:

is.null(kinships_nosubset$kinship_targetN)
[1] TRUE

Kinship for a subset of the population

Setting subset = TRUE and a targetN also returns individual-level mean kinship for the targetN least closely-related colonies. Colonies are removed one at a time (the highest-mean-kinship colony each round) until targetN remain.

Because kinship must be recomputed every time a colony is removed, the smaller targetN is relative to the starting number of colonies, the longer this runs. runKinship with many colonies and a small targetN can take a long time.

Here we target 50% of the colonies (targetN must be ≥ 2):

(N <- round(0.50 * nrow(check)))
[1] 10
kinships_subset <- runKinship(subset = TRUE, targetN = N)
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID

Population-level values (PopAvgMKGD) and full individual-level kinship (kinship_init) are returned as before. The new piece is kinship_targetN: individual-level mean kinship for the retained targetN colonies, written to Results with kinship_targetN in the file name.

 

Working in memory, and clone-corrected diversity

If your genotypes are already in R, you can skip the folders entirely. computeGenets and computeKinship are the in-memory counterparts of runGenets and runKinship: they take a data frame (same format as the input CSV) and return the same results, writing nothing.

# genet assignments, straight from a data frame
ga_mem <- computeGenets(check, PctMatchThreshold = 99, PctNotNullThreshold = 78)

# eligible-pool kinship (all colonies, including any clonal replicates)
computeKinship(check)$PopAvgMKGD
# A tibble: 1 × 2
  PopAvgMK PopAvgGD
     <dbl>    <dbl>
1    0.688    0.312

By default, kinship and gene diversity are computed over all colonies, including multiple ramets (physical fragments) of the same genet — appropriate when every colony is an eligible candidate for breeding or outplanting. If instead you want a clone-corrected estimate (i.e., one colony per genet), which might be more in line with the standard population-genetics view, use collapseToGenets to reduce the data to a single representative per genet (the colony with the most data), then run kinship on that:

reduced <- collapseToGenets(check, ga_mem)   # one colony per genet
computeKinship(reduced)$PopAvgMKGD            # clone-corrected diversity
# A tibble: 1 × 2
  PopAvgMK PopAvgGD
     <dbl>    <dbl>
1    0.686    0.314

The whole clone-corrected workflow — computeGenets → collapseToGenets → computeKinship — runs entirely in memory, with no Data/Results folders required.

Data checks and warnings

Argument errors in runKinship. The calls below stop with an error because the arguments are contradictory or invalid:

Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID
Error : subset = TRUE requires a targetN value.
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID
Error : subset = FALSE but a targetN value was supplied. Set subset = TRUE or leave targetN = NULL.
Columns that do not match the required IUPAC base-pair format (e.g., a site-name column or an invalid base pair) were removed prior to genet and/or kinship calculations.
Offending columns:
Site.Name
Region
Latitude
Longitude
Origin
Facility
Status
Sample.ID
Error : When subset = TRUE, targetN must be >= 2.

Input-column contract. Unlike earlier versions, the package does not silently rename columns. The first column must be named Coral_ID, the second must be named MatchMaker_Index and hold whole numbers. If not, the functions stop with a message identifying the problem (e.g., “Column 2 must be named ‘MatchMaker_Index’ but is named ‘Longitude’…” or “‘MatchMaker_Index’ must contain only whole numbers…”). Fix the input and re-run.

Extraneous / invalid columns. Any column other than the two keys that does not contain valid paired-allele data (e.g., a site-name column, or an invalid base pair) is reported and dropped before analysis, with a message listing the offending columns. If a locus you expected shows up as offending, check that column in your data.

Duplicated colonies. If the same Coral_ID appears in more than one row, runGenets and runKinship stop with an error listing the duplicated colonies. Every colony must be uniquely named.

Colonies with no SNP data. Colonies that are NA/blank/? at every locus are flagged and added to the genet assignment with genet = XXXX_NA (where XXXX is the 4-letter species code), pctNull = 100, and AdequateData = No, so they can be filtered out easily. A message lists any such colonies.

Tips for one or multiple species

Keeping a small R script in each species’ folder makes batch runs easy: set the working directory to that folder and call runGenets / runKinship. In RStudio you can run each as a background job (Source → Source as Background Job), which runs in a clean session with that folder as the working directory, so you can process several species at once without cross-contaminating results.