In this tutorial, we will go step by step through the basic steps of performing linkage analysis in polyploids using polymapR. Most of the theory behind the functions described here can be found in the Bioinformatics publication of Bourke et al. (2018) at https://doi.org/10.1093/bioinformatics/bty371. The example dataset is a simulated dataset of an hypothetical tetraploid outcrossing plant species with five chromosomes. Chromosomes have random pairing behaviour. The population consists of 207 individuals originating from an F1 cross between two non-inbred parents. Frequently occurring data imperfections like missing values, genotyping errors and duplicated individuals are taken into account.
7. Quality checks on marker data
A number of quality checks are recommended before we begin mapping. For example, we do not want to include markers which have too many missing values (“NA” scores) as this may be an indication that the marker has performed poorly and may have many errors associated with it. We would also like to screen out individuals that have many missing values as this might be an indication that the DNA quality of that sample was poor, or duplicate individuals which give no extra information and should probably also be removed / combined. Because there are no fixed rules for how stringent we should be with quality checks (5% missing values ok? 10% missing values ok? 20%?), some visual aids are produced in order to help with these decisions.
7.1 Missing value rate per marker
We first screen the marker data for missing values using screen_for_NA_values. This function allows the user to remove rows or columns with a certain rate of missing values from the dataset. The dosage matrix, which is specified as segregating_data contains markers in rows and individuals in columns. In R, rows are usually specified by the number 1, and columns by 2. The argument margin lets you specify by margin number whether you want to screen markers (margin = 1) or individuals (margin = 2) If the argument cutoff is set to NULL, it takes user input on the rate of markers to be screened. In the example below cutoff is specified, but because decisions are made after data inspection, we recommend to decide on the cutoff value after inspection of the data by specifying cutoff = NULL. A threshold of 10% missing values is sensible, as a compromise between high-confidence marker data versus keeping a large proportion of the markers for mapping and QTL analysis. In general, the most problematic markers tend to have higher rates of missing values so any reasonable threshold here should be fine. If problems are encountered later in the mapping, it might be worthwhile to re-run this step with a more stringent (e.g. 5%) threshold.
screened_data <- screen_for_NA_values(dosage_matrix = segregating_data,
margin = 1, # margin 1 means markers
cutoff = 0.10,
print.removed = FALSE)
##
## Number of removed markers : 1458
##
##
## There are now 1415 markers leftover in the dataset.
7.2 Missing value rate per individual
As in step 7.1 we can also visualise the rate of missing values per individual, which can help us to make a decision about whether certain individuals should be removed from the dataset before commencing the mapping. In order to screen the individuals for missing values we enter:
screened_data2 <- screen_for_NA_values(dosage_matrix = screened_data,
cutoff = 0.1,
margin = 2, #margin = 2 means columns
print.removed = FALSE)
##
## Number of removed individuals : 2
##
##
## There are now 205 individuals leftover in the dataset.
Again, it is not clear what an acceptable threshold for the rate of missing values should be, but we recommend a cut-off rate of 0.1 – 0.15 (10-15%) missing values allowed.
7.3 Duplicate individuals
It is also desirable to be able to check whether there are any duplicated individuals in the population as can sometimes happen (or due to mix-ups in the DNA preparation and genotyping). The function screen_for_duplicate_individuals does this for you. Note that a duplicate can remain in the dataset, but they add no further information for mapping or QTL analysis, and may even distort the analysis (by giving an inaccurate population size). As for the function screen_for_NA_values, cutoff can be set to NULL to accept user input. Especially for screen_for_duplicate_individuals it makes sense to make a decision on the cut off after inspection of the data and figures. To screen the data for duplicate individuals, enter:
screened_data3 <- screen_for_duplicate_individuals(dosage_matrix = screened_data2,
cutoff = 0.95,
plot_cor = T)
##
## Combining F1_060 & F1_081 into F1_060
##
## Combining F1_046 & F1_085 into F1_046
## Warning in screen_for_duplicate_individuals(dosage_matrix = screened_data2, :
## Multiple duplicates of single genotype identified at this threshold. Attempting
## to merge...
##
## Combining F1_032 & F1_100 into F1_032
##
## Combining F1_032 & F1_114 into F1_032
##
## Combining F1_079 & F1_101 into F1_079
##
## Combining F1_106 & F1_113 into F1_106
##
## Combining F1_112 & F1_199 into F1_112
##
## #### 7 individuals removed:
##
## _ _ _ _
## [1,] "F1_081" "F1_085" "F1_100" "F1_114"
## [2,] "F1_101" "F1_113" "F1_199" ""
## |_ |_ |_ |_ |
## |:------|:------|:------|:------|
## |F1_081 |F1_085 |F1_100 |F1_114 |
## |F1_101 |F1_113 |F1_199 | |
Here we see that there were a number of duplicate pairs of individuals identified (red points in the output graph). Each duplicate pair is merged into a consensus individual (conflicts are made missing), keeping the name of the first individual. This is summarised in the output as well.
7.4 Duplicated markers
Finally, we can also check for duplicated markers (those which give us no extra information about a locus) and remove these from the dataset before mapping, using the screen_for_duplicate_markers function. This function returns a list with 2 items - the filtered_dosage_matrix which we are interested in now, and bin_list, which we are possibly interested in again later. It is a good idea to save the output of this function, as you might like to add back the duplicate markers after mapping.
screened_data4 <- screen_for_duplicate_markers(dosage_matrix = screened_data3)
##
|
| | 0%
|
|======================================================================| 100%
## Marked and merged 3 duplicated markers
Note that increasingly, there are far too many markers generated in comparison to the population size. It is possible to estimate the genetic resolution of the population you are working with, and from that estimate the expected bin-size in your dataset, assuming your markers are approximately randomly-distributed across the genome (this may not be a good assumption in your case, but as a rough guide this function may still be useful). Considering only bivalent pairing, in a tetrasomic species we can assume approximately four cross-over recombinations occurred, taking the rule-of-thumb guide of one cross-over recombination per bivalent, which has been observed in autotetraploid Arabidopsis arenosa for example [2].
We can try to estimate the average bin size if we run the previous function and specifying the argument estimate_bin_size = TRUE. The function then does a quick calculation on how many markers would be expected in an average marker bin (under the tenuous assumptions of uniform distributions of both markers and cross-over break-points). In a tetraploid, each offspring would be expected to have 4 cross-over break-points per chromosomal linkage group (2 maternal and 2 paternal), so across a population of N individuals there are approximately 4N breaks per linkage group, and 4NL such breaks across L linkage groups. Supposing there are in total y distinct marker types (like 1x0, 1x1, 2x0, 0x1 etc), then your M markers can be split into M/4NLy segregation patterns (assuming no errors), which is an approximate estimate for the average bin size. In datasets with very large numbers of markers (e.g. genotyping using NGS), it may be worthwhile to only include markers with a minimum number of duplicates present, and that minimum number can be informed by the calculation above. In this case, supposing we demand at least 7 duplicates (so a minimum bin size of 6), we would do something like the following:
reliable.markers <- names(which(sapply(screened_data4$bin_list,length) >= 6))
reliable_data <- screened_data4$filtered_dosage_matrix[reliable.markers,]
Of course, this assumes we have a very large number of markers to choose from, as this is a very severe demand in a more traditional, well-balanced dataset. It also presumes there are no genotyping errors. A genotyping error will make duplicate markers appear to be different even when they are not (from a genetic mapping perspective). For now we are ignoring this issue, but may return to it in future updates.
In the current data-set there is quite a good balance between the number of markers and the number of individuals in the mapping population, and so we extract the filtered_dosage_matrix that has the duplicate markers merged:
filtered_data <- screened_data4$filtered_dosage_matrix
Before we start the mapping, it is a good idea to run the summary function again, this time on the filtered_data object, so we have a clear record of the numbers of markers going forward to the mapping stage.
pq_screened_data <- parental_quantities(dosage_matrix = filtered_data)
## | | 1x0| 2x0| 0x1| 1x1| 2x1| 0x2| 1x2| 2x2| 1x3|
## |:---------|---:|---:|---:|---:|---:|---:|---:|---:|---:|
## |frequency | 175| 111| 194| 301| 183| 101| 198| 88| 61|
8. Simplex x nulliplex markers – defining chromosomes and homologues
The first step in actually mapping the markers is to define linkage groups for the markers. We use the simplex x nulliplex markers to define the linkage groups (chromosomes) and following that, the homologues (4 per linkage group are expected in a tetraploid). We use the LOD scores to cluster markers, which increases dramatically for tight coupling-phase estimates.
To estimate the recombination frequencies, LOD and phasing between all marker pairs of P1 enter:
SN_SN_P1 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
target_parent = "P1",
other_parent = "P2",
ploidy = 4,
pairing = "random"
)
A note on arguments and multi-core processing
There is much more to the linkage function. Checkout ?linkage to see the arguments that can be passed to it. One important feature is that it can run in parallel on multiple cores. This can be very time-efficient if you are working with a large number of markers. The number of cores to use can be specified using the ncores argument. Before doing so, first check the number of cores on your system using parallel::detectCores(). Use at least one core less than you have available to prevent overcommitment.
It is a good idea to first visualise the relationship between the recombination frequency and LOD score, which gives us some insight into the differences between the information content of the different phases, as well as highlighting possible problems in the data. This can be done by using the plot_linkage_df function:
plot_linkage_df(SN_SN_P1, r_max = 0.5)
In some situations, it is possible to identify outlying datapoints on this plot, which do not fit the expected relationship. We have implemented a function to check this for the case of pairs of simplex x nulliplex markers, as these are the markers that we use in the subsequent steps of clustering and identifying chromosome and homologue linkage groups. If you are unsure, you can run the SNSN_LOD_deviation function anyway to see whether your data fits the expected pattern:
P1deviations <- SNSN_LOD_deviations(linkage_df = SN_SN_P1,
ploidy = 4,
N = ncol(filtered_data) - 2, #The F1 population size
alpha = c(0.05,0.2),
plot_expected = TRUE,
phase="coupling")
Here we see there are a few “outliers” in the coupling pairwise data (identified by the stars) from the expected relationship between r and LOD (the bounds of which are shown by the dotted lines - these can be widened as necessary using the alpha argument for the upper and lower tolerances around the “true” line, usually in a trial and error approach). In our example dataset, these highlighted points are hardly outliers and we can proceed. In other datasets, you may want to remove these starred pairwise estimates before continuing with marker clustering, as it is possible they may cause unrelated homologue clusters to clump together and prevent a clear clustering structure from being defined. To do so, remove them from the linkage data frame as follows:
SN_SN_P1.1 <- SN_SN_P1[SN_SN_P1$phase == "coupling",][-which(P1deviations > 0.2),]
We are now in a position to cluster the marker data into chromosomes using the cluster_SN_markers function. This function clusters markers over a range of LOD thresholds. In the example below this is over a range from LOD 3 to 10 in steps of 1 (LOD_sequence = c(3:10)). If the argument plot_network is set to TRUE this function will plot the network of the lowest LOD threshold (LOD 3 in this case), and overlays the clusters of the higher LOD scores. This setting is not recommended with large number of marker pairs. Your computer will most likely run out of memory.
P1_homologues <- cluster_SN_markers(linkage_df = SN_SN_P1,
LOD_sequence = seq(3, 10, 0.5),
LG_number = 5,
ploidy = 4,
parentname = "P1",
plot_network = FALSE,
plot_clust_size = FALSE)
## Total number of edges: 894
The function cluster_SN_markers returns a list of cluster data frames, one data frame for each LOD threshold. The first plot shows two things - on the normal y-axis, the number of clusters over a range of LOD thresholds (x-axis) - and on the right-hand y-axis, the number of markers that do not form part of any cluster, which is shown by the blue dashed line. In the cluster_SN_markers function, we specify the ploidy (4 for a tetraploid) and the expected number of chromosomes (LG_number = 5) and therefore there are 4 x 5 = 20 homologue clusters expected. If we examine the second plot, we see how these clusters stay together or fragment as the linkage stringency (LOD score) increases. Here, we see clearly that at LOD 3.5, there are 5 clusters which all eventually divide into 4 sub-clusters each. This is the ideal scenario - we have identified both the chromosomal linkage groups and the homologues using only the simplex x nulliplex marker data. Note that it is the coupling-phase estimates which define the homologue clusters and the repulsion-phase estimates which provide the associations between these coupling-clusters (homologues). In practice, things might not always be so orderly. In noisy datasets with higher numbers of genotyping errors, or at higher ploidy levels, the repulsion-phase linkages between homologues will be lost among false positive linkages between un-linked markers. This means that all markers will likely clump together at lower LOD scores.
You can manually check how many clusters there are at a specific LOD score, and how many markers are in a cluster:
P1_hom_LOD3.5 <- P1_homologues[["3.5"]]
t <- table(P1_hom_LOD3.5$cluster)
print(paste("Number of clusters:",length(t)))
## [1] "Number of clusters: 5"
t[order(as.numeric(names(t)))]
##
## 1 2 3 4 5
## 37 36 29 31 41
P1_hom_LOD5 <- P1_homologues[["5"]]
t <- table(P1_hom_LOD5$cluster)
print(paste("Number of clusters:",length(t)))
## [1] "Number of clusters: 20"
t[order(as.numeric(names(t)))]
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 13 9 10 4 12 7 9 8 3 13 6 7 4 10 6 11 17 6 6 12
Based on the plots (or the output as shown) you should be able to decide which LOD threshold best separates the data. In this example, at a LOD of 3.5 we identify chromosomes and anywhere between LOD 4.5 and LOD 7 we have identified our homologues (4 x 5 = 20 clusters). In this parent (P1) we have therefore solved the marker clustering problem without any difficulty. It remains to gather this into a structure that can be used at later stages. To do this, we need two LOD values which define chromosomal groupings (here, LOD_chm = 3.5) and homologue groupings (here, LOD_hom = 5). We can then use the define_LG_structure function with these arguments as follows (note the argument LG_number refers to the expected number of chromosomes for this species):
LGHomDf_P1 <- define_LG_structure(cluster_list = P1_homologues,
LOD_chm = 3.5,
LOD_hom = 5,
LG_number = 5)
##
## #### SxN Marker(s) lost in clustering step at LOD 5:
##
## |_ |
## |:-----------|
## |Ls_rs099908 |
## SxN_Marker LG homologue
## 1 Ac_ns035629 1 1
## 3 Ac_ns045766 1 1
## 7 Ac_ns096474 1 1
## 22 Ac_ts072686 1 1
## 24 Ac_ws013795 1 1
## 42 Ap_rs075640 1 1
However, you may find that it is not possible to identify chromosomal groups as simply as this. We therefore recommend that homologue groups be first identified, and then another marker type be used to provide associations between these clusters, and hence help form chromosomal linkage groups. This can be done as follows:
Generate homologue clusters using the cluster_SN_markers function at a high LOD.
Use the Duplex x Nulliplex or Simplex x Simplex marker information to put the clusters back together into chromosomal groups.
It makes sense to use only SxN marker pairs in coupling phase for the homologue clustering (since these associations define the homologues). Note that the most likely phase per marker pair (coupling or repulsion) was already estimated by our previous call to the linkage function.
SN_SN_P1_coupl <- SN_SN_P1[SN_SN_P1$phase == "coupling",] # select only markerpairs in coupling
P1_homologues_1 <- cluster_SN_markers(linkage_df = SN_SN_P1_coupl,
LOD_sequence = c(3:12),
LG_number = 5,
ploidy = 4,
parentname = "P1",
plot_network = FALSE,
plot_clust_size = FALSE)
## Total number of edges: 682
In the dot plot above, we see that all chromosomal associations have now been removed with the repulsion information.
A note on nested lists
In the polymapR package, output is often given as a nested list. Some basic information on lists is available on e.g. r-tutor. We use nested lists because we are often dealing with hierarchical data: marker data > homologue > linkage group > organism. Some basic R knowledge is needed to retrieve, visualize or filter data from such a nested list. However, if you would like to visualize it without R (e.g. in your favourite spreadsheet program), you can choose to export a nested list using write_nested_list. It will keep its hierarchical structure by writing it to directories that have the same structure as the the nested list. Here’s an example: write_nested_list(list = P1_homologues, directory = "P1_homologues")
Next, we use linkage to the DxN markers to provide bridging linkages between homologue clusters, using the bridgeHomologues function. The first time we run this function, we use a LOD threshold of 4 for evidence of linkage between a SxN and a DxN marker. However, there may also be some false positive linkages at this level. It might therefore be necessary to run the function again at a higher LOD (say LOD 8) to avoid such false positives. In doing so, we might lose true linkage information. Therefore, it’s a good idea to run the function at LOD 4 first, then increase to a higher LOD if necessary, and build up a picture of which clusters form part of the same chromosome. First calculate linkage between SxN and DxN markers:
SN_DN_P1 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
markertype2 = c(2,0),
target_parent = "P1",
other_parent = "P2",
ploidy = 4,
pairing = "random")
Then, use bridgeHomologues to use these SxN - DxN linkages to find associations (“bridges”) between homologue clusters. Note that we must provide a particular division of the Simplex x Nulliplex data - in this example, we use the division at LOD 6 which gave us 20 clusters (our 20 putative homologues):
LGHomDf_P1_1 <- bridgeHomologues(cluster_stack = P1_homologues_1[["6"]],
linkage_df = SN_DN_P1,
LOD_threshold = 4,
automatic_clustering = TRUE,
LG_number = 5,
parentname = "P1")
This results in an ideal situation of 5 linkage groups each with exactly 4 homologues:
table(LGHomDf_P1_1$LG, LGHomDf_P1_1$homologue)
##
## 1 2 3 4
## 1 13 9 4 10
## 2 12 7 9 8
## 3 3 13 5 7
## 4 4 10 6 11
## 5 17 12 6 6
However, clustering of P2 is a bit more difficult. The following re-does everything for P2:
SN_SN_P2 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
target_parent = "P2",
other_parent = "P1",
ploidy = 4,
pairing = "random"
)
SN_SN_P2_coupl <- SN_SN_P2[SN_SN_P2$phase == "coupling",] # get only markerpairs in coupling
P2_homologues <- cluster_SN_markers(linkage_df = SN_SN_P2_coupl,
LOD_sequence = c(3:12),
LG_number = 5,
ploidy = 4,
parentname = "P2",
plot_network = FALSE,
plot_clust_size = FALSE)
## Total number of edges: 743
SN_DN_P2 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
markertype2 = c(2,0),
target_parent = "P2",
other_parent = "P1",
ploidy = 4,
pairing = "random")
LGHomDf_P2 <- bridgeHomologues(cluster_stack = P2_homologues[["6"]],
linkage_df = SN_DN_P2,
LOD_threshold = 4,
automatic_clustering = TRUE,
LG_number = 5,
parentname = "P2")
table(LGHomDf_P2$LG,LGHomDf_P2$homologue)
##
## 1 2 3 4 5
## 1 11 9 3 5 0
## 2 13 13 7 11 0
## 3 13 11 11 8 2
## 4 9 9 13 11 0
## 5 10 10 6 9 0
Occasionally linkage groups possesses fewer than four “homologues”. This can result in a loss of phase information in later steps, but will not prevent us from proceeding with the mapping. However, it is also possible that more than 4 homologues are grouped together, likely the result of having a number of fragmented homologues. However, as we don’t know which clusters are actual homologues, and which are fragments, we must examine the affected linkage group in more detail. polymapR currently offers two methods for this.
8.1 Using the function overviewSNlinks
Run overviewSNlinks as follows:
overviewSNlinks(linkage_df = SN_SN_P2,
LG_hom_stack = LGHomDf_P2,
LG = 3,
LOD_threshold = 0)
In cases of fragmented homologues, it may be possible to identify which fragments belong together by examining the phase of intra-cluster marker pairs. Coupling-phase (green) linkage between fragments can be an indication that separate homologue fragments belong together. This can be then be used to merge these homologues (specified in the argument mergeList - use ?merge_homologues for more information). For example, if we wanted to merge homologues 4 and 5 of linkage group 3 as the above data suggests, we would use the following:
LGHomDf_P2_1 <- merge_homologues(LG_hom_stack = LGHomDf_P2,
ploidy = 4,
LG = 3,
mergeList = list(c(4,5)))
## #### Number of markers per homologue:
##
## | | 1| 2| 3| 4| 5|
## |:--|--:|--:|--:|--:|--:|
## |1 | 11| 13| 13| 9| 10|
## |2 | 9| 13| 11| 9| 10|
## |3 | 3| 7| 11| 13| 6|
## |4 | 5| 11| 10| 11| 9|
8.2 Using the function cluster_per_LG
cluster_per_LG does not rely on markers in repulsion, but uses variation in LOD thresholds between markers in coupling. Run it as follows:
cluster_per_LG(LG = 3,
linkage_df = SN_SN_P2[SN_SN_P2$phase == "coupling",],
LG_hom_stack = LGHomDf_P2,
LOD_sequence = c(3:10), # The first element is used for network layout
modify_LG_hom_stack = FALSE,
network.layout = "stacked",
nclust_out = 4,
label.offset=1.2)
## Total number of edges: 185
The network layout is based on the first LOD score you provide in LOD_sequence (which is LOD = 3 in this example). The darker the background, the more consistent the clustering is over the range of LOD scores. Change the LOD_sequence if the network is not as required (4 clusters in the case of a tetraploid) and re-run. If the network looks as you would like it to, use the cluster_per_LG function but with a LOD_sequence of length 1 (the first LOD score, so again LOD = 3 in this example), and specify modify_LG_hom_stack = TRUE to make a modification to the cluster numbering (the output is then returned and can be saved as a new LG_hom_stack object):
LGHomDf_P2_1 <- cluster_per_LG(LG = 3,
linkage_df = SN_SN_P2[SN_SN_P2$phase == "coupling",],
LG_hom_stack = LGHomDf_P2,
LOD_sequence = 3,
modify_LG_hom_stack = TRUE,
network.layout = "n",
nclust_out = 4)
## Total number of edges: 185
table(LGHomDf_P2_1$homologue, LGHomDf_P2_1$LG)
##
## 1 2 3 4 5
## 1 11 13 10 9 10
## 2 9 13 13 9 10
## 3 3 7 11 13 6
## 4 5 11 11 11 9
8.3 Cross-ploidy populations (e.g. tetraploid x diploid to give triploid F1)
polymapR is also able to map triploid populations, and may be extended to other cross-ploidy populations if there is demand for that in future. However, because of differences in parental ploidy levels, we have implemented a mapping approach that assumes that the parent of higher ploidy level is parent 1. Usually in polymapR a cross is understood to be between “maternal parent” x “paternal parent”. However this is just a convenient convention, and parent numbering can be reversed if so desired. In the case of triploid populations, there is no flexibility anymore, parent 1 must be tetraploid and parent 2 must be diploid.
The next point to be aware of is that the segregation of markers in the diploid parent is assumed to follow disomic inheritance (there is no other option in a diploid), whereas the segregation of the tetraploid parent is assumed to follow tetrasomic inheritance. Mixtures with preferential pairing in the tetraploid parent are currently not implemented for triploid populations.
Because of this disomic behaviour in the diploid parent, we are unable to use the normal clustering function (cluster_SN_markers) for identifying the homologues. Proceed as normal with the steps described above for parent 1. However in order to proceed correctly for parent 2 (diploid), you will need to provide the markertype as c(0,1) and not c(1,0). This is easy to understand since a 1x0 marker in a triploid population is Aaaa x aa while a 0x1 marker is aaaa x Aa. These are clearly different marker types and we cannot reverse them. To access the correct linkage functions for the diploid parent, therefore, you will need to think in terms of 1x0 markers segregating from parent 1, and 0x1 markers segregating from parent 2.
If you are confused about what marker types are possible in a triploid population, run the following, which gives all possible marker combinations (and the segregation ratios in the F1):
get("seg_p3_random",envir=getNamespace("polymapR"))
## dosage1 dosage2 f1_0 f1_1 f1_2 f1_3
## 1 0 0 1 0 0 0
## 2 0 1 1 1 0 0
## 3 0 2 0 1 0 0
## 4 1 0 1 1 0 0
## 5 1 1 1 2 1 0
## 6 1 2 0 1 1 0
## 7 2 0 1 4 1 0
## 8 2 1 1 5 5 1
## 9 2 2 0 1 4 1
## 10 3 0 0 1 1 0
## 11 3 1 0 1 2 1
## 12 3 2 0 0 1 1
## 13 4 0 0 0 1 0
## 14 4 1 0 0 1 1
## 15 4 2 0 0 0 1
To run the preliminary linkage analysis of 0x1 markers in parent 2, we also need to specify the ploidy of parent 2, in this case ploidy2 = 2. Note that the target parent is P2 as this is the parent we are now interested in, but markertype1 = c(0,1). For this we need a special dataset with triploid data, a sample of which is available from the polymapR package also:
data("TRI_dosages")
SN_SN_P2.tri <- linkage(dosage_matrix = TRI_dosages,
markertype1 = c(0,1),
target_parent = "P2", #Note the target parent is P2
other_parent = "P1",
ploidy = 4,
ploidy2 = 2,
pairing = "random"
)
Have a look at the plot of r versus LOD score - here you should notice that both the coupling and repulsion phases are equally informative now and cannot therefore be separated by LOD score as before:
plot_linkage_df(SN_SN_P2.tri)
If we try to cluster using the old function, we get the following result - no separation of the homologues, but the chromosomal linkage groups are identified:
P2_homologues.tri <- cluster_SN_markers(linkage_df = SN_SN_P2.tri,
LOD_sequence = seq(3, 10, 1),
LG_number = 3,
ploidy = 2, #because P2 is diploid..
parentname = "P2",
plot_network = FALSE,
plot_clust_size = FALSE)
## Total number of edges: 345
To separate the 0x1 markers into homologues, we use the phase information from the linkage analysis directly, using the phase_SN_diploid function. To learn more about how this function works, input ?phase_SN_diploid first, to get an idea of the arguments used.
LGHomDf_P2.tri <- phase_SN_diploid(linkage_df = SN_SN_P2.tri,
cluster_list = P2_homologues.tri,
LOD_chm = 4, #LOD at which chromosomes are identified
LG_number = 3) #number of linkage groups
## Total number of edges: 161
## Complete phase assignment possible using only coupling information at LOD 4
The rest of the analysis should be pretty much the same as for tetraploids or hexaploids, although you should check each time you use a new function whether that function has a ploidy2 argument, and if it does, use it!
9. Assigning SxS and DxN markers and consensus linkage group (LG) names
Assigning simplex x simplex markers to a homologue and linkage group is done by calculating linkage between SxS and SxN markers and after that assigning them based on this linkage using assign_linkage_group:
SN_SS_P1 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
markertype2 = c(1,1),
target_parent = "P1",
other_parent = "P2",
ploidy = 4,
pairing = "random")
P1_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P1,
LG_hom_stack = LGHomDf_P1,
SN_colname = "marker_a",
unassigned_marker_name = "marker_b",
phase_considered = "coupling",
LG_number = 5,
LOD_threshold = 3,
ploidy = 4)
##
## #### Marker(s) showing ambiguous linkage to more than one LG:
##
## |_ |
## |:-----------|
## |Ap_ts069042 |
##
## In total, 296 out of 301 markers were assigned.
##
## #### Marker(s) not assigned:
##
## |_ |_ |_ |_ |
## |:-----------|:-----------|:-----------|:-----------|
## |Ac_ws033297 |Ap_ts089267 |St_ns048176 |St_ts030296 |
## |Zm_ws010615 | | | |
## Assigned_LG LG1 LG2 LG3 LG4 LG5 Hom1 Hom2 Hom3 Hom4 Assigned_hom1
## Ac_ns002135 4 0 0 0 7 0 0 7 0 0 2
## Ac_ns002510 2 0 7 0 0 0 0 0 0 7 4
## Ac_ns024650 2 0 11 0 0 0 11 0 0 0 1
## Ac_ns028519 5 0 0 0 0 5 0 5 0 0 2
## Ac_ns028533 5 0 0 0 0 3 0 0 3 0 3
## Ac_ns029178 3 0 0 12 0 0 0 12 0 0 2
## Assigned_hom2 Assigned_hom3 Assigned_hom4
## Ac_ns002135 NA NA NA
## Ac_ns002510 NA NA NA
## Ac_ns024650 NA NA NA
## Ac_ns028519 NA NA NA
## Ac_ns028533 NA NA NA
## Ac_ns029178 NA NA NA
SN_SS_P2 <- linkage(dosage_matrix = filtered_data,
markertype1 = c(1,0),
markertype2 = c(1,1),
target_parent = "P2",
other_parent = "P1",
ploidy = 4,
pairing = "random")
P2_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P2,
LG_hom_stack = LGHomDf_P2_1,
SN_colname = "marker_a",
unassigned_marker_name = "marker_b",
phase_considered = "coupling",
LG_number = 5,
LOD_threshold = 3,
ploidy = 4)
##
## In total, 300 out of 300 markers were assigned.
As simplex x simplex markers are present in both parents, we can define which linkage groups correspond with each other between parents. After this, one of the linkage groups of the parents should be renamed and the SxS markers assigned again according to the new linkage group names:
LGHomDf_P2_2 <- consensus_LG_names(modify_LG = LGHomDf_P2_1,
template_SxS = P1_SxS_Assigned,
modify_SxS = P2_SxS_Assigned)
##
## #### Original LG names
##
## | | 1| 2| 3| 4| 5|
## |:--|--:|--:|--:|--:|--:|
## |1 | 0| 0| 0| 47| 0|
## |2 | 64| 0| 0| 0| 0|
## |3 | 0| 0| 58| 0| 0|
## |4 | 0| 58| 0| 0| 0|
## |5 | 0| 0| 0| 0| 68|
##
## #### Modified LG names
##
## | | 1| 2| 3| 4| 5|
## |:--|--:|--:|--:|--:|--:|
## |1 | 47| 0| 0| 0| 0|
## |2 | 0| 64| 0| 0| 0|
## |3 | 0| 0| 58| 0| 0|
## |4 | 0| 0| 0| 58| 0|
## |5 | 0| 0| 0| 0| 68|
P2_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P2,
LG_hom_stack = LGHomDf_P2_2,
SN_colname = "marker_a",
unassigned_marker_name = "marker_b",
phase_considered = "coupling",
LG_number = 5,
LOD_threshold = 3,
ploidy = 4)
##
## In total, 300 out of 300 markers were assigned.
Since we now have a consistent linkage group numbering, we can also assign the DxN markers:
P1_DxN_Assigned <- assign_linkage_group(linkage_df = SN_DN_P1,
LG_hom_stack = LGHomDf_P1,
SN_colname = "marker_a",
unassigned_marker_name = "marker_b",
phase_considered = "coupling",
LG_number = 5,
LOD_threshold = 3,
ploidy = 4)
##
## In total, 111 out of 111 markers were assigned.
P2_DxN_Assigned <- assign_linkage_group(linkage_df = SN_DN_P2,
LG_hom_stack = LGHomDf_P2_2,
SN_colname = "marker_a",
unassigned_marker_name = "marker_b",
phase_considered = "coupling",
LG_number = 5,
LOD_threshold = 3,
ploidy = 4)
##
## In total, 101 out of 101 markers were assigned.
15. Preferential pairing
Up to now, we have been working under the assumption of random bivalent pairing in the parents. However, it has been shown that in certain species the pairing behaviour is neither completely random nor completely preferential (as we have in allopolyploids) but something in-between, a condition termed segmental allopolyploidy [5]. A comprehensive study on this topic was recently performed in rose [6], which found evidence of disomic behaviour in certain regions of the genome, with tetrasomic behaviour everywhere else. This sort of mixed inheritance can complicate the analysis and in cases where this effect is particularly pronounced, it is probably unwise to ignore [7].
The polymapR package can currently accomodate preferential pairing in tetraploid species only, by correcting the recombination frequency estimates once the level of preferential pairing in a parental chromosome is known. In order to determine whether this is necessary, we can run some diagnostic tests using the marker data after an initial map has been produced which can inform our choice regarding whether to include a preferential pairing parameter or not (in a subsequent re-analysis).
First, a word of warning. There is a function to test for preferential pairing among pairs of closely-linked simplex x nulliplex markers in repulsion phase within polymapR. However, it is certainly also advised to use multi-point methods to test for preferential pairing (if possible). The only possible drawback is that multi-point methods (using all markers across a chromosome) generally assume uniform pairing behaviour across the length of a chromosome, which may not always occur [6]). A robust multi-point approach as described in [6] used identity-by-descent probabilities for the population, estimated using TetraOrigin [8], which estimates the most likely pairing behaviour per individual and can therefore reveal whether there are deviations from the assumption of random chromosomal pairing and recombination across the population.
The function test_prefpairing examines closely-linked repulsion-phase simplex x nulliplex marker pairs only, and is therefore the results of this function may not be as robust as a multi-point whole-chromosome approach. Be that as it may, it is probably the simplest check we can perform post-mapping for the possibility of unusual pairing behaviour. The method used here is described also in Bourke et al.[6], which is a development on ideas originally described by Qu and Hancock[9]. The idea is to look for signatures of preferential pairing in repulsion-phase pairs which map to the same locus.
We first need to define a minimum distance by which to consider markers essentially at the same locus. This is somewhat arbitrary - in our study of tetraploid rose, we used pairs of duplicated individuals to determine an approximate error rate which, along with considerations of the size of the population and rate of missing values, led us to conclude that distances less than ~1 cM were close enough to be informative (and for which we could assume 0 recombinations). The smaller the distance we choose here the better, although by choosing a very small distance we are limiting the number of repulsion-pairs available. By default in test_prefpairing we use a value of 0.5 cM. This can be increased or decreased as required using the min_cM argument.
P1.prefPairing <- test_prefpairing(dosage_matrix = ALL_dosages,
maplist = integrated.maplist,
LG_hom_stack = LGHomDf_P1_1,
min_cM = 1, #changed from default of 0.5 cM
ploidy = 4)
head(P1.prefPairing)
The output shows the repulsion-phase marker pairs tested. The most important column is probably the column P_value.adj, which gives the FDR-adjusted P-values from a Binomial test of the hypothesis that the disomic repulsion estimate of recombination frequency differs significantly from \(\frac{1}{3}\) (actually it is a one-sided test, so we are testing that r disom falls below \(\frac{1}{3}\)). In cases where this is significant (with P_value.adj < 0.01, say), we might have some evidence for preferential pairing (but from that marker pair only!). The output of test_prefpairing also lists the markers, their positions and on which homologues they have been mapped. There is also an estimate of the pairwise preferential pairing parameter pref.p. Note that there are two possible (maximum-likelihood) estimators of a preferential pairing parameter in the case of a pair of simplex x nulliplex markers, depending on whether the marker alleles reside on homologous chromosomes, or homoeologous chromosomes (i.e. whether the pair are contained within a “subgenome”, or are straddled across “subgenomes”, to borrow the language of allopolyploids for convenience). If there is strong evidence of preferential pairing from multiple marker pairs, the approach we recommend is to take the average of the pref.p values for a particular chromosomal linkage group, which becomes one of the parameters p1 or p2 in the linkage function.
Here, we have purely random pairing in our parents, so the estimates of pref.p should not be considered meaningful. The preferential parameter p is defined such that \(0 < p < \frac{2}{3}\). Supposing this were not the case, and we wanted to estimate the strength of preferential pairing on LG1 of parent 1, we would do something like the following:
mean(P1.prefPairing[P1.prefPairing$P_value.adj < 0.01 & P1.prefPairing$LG_a == 1,]$pref.p)
Note that again we require first that there be significant evidence of disomic behaviour before including the estimate. If there is evidence of preferential pairing, it is probably a good idea to go back and reduce the min_cM argument even further, since we really want this value to be as close to 0 as possible for an accurate estimate of a preferential pairing parameter.
Finally, to re-run the linkage analysis with this estimate for parent 1 of linkage group 1, we would first have to subset our marker dosage data corresponding to that linkage group only (since we do not need to re-analyse any other chromosome!), and then re-run the linkage analysis. Note that in cases with extreme levels of preferential pairing, the clustering into separate homologues can also be affected. Supposing we had estimated a value of 0.25 for the preferential pairing parameter we could proceed as follows:
lg1_markers <- unique(c(rownames(marker_assignments_P1[marker_assignments_P1[,"Assigned_LG"] == 1,]),
rownames(marker_assignments_P2[marker_assignments_P2[,"Assigned_LG"] == 1,])))
all_linkages_list_P1_lg1 <- finish_linkage_analysis(marker_assignment = marker_assignments$P1[lg1_markers,],
dosage_matrix = filtered_data[lg1_markers,],
target_parent = "P1",
other_parent = "P2",
convert_palindrome_markers = FALSE,
ploidy = 4,
pairing = "preferential",
prefPars = c(0.25,0), #just for example!
LG_number = 1 #interested in just 1 chm.
)
