This is the documentation of SelHap, a pipeline to help selection of new accessions for pangenome projects. It was developed in the research group Domestication Genomics at the Leibniz Institute of Plant Genetics and Crop Research (IPK) Gatersleben.
Selecting N new genomes that maximize haplotype novelty relative to an existing pangenome.
Given a set of candidate samples and a background dataset, SelHap ranks accessions based on their contribution of novel haplotypes.
SelHap is a pipeline designed to identify new accessions that contribute novel haplotypes to an existing pangenome. As pangenomes mature, selecting additional genomes that add genuinely new sequence information becomes increasingly challenging. SelHap addresses this problem by prioritizing samples based on haplotype novelty relative to an existing pangenome dataset.
SelHap builds on the output of IntroBlocker (Wang et al., 2022), a tool that assigns ancestral haplotype groups (AHGs) to genomic bins to investigate fine-scale ancestry and evolutionary history. Using these haplotype clusters, SelHap identifies genomic segments that are not yet represented in a background dataset and ranks candidate samples according to their contribution of novel haplotypes. This differs from traditional approaches that maximize global genetic diversity rather than conditional novelty.
The pipeline requires a VCF file derived from whole-genome shotgun (WGS) data that includes both the existing pangenome accessions (background) and candidate samples. IntroBlocker is first used to cluster haplotypes across genomic bins. SelHap then processes these clusters to generate a ranked list of samples carrying novel haplotypes relative to the background dataset. In addition, SelHap can generate visualizations showing where these novel haplotypes occur across the genome.
Beyond pangenome expansion, SelHap can also be used to explore haplotype frequencies between groups of samples, for example when comparing domesticated and wild germplasm.
SelHap requires a VCF file derived WGS data that includes both the background samples (e.g. existing pangenome accessions) and the candidate samples to be evaluated. All samples should be aligned to the same reference genome before variant calling to ensure consistent genomic coordinates.
The VCF should contain high-quality SNPs distributed across the genome, as these markers are used by IntroBlocker to define haplotype clusters. Standard filtering to remove low-quality variants and sites with excessive missing data is recommended.
The filtered VCF can then be used as input for IntroBlocker.
Instructions on how to run IntroBlocker are available here. Choose window size and SNP threshold according to the species genome size. For our barley project we used 2 Mb windows, cluster_thr=500 and smooth_thr=400. We recommend changing sample IDs to remove any special character other than "-".
SelHap should be run in R, check library requirements above. Just load SelHap_v.1.0.R script to use it.
Commands can be run with test data available on demo folder and further adapted to each use case, mainly changing input files, chrtotal, and number_of_samples.
The following step loads the input files and generates the processed haplotype tables used by SelHap. These tables are derived from the IntroBlocker output and serve as the main input for the three SelHap modes (bin-haps, freq-haps, and select-by-region).
# load SelHap
source("SelHap_v1.0.R")
# input files
chr_list <- fread("chr_list.txt", header=F)
background <- fread("all_samples_pangenome.txt", header=F)[[1]]
tables <- read_tables(path=".")
# generate tables to be used as input for the other functions.
tables <- process_tables(tables, background, chrtotal=7)chr_list is a text file with the list of chromosome names, such as "chr1", "1", or "chr1H", and so on.
background is a text file with sample names as they are in the VCF and IntroBlocker output, one per line. It will define which samples will be used as the background for the analysis. This could be the pangenome accessions, for example.
read_tables() will scan the selected path and read one file per chromosome with matching name like "chrName*.grp", or specifically "chr1H.grp".
Binary-haps mode will select samples based on their presence or absence in the background. This mode is run iteratively, so the chosen sample from iteration 1 will become part of the background until it reaches number of selected iterations (number_of_samples).
selected_samples <- compute_bin_haps(tables = tables,
background = background,
chrtotal = 7, # total number of chromosomes
out_prefix = "example_run_bin", # prefix for output file, optional
number_of_samples = 5, # number of samples to selected, default=1
chr_subset = NULL)chr_subset can be NULL (all chromosomes), a vector of chromosome names from chr_list, or chromosome numbers/indices only.
Files are written using the prefix defined in out_prefix.
In most cases, users will only need the out_prefix_selected_samples_<n>_iterations.txt file, which is also saved to the object selected_samples (table below).
For plotting, the object bin_tables, saved also to out_prefix_binary_tables_first_iteration.Rds, should be used.
| Sample | chr1H | chr2H | chr3H | chr4H | chr5H | chr6H | chr7H | sum_stretches |
|---|---|---|---|---|---|---|---|---|
core150-HOR-15254 |
180 |
50 |
736 |
8250 |
286 |
2941 |
2455 |
14898 |
core150-HOR-16149 |
312 |
202 |
71 |
480 |
198 |
11 |
4918 |
6192 |
core150-HOR-18300 |
3 |
66 |
45 |
8 |
22 |
5855 |
51 |
6050 |
core150-HOR-1048 |
51 |
50 |
894 |
58 |
91 |
3529 |
33 |
4706 |
core150-HOR-14557 |
124 |
2593 |
325 |
5 |
123 |
52 |
227 |
3449 |
Sum of stretches are calculated according to haplotype novelty. Each haplotype is compared to the background and gets 1 for absent in background and 0 for present in background. Missing values are not computed. To account for long novel haplotypes, we apply run-length encoding (RLE) and square the resulting segment lengths to weight longer haplotypes more strongly, so samples with longer novel haplotypes receive a higher score.
Scores for samples that were not selected tend to diminish across iterations, since the background is being fed with novel haplotypes.
- selected_samples
-
Variable defined when calling
compute_bin_haps(). A table containing the samples selected by SelHap and their total number of novel haplotype stretches per chromosome. This object is also written to disk as a text fileout_prefix_selected_samples_<n>_iterations.txt. - bin_tables
-
Binary haplotype tables corresponding to the first iteration. This is the object used for plotting. It is also written to disk as a Rds file
out_prefix_binary_tables_first_iteration.Rds.
- out_prefix_selected_samples_<n>_iterations.txt
-
List of samples selected and per-chromosome stretches after n iterations (same as
selected_samples). - out_prefix_binary_tables_first_iteration.Rds
-
Binary haplotype tables from the first iteration, where all samples are compared to the original background.
- out_prefix_binary_tables_last_iteration.Rds
-
Binary haplotype tables from the final iteration, showing remaining samples. The background includes all previously selected samples.
- out_prefix_all_stretches.txt
-
Table containing the number of novel haplotype stretches per sample across all chromosomes comparing to the initial background.
- out_prefix_stretches_iteration_<n>.txt
-
Table containing the number of novel haplotype stretches per sample across all chromosomes for the remaining samples after the final iteration (i.e. relative to the final background).
When bin-haps is run iteratively, the background changes at each iteration as newly selected samples are added. Because of this, binary haplotype tables generated in later iterations are calculated relative to different backgrounds. To ensure that all samples are plotted relative to the same background, SelHap saves the binary haplotype tables (bin_tables) from the first iteration, where all samples are compared to the original background. These tables should be used for plotting.
# plotting all samples from 1st iteration
plot_results(bin_tables = bin_tables,
plot_selected = F,
out_prefix = "all_samples")
# plotting the 20 selected samples relative to the initial background (1st iteration table)
samples_to_plot <- selected_samples[, .(Sample)]
plot_results(bin_tables = bin_tables,
plot_selected = T,
selected = samples_to_plot,
out_prefix = "20_selected_samples")plot_selected is a logical value to plot all samples (F) from the input table or just a subset (T).
selected must be declared if plot_selected=T. List of sample names to be plotted (default=NULL).
Frequency-haps mode assigns a score from 0 to 10 representing how frequent a haplotype is in the background dataset. Haplotypes absent from the background receive a score of 0, whereas more common haplotypes receive higher values up to 10. These scores are summed and ranked at the end. The higher the final score, the more novel/exclusive haplotypes the sample has.
Unlike bin-haps, this mode is not run iteratively and therefore the ranking and ordering of samples may differ from the iterative selection.
This mode is useful for exploring haplotype frequency patterns between groups of samples, for example when comparing domesticated and wild germplasm or different breeding pools.
freq_haps <- compute_freq_haps(tables = tables,
background = background,
chrtotal = 7,
out_prefix = "example_run_freq",
chr_subset = NULL)Running freq-haps produces two objects in the R environment, which are also written to disk:
- freq_haps
-
table containing all samples and their haplotype frequency scores. Saved as
out_prefix_all_sample_scores_freq_haps.txt. - freq_tables
-
list of per-chromosome frequency tables used to compute the scores. Saved as
out_prefix_freq_tables.Rds.
| Sample | chr1H | chr2H | chr3H | chr4H | chr5H | chr6H | chr7H | sum_scores |
|---|---|---|---|---|---|---|---|---|
core150-HOR-10702 |
972 |
1011 |
768 |
1199 |
775 |
547 |
948 |
6220 |
core50-HOR-14439 |
992 |
998 |
776 |
1180 |
781 |
518 |
963 |
6208 |
core150-HOR-16358 |
995 |
1012 |
789 |
1186 |
592 |
509 |
992 |
6075 |
core50-HOR-12184 |
979 |
1014 |
697 |
1163 |
746 |
540 |
910 |
6049 |
core50-HOR-17077 |
815 |
1000 |
764 |
1180 |
776 |
513 |
971 |
6019 |
# plot all samples
plot_results_freq(freq_tables = freq_tables,
plot_selected = F,
out_prefix = "all")
# plot specific samples
samples_to_plot <- freq[, .(Sample)]
plot_results_freq(freq_tables = freq_tables,
plot_selected = T,
selected = samples_to_plot,
out_prefix = "selected_freq")freq_tables is generated by compute_freq_haps.
plot_selected is a logical value to plot all samples (F) from the input table or just a subset (T).
selected must be declared if plot_selected=T. List of sample names to be plotted (default=NULL).
This mode identifies samples carrying novel haplotypes within a specific genomic region using the bin_tables object generated by compute_bin_haps(). Because it is not iterative, the tables should be generated with only one iteration, ensuring that all samples are compared to the same background.
# for the same region in all chromosomes
regall <- select_by_region(bin_table = bin_tables,
coordmin = 180, # start coordinate of the region
coordmax = 300, # end coordinate of the region
prop = 0.6) # minimum proportion of bins in the selected region that must contain novel haplotypes
# chromosome specific regions
reg1 <- select_by_region(bin_table = bin_tables,
coordmin = 126,
coordmax = 286,
prop = 0.5,
chr_subset = 1) # evaluated chromosome
reg2 <- select_by_region(bin_table = bin_tables,
coordmin = 221,
coordmax = 381,
prop = 0.5,
chr_subset = 2)
merged_regions <- rbind(reg1,reg2)
unique(merged_regions$Sample)coordmin = start coordinate of the region.
coordmax = end coordinate of the region.
prop = minimum proportion of bins in the selected region that must contain novel haplotypes.
chr_subset = chromosome(s) where the region is evaluated. If omitted, the region is evaluated across all chromosomes.