Performing polyploid QTL analysis using polyqtlR

Data structures

As the IBD data are the most central objects in this package it is worth spending a moment describing them. In polyqtlR, IBD probabilities are organised in nested list form. The first level of structure are the linkage groups. In our example dataset, the list has two elements corresponding to the two linkage groups, each of which is itself a list containing the following elements:

• IBDtype The type of IBD probability, either “genotypeIBD” or “haplotypeIBD”

• IBDarray A 3d array of IBD probabilities, with dimensions “locus”,“genotype class”, “individual”

• map The integrated linkage map (marker names and cM positions, in order)

• parental_phase The phasing of the markers in the parental map

• marginal.likelihoods A list of marginal likelihoods of different valencies if method “hmm” was used, otherwise NULL

• valency The predicted valency that maximised the marginal likelihood, per offspring. For method “heur”, NULL

• offspring The offspring names

• biv_dec Whether bivalents only (TRUE) or also multivalents were allowed (FALSE) in the procedure

• gap Here NULL, but later this can hold a numeric value (e.g. 1 cM) if IBDs are ‘splined’ or interpolated.

• genocodes Ordered list of genotype codes used to represent the different genotype classes

• pairing Log likelihoods of each of the different pairing scenarios considered

• ploidy The ploidy of parent 1

• ploidy2 The ploidy of parent 2

• method The method used, either “hmm” (Hidden Markov Model), “heur” (heuristic) or “hmm_TO” (Hidden Markov Model TetraOrigin)

• error The offspring error prior (eps) used in the calculation.

All functions within polyqtlR for estimating or importing IBD probabilities will automatically return them in this form.

HMM for IBD estimation

The estimate_IBD function of polyqtlR with the default setting method = "hmm" estimates offspring IBD probabilities using a Hidden Markov Model (HMM) developed by Zheng et al (Zheng et al. 2016) in the TetraOrigin package but generalised in polyqtlR to multiple ploidy levels. Currently, diploid (2x), triploid (3x = 4x × 2x), tetraploid (4x) and hexaploid (6x) populations are handled. genotypes can either be discrete (i.e. integer marker dosage scores), or probabilistic (i.e. the probabilities of each of the ploidy + 1 possible scores). For triploids, tetraploids and hexaploids, the possibility of incorporating double reduction (Peter M. Bourke et al. 2015) (i.e. allowing for the phenomenon of multivalent pairing) is available, although this can have serious performance implications for hexaploids (and in hexaploids, the option of allowing multivalent pairing in both parents for the same chromosome is by default disabled as it requires very high RAM requirements (> 32 Gb). Use argument full_multivalent_hexa = TRUE to allow multivalents simultaneously at the same position from both hexaploid parents).

Some of the code used to estimate these IBD probabilities has been programmed in C++ to improve performance, and relies on both the Rcpp and RcppArmadillo packages to provide the interface between R, C++ and the Armadillo C++ library.

The expected format of input files is that used by the mapping software polymapR (Peter M. Bourke, Geest, et al. 2018). If you have used other software for map construction and/or genotype calling, you will need to convert your input files to the correct format. There are already some tools available to do this (although it should be quite straightforward if you look at the expected format using the example data files provided here). For example, if you generated the phased linkage map using mappoly (Mollinari and Garcia 2019), you can convert your map to the required format using the convert_mappoly_to_phased.maplist function (check out the help using ?convert_mappoly_to_phased.maplist for details).

The genotypes can be either discrete or probabilistic. If the genotypes are discrete, they must be provided as a matrix of marker dosage scores (for example as provided in this tutorial in SNP_dosages.4x) with marker names as row-names, and individual names (including the two parents) as column-names. Checks are performed and non-numeric data is converted to numeric data if needed. Alternatively, probabilistic genotypes can be provided which can either be the direct output of the fitpoly (Voorrips, Gort, and Vosman 2011) function saveMarkerModels, or from some other polyploid-appropriate genotype-calling software. For example, if polyRAD (Clark, Lipka, and Sacks 2019) or updog (Gerard et al. 2018) were used for genotype calling, a conversion step is needed. Functions for doing this are provided by polymapR. Check out ?polymapR::convert_polyRAD or ?polymapR::convert_updog for details.

We run the function estimate_IBD using the phased linkage map phased_maplist.4x as generated by polymapR, in this case allowing for multivalents (bivalent_decoding = FALSE), as follows:

IBD_4x <- estimate_IBD(phased_maplist = phased_maplist.4x,
                       genotypes = SNP_dosages.4x,
                       ploidy = 4,
                       bivalent_decoding = FALSE,
                       ncores = 4)

Note that the default setting of bivalent_decoding is TRUE, as the function evaluates faster if only bivalents are considered (the difference in run-time becomes more pronounced for hexaploids). However, allowing multivalents in the model gives a more complete and correct picture.

Here we chose to enable parallel processing to speed up the calculations. Note that jobs are split across linkage groups, so with 5 linkage groups it would have made more sense to use 5 cores if they were available. Use parallel::detectCores() to display how many CPU cores are available, and use 1 or 2 less than this at the very most. In this dataset there are only 2 linkage groups, so no more than 2 cores are needed.

nc <- parallel::detectCores() - 1

IBD_4x <- import_IBD(method = "TO", #TO for TetraOrigin
                     folder = "TetraOrigin",
                     filename = paste0("TetraOrigin_Output_bivs_LinkageGroup",1:5,"_Summary"),
                     bivalent_decoding = TRUE)

Importing map data under description inferTetraOrigin-Summary,Genetic map of biallelic markers
Importing parental phasing under description inferTetraOrigin-Summary,MAP of parental haplotypes
Importing IBD data under description inferTetraOrigin-Summary,Conditonal genotype probability

IBD_4x <- import_IBD(method = "PO", #PO for PolyOrigin
                     folder = "PolyOrigin",
                     filename = "PolyOrigin_Output")

“Heuristic” method for IBD estimation

In cases where the results are needed quickly, or where there are very large numbers of markers, or for ploidy levels above 6, it is convenient to use the heuristic approach to IBD estimation. We do this using the estimate_IBD function as follows:

IBD_4x <- estimate_IBD(phased_maplist = phased_maplist.4x,
                       genotypes = SNP_dosages.4x,
                       method = "heur",
                       ploidy = 4)

Note that the attribute IBDtype is now haplotypeIBD, referring to the fact that these IBDs are the probabilities of inheriting each parental haplotype at a locus, as opposed to the probabilities of inheriting particular combinations of parental alleles at a locus (genotypeIBD). Although similar, certain downstream functions such as exploreQTL can only work with the former (although you will still be able to visualise QTL allele effects with the visualiseQTL function). The accuracy of these IBDs is generally lower than if method = "hmm" had been used (and therefore the power of subsequent QTL analyses will be somewhat reduced).

High marker densities

Particularly at higher ploidy levels (6x), it becomes computationally expensive to estimate IBDs using the hmm method. Our experience has shown that for hexaploids with a mapping population size of 400 or more individuals running on 4 cores, about 200 - 300 markers per chromosome can be accommodated on an “above-average” desktop computer (16 GB RAM). Running on a single core will reduce the committed memory, but the evaluation time will therefore be much longer. 250 markers per chromosome should be already enough to estimate IBDs with high accuracy - including extra markers over this amount will add incrementally less information (following the law of diminishing returns).

The function thinmap can be used to make a selection of markers that tries to maximise their distribution across the genome and across parental homologues. The argument bin_size is used to specify the size of the bins in which markers are selected - increasing this results in fewer markers being used. The ideal bin_size is 1 cM, although at higher ploidies and with higher marker densities, wider bins may be needed. The function is called as follows:

thinned_maplist.4x <- thinmap(maplist = phased_maplist.4x,
                              dosage_matrix = SNP_dosages.4x)

## 87 markers from a possible 93 on LG1 were included.
## 
## 89 markers from a possible 93 on LG2 were included.

The object thinned_maplist.4x can then be used in place of phased_maplist.4x in a call to estimate_IBD.

Interpolating IBDs

Regardless of how they were generated, IBD probabilities are estimated at all marker positions provided. When we perform an initial scan for QTL, it is often more efficient to look for QTL at a grid of evenly-spaced positions (such as at every 1 cM). This is because the genotype data has been replaced with multi-point IBD estimates at the marker positions. Information at e.g. ten different markers within a 0.5 cM window is virtually identical and therefore just one representative for this locus should approximate the full information, while reducing the number of statistical tests. This becomes particularly relevant when we generate significance thresholds using permutation tests later, which are computationally quite demanding.

It is therefore recommended to interpolate the probabilities using the spline_IBD function as follows:

IBD_4x.spl <- spline_IBD(IBD_list = IBD_4x, 
                         gap = 1) #grid at 1 cM spacing

Including experimental blocks

Note also that there is a separate column for “year” in this dataset. In general, only 1 level of experimental blocking is allowed. If there were replicated observations within blocks, then a more complicated nested block analysis would be needed. In such instances, it would probably be best to first derive best linear unbiased estimates (BLUEs) for the trait, with the full nested block structure provided. See the section on BLUEs below for more details. However, the QTL functions provided within polyqtlR are unable to inherit variance components from such an approach, and can only deal with the means in this case.

Running an initial QTL scan

Assuming we have a simple blocking structure as is the case here (single measurements per genotype repeated over 3 years) we can proceed with a preliminary QTL scan as follows:

qtl_LODs.4x <- QTLscan(IBD_list = IBD_4x,
                       Phenotype.df = Phenotypes_4x,
                       genotype.ID = "geno",
                       trait.ID = "pheno",
                       block = "year")

QTLscan by default uses an additive-effects model (Peter M. Bourke et al. 2021). Since package v.0.1.0 it is also possible to run a QTL scan using the full model as described in (Hackett, Bradshaw, and Bryan 2014), which uses the genotypic IBD probabilties rather than the haplotype IBD probabilities as regressors in the detection model. Although this risks over-fitting (as the model is over-parametrised, particularly if genotype probabilities were estimated allowing for double reduction which adds many extra genotype classes), it may be useful to also check and compare the results after significance thresholds have been calculated (see below). The full model can be accessed in the QTLscan function by setting argument allelic_interaction = TRUE.

The function gives a summary table of the trait over the different blocks, which can be useful for identifying possibly outlying observations. Because the function calls on the lm function from base R, missing values are not allowed in the trait scores. QTLscan performs some imputation if blocks are present, by estimating block effects and using them and the observed values to predict a missing phenotype. If this behaviour is not desired, then manually removing these datapoints beforehand is necessary. QTLscan can only work if the genotype identifiers (the offspring names) match completely between the trait data and the IBD data. If this is not the case, an error will result. The function reports the number of offspring with matching genotype and phenotype identifiers, so the user can record the population size that was used in the analysis.

Visualising QTL results

The object returned by QTLscan is a list with a number of items, including:

• QTL.res A matrix of QTL results

• Perm.res The results of a permutation test, if run (see below). In our example here, this field is NULL, i.e empty.

• Residuals Vector of residual phenotypes after blocks or co-factors have been fitted.

• Map The linkage map upon which the IBD probabilities were based.

More details on the output can be found by running ?QTLscan. There are a number of functions in the package for visualising the results. One option is to plot per linkage group, which can be achieved using the plotQTL function and grid layout (layout = 'g'). An example call would be as follows:

plotQTL(LOD_data = qtl_LODs.4x,
        layout = "g",
        colour = "darkgreen")

As before, details of the options available for this function can be accessed using ?plotQTL.

An alternative is to plot the results of the whole-genome scan on a single track, which since package v.0.1.0 is the default plotting behaviour of the plotQTL function:

plotQTL(LOD_data = qtl_LODs.4x,
        colour = "dodgerblue")

plotQTL can also visualise the results of a number of analysis together. By default, the superimposed plots are re-scaled so that threshold lines overlap (argument rescale = TRUE). This is covered in more detail in the section on co-factors.

Significance thresholds

Although it already appears like we have an interesting peak on LG 1, we are not sure whether this region is truly associated with our trait or not. One method of establishing whether the association is significant is to run a permutation test. Currently the package has two options. The first uses the full structure of your data as provided. However, it is much faster to run a permutation test when your data has a simple structure, with no experimental blocks or no genetic co-factors. To deal with experimental blocks, you should estimate BLUEs first before running the permutation test. For genetic co-factors, you can first run an analysis with the co-factor(s) included, and then use the residuals as your phenotype. However, it is also possible (but not recommended) to simply run a permutation test without these steps, using the QTLscan function as follows:

qtl_LODs.4x <- QTLscan(IBD_list = IBD_4x,
                       Phenotype.df = Phenotypes_4x,
                       genotype.ID = "geno",
                       trait.ID = "pheno",
                       block = "year",
                       perm_test = TRUE,
                       ncores = nc)

If we now plot the results, the threshold line is automatically included:

plotQTL(LOD_data = qtl_LODs.4x,
        colour = "dodgerblue")

## LG_names not detected in LOD_data. Re-creating LG names for plot...

However, if you ran this step you will have noticed that it is quite slow. Therefore, it is recommended to first estimate a mean phenotypic value per individual over years before running the permutation test, which we do here using best linear unbiased estimates in the nlme package which fits a linear mixed model with genotype as fixed effect.

BLUEs

We can estimate the best linear unbiased estimates using the BLUE function as follows:

blues <- BLUE(data = Phenotypes_4x,
              model = pheno~geno,
              random = ~1|year,
              genotype.ID = "geno")

The resulting object is a data-frame with columns “geno” and “blue”, which we need in our revised call to QTLscan:

qtl_LODs.4x <- QTLscan(IBD_list = IBD_4x,
                       Phenotype.df = blues,
                       genotype.ID = "geno",
                       trait.ID = "blue",
                       perm_test = TRUE,
                       ncores = nc)

If we now plot the results, the threshold line is automatically included:

plotQTL(LOD_data = qtl_LODs.4x,
        colour = "dodgerblue")

Note that the LOD scores have now increased, but the significance threshold has also increased. This has to do with how the LOD score is estimated, using a ratio of the RSS from your fitted model at a putative QTL position with that of the NULL model with no QTL. In the case of blocks, the NULL model already contains your blocking term, and so a portion of the phenotypic variance is already accounted for when calculating the LOD score. When using BLUEs as your phenotype, you first correct your phenotypes for block effects, but the NULL model now has no explanatory variables, which artificially inflates your LODs. However, a permutation test will account for these sorts of discrepancies, and the region of the genome that lies above the threshold should be the same.

Adding genetic co-factors

In many situations we may wish to include already-identified QTL positions as fixed genetic effects in a QTL search. This can help to reduce the amount of unwanted background variance in the phenotypic data which can lead to greater power to detect QTL at other positions. Another possible reason could be to test whether multiple peaks on a single linkage group are independent. There are two options for including co-factors in the package. We could use the automated function check_cofactors which performs all steps automatically, or we can run a manual co-factor analysis using QTLscan by specifying the cofactor_df argument.

findPeak(qtl_LODs.4x, linkage_group = 1)

##   chromosome          marker position      LOD  thresh
## 9          1 LG1_12.3_0x1_h5     12.3 12.29909 5.27213

## [1] 12.3

qtl_LODs.4x_cofactor <- QTLscan(IBD_list = IBD_4x,
                                Phenotype.df = Phenotypes_4x,
                                genotype.ID = "geno",
                                trait.ID = "pheno",
                                block = "year",
                                cofactor_df = data.frame("LG" = 1,
                                                         "cM" = 12.3),
                                perm_test = FALSE,
                                ncores = nc)#nc is the number of cores, defined earlier

new_pheno <- qtl_LODs.4x_cofactor$Residuals
head(new_pheno)

##    geno       Pheno Block
## 1 F1_01  3.06894507  2014
## 2 F1_02  4.80132691  2014
## 3 F1_03  0.09938343  2014
## 4 F1_04 -1.65306259  2014
## 5 F1_05  0.96184614  2014
## 6 F1_06 -5.48341785  2014

colnames(new_pheno)

## [1] "geno"  "Pheno" "Block"

blues <- BLUE(data = new_pheno,
              model = Pheno~geno,
              random = ~1|Block,
              genotype.ID = "geno")

qtl_LODs.4x_cofactor <- QTLscan(IBD_list = IBD_4x,
                                Phenotype.df = blues,
                                genotype.ID = "geno",
                                trait.ID = "blue",
                                perm_test = TRUE,
                                ncores = nc)

plotQTL(LOD_data = qtl_LODs.4x_cofactor,
        colour = "red")

blues <- BLUE(data = Phenotypes_4x,
              model = pheno~geno,
              random = ~1|year,
              genotype.ID = "geno")

QTLmodel <- check_cofactors(IBD_list = IBD_4x,
                            Phenotype.df = blues,
                            genotype.ID = "geno",
                            trait.ID = "blue",
                            LOD_data = qtl_LODs.4x,
                            ncores = nc)


QTLmodel

##   LG   cM deltaLOD CofactorID
## 1  1 12.3 7.026958         NA

PVE(IBD_list = IBD_4x,
    Phenotype.df = blues,
    genotype.ID = "geno",
    trait.ID = "blue",
    QTL_df = QTLmodel)

## Warning in PVE(IBD_list = IBD_4x, Phenotype.df = blues, genotype.ID = "geno", : Maximal QTL model has a single QTL!
## A single QTL scan using function QTLscan provides comparable results!

## $`1 QTL model`
##    Q1 
## 67.79

Comparing multiple analyses

It is often convenient to plot the results of different analyses together. For example comparing the results of a single marker analysis with an IBD-based approach, or plotting an initial QTL scan versus the results when a set of genetic co-factors were also included. In cases where significance thresholds were determined in separate analyses, this will soon lead to cluttered-looking plots with multiple significance thresholds on top of each other. One possibility is to disable the visualisation of significance thresholds when rendering each plot. However, it is often of interest to compare the results in terms of the presence (or absence) of significant QTL peaks. For this reason the default behaviour of plotQTL when multiple results are passed to LOD_data is to re-scale the overlaid plots so that the significance thresholds for each set of results overlaps perfectly. In this case, LOD_data should be a list of separate QTL results, each of which is the output of a separate call to the function QTLscan or singleMarkerRegression. Re-scaling can be disabled by setting rescale = FALSE.

We can compare the results obtained earlier from the original and co-factor analyses as follows:

plotQTL(LOD_data = list(qtl_LODs.4x,
                        qtl_LODs.4x_cofactor))

Compare QTL results with offspring

We might also like to see whether the QTL results tally with the allelic conformation of the population if sorted by phenotype. For this, we use the display_by = "phenotype" argument in the visualiseHaplo function. First, we look at the distribution of (BLUE) phenotypes of the population:

hist(x = blues$blue)

Suppose we are only interested in offspring with a trait value higher than 44. We can choose to only display the haplotypes of these offspring as follows:

visualiseHaplo(IBD_list = IBD_4x,
               display_by = "phenotype",
               Phenotype.df = blues,
               genotype.ID = "geno",
               trait.ID = "blue",
               pheno_range = c(44,max(blues$blue)),
               linkage_group = 1,
               multiplot = c(3,3))

All these offspring seem to carry homologue 1 as anticipated, although certain individuals (like F1_24) possess recombinations in the region of the QTL, which might be of interest for further fine-mapping work. However, the main point here is that the results seem to be consistent with the predictions from the BIC analysis.

Meiosis report

It is interesting to look at what went on in parental meiosis in the population, as it may give an indication on e.g. the pairing behaviour of chromosomes during meiosis. For many polyploid crops this is not necessarily a well-understood process, so we provide some tools to help give insights into whether your species of interest is autopolyploid, allopolyploid or something in between (“segmental allopolyploid”).

To this end, we run the function meiosis_report as follows:

mr <- meiosis_report(IBD_list = IBD_4x)

Preferential pairing

Apart from visualising the deviations (in counts) from a purely random bivalent pairing model, the function meiosis_report also analyses all predicted bivalent pairings and performs an uncorrected χ² test to detect deviations from a polyomic model. If multivalents were also included during IBD estimation these are also counted and recorded, but do not contribute to the χ² test or output plots. Note that the colour of the output plot conveys information about the results of the χ² test: green signifying no strong deviation from a polysomic model, while red colours indicating a more severe deviation from a polysomic model. In this example, there is no suggestion of anything but polysomic inheritance, associated with autopolyploidy.

An alternative approach to visualise the pairing behaviour of parental homologues is provided by the function visualisePairing. This function displays deviations from random expectations as before, but using line thickness to represent the magnitude of the deviation (by default, an excess of pairing counts is coloured red and a lack of counts coloured blue). One argument needed in order for the function to work is datawidemax, which is the data-wide maximum deviation. This is then given the line width as specified by max.lwd, and all other deviations are scaled accordingly. See ?visualisePairing for more details.

The function can be run using the output of meiosis_report as follows:

par(mfrow = c(1,2))
visualisePairing(meiosis_report.ls = mr,
                 parent = "P1",
                 datawidemax = 5)

Recombination landscape

It can be very informative to perform an analysis of all the predicted recombinations (from cross-over events), as these can indicate issues such as poor map order, problematic individuals (carrying too many predicted recombinations) or errors in marker phase. The position of these can also be biologically relevant, indicating recombination hot-spots or cold-spots. They also provide a means to correct genotyping errors, if we take the view that genotyping errors are a more likely cause of multiple recombination events in an individual rather than cross-overs (which they usually are). In short, an analysis of the recombination landscape offers a lot of scope for improved understanding on many levels.

Within polyqtlR a number of tools are provided for this. We start by looking at the function count_recombinations, which generates a summary of the recombination events per linkage group and per individual:

recom.ls <- count_recombinations(IBD_4x)

Note that this function has only been implemented in the context of bivalent pairing - a likely future update is to also consider and count recombination events from multivalents.

The output is a rather long list giving, per linkage group and per individual the probability of the most likely pairing configuration and the most likely positions of the recombination break-point (taken as the midpoint of the two bounding positions which define the interval). Although we do not know the precise break-point around cross-over events, this value can be taken as indicative of the likely location. To then make sense of this list we call on the visualisation function plotRecLS (plot the Recombination LandScape) to plot the results and return some summary information from the plots:

layout(matrix(c(1,3,2,3),nrow=2)) #make tidy layout
RLS_summary <- plotRecLS(recom.ls) #capture the function output as well

The green lines show a fitted spline to the recombination count data - which are also returned by the function. Indeed, if you explore the output saved in RLS_summary you should see that there are two items: $per_LG and per_individual. In the per-linkage group item, recombination rates that have been interpolated using the splines (rec.rates) are given (per linkage group). In the per-individual item, a named vector of the genome-wide counts of recombinations is given, which can be useful in identifying offspring that perhaps should not have been included in the mapping.

It is also of interest to overlay the predicted recombinations with the visualised haplotypes. This can be done with the visualiseHaplo function, similar to what we did earlier but now adding the argument recombination_data to the call:

visualiseHaplo(IBD_list = IBD_4x,
               display_by = "name",
               linkage_group = 1,
               select_offspring = colnames(SNP_dosages.4x)[3:11], #cols 1+2 are parents
               multiplot = c(3,3), #plot layout in 3x3 grid
               recombination_data = recom.ls) 

It is interesting to note the problems with reconstructing the inheritance probabilities of individual F1_01 - compare this to the reconstruction earlier when we also allowed for multivalents and double reduction. F1_04 has also been reconstructed differently, but here the bivalent-only based prediction does not seem to have caused serious conflicts.

Searching for recombinant individuals

We have analysed the recombination landscape, but suppose for the sake of argument we were interested in all offspring that had inherited homologues 1 and 3 from linkage group 1 in the region of our detected QTL. Although we only found evidence for a positive effect of homologue 1, we could imagine a different dataset where there was a second QTL on homologue 3 some distance away. If we wanted to stack these QTL, then recombinant offspring could offer the possibility of stacking these QTL in coupling linkage. We can search for these individuals using the visualiseHaplo function.

We display the offspring by “name” as before. However, we now also provide details of the recombinant_scan homologues that we are interested in. The output of this initial function call will be used in the subsequent step, so we save the output as “visHap1”. We will also use multi-plotting for convenience. Note that the plots have been suppressed here (as there are too many of them to begin with):

visHap1 <- visualiseHaplo(IBD_list = IBD_4x,
                          display_by = "name",
                          select_offspring = "all",
                          linkage_group = 1,
                          cM_range = c(1.56,50),
                          recombinant_scan = c(1,3),
                          multiplot = c(4,2))

We now look at the output of the function:

visHap1

## $recombinants
## [1] "F1_02" "F1_03" "F1_16" "F1_22" "F1_26" "F1_41"
## 
## $allele_fish
## NULL

Six offspring have been identified as having a recombination in this region. There may be others, as this function is only searching among the IBD results for potential recombinants and may miss some with less clear features. However, six is already enough to continue with. We pass the names of these individuals to the function again (in place of “all” offspring) using the select_offspring argument, and specifying the individuals mentioned in visHap1$recombinants, and specifying the desired region of linkage_group 1 using cM_range as follows:

visualiseHaplo(IBD_list = IBD_4x,
               display_by = "name",
               select_offspring = visHap1$recombinants,
               linkage_group = 1,
               cM_range = c(1.56,50),#note: start position = first marker @1.56 cM
               recombinant_scan = c(1,3),
               multiplot = c(3,2),
               recombination_data = recom.ls) # we can add this track for clarity

In each of these individuals we can see evidence for a recombination between homologues 1 and 3 in the specified region.

Note - currently this process takes two steps, which could be confusing. A planned update will allow this to be done in a single step, using the recombination_data instead.

Fishing for specific alleles

The other output associated with the function is $allele_fish, which allows us to fish out individuals that carry a particular allele or combinations of alleles around a particular locus. Suppose we were just interested in the offspring that carried homologue 1 of chromosome 1 around the QTL, so in the region up to 15 cM say. We can “fish out” these individuals as follows:

visHap2 <- visualiseHaplo(IBD_list = IBD_4x,
                          display_by = "name",
                          select_offspring = "all",
                          linkage_group = 1,
                          cM_range = c(1.56,15),
                          allele_fish = 1,
                          multiplot = c(4,2)) 

Note that these individuals must carry homologue 1 for almost all of the specified segment - therefore specifying a longer segment could also result in fewer individuals (as they are more likely to carry recombined chromosomes). As before, we can inspect the output of this function call as follows:

visHap2

## $recombinants
## NULL
## 
## $allele_fish
## $allele_fish$h1
##  [1] "F1_02" "F1_08" "F1_10" "F1_15" "F1_17" "F1_19" "F1_23" "F1_31" "F1_32"
## [10] "F1_34" "F1_35" "F1_40" "F1_43" "F1_44" "F1_46" "F1_49"

Here there are 16 offspring carrying that particular segment of homologue 1. We can visualise them using a call to the visualiseHaplo function as before:

visualiseHaplo(IBD_list = IBD_4x,
               display_by = "name",
               select_offspring = visHap2$allele_fish$h1,
               linkage_group = 1,
               cM_range = c(1.56,15),
               multiplot = c(4,4))

If all went well, these offspring should all possess the homologue 1 allele in this particular region. Examining the plots shows this to be the case.