The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
The following example is intended to show the usage of
OneMap
functions for linkage mapping in outcrossing
(non-inbred) populations. With basic knowledge of R syntax, one should
have no problems using it. If you are not familiar with R, we recommend
reading the vignette Introduction
to R.
Hopefully, these examples will be clear enough to help any user to understand its functionality and start using it. You do not need to be an expert in R to build your linkage map, but some concepts are necessary and will help you through the process.
There is a GitHub OneMap
version which is continuously
improved, we strongly recommend all users to try this version. On
augusto-garcia/onemap
GitHub page you can find instructions
to install the package from GitHub and also more fancy tutorials.
This step may be quite difficult because the data file is not very simple, and some errors can occur while reading it. The input file format is similar to that used by MAPMAKER/EXP (Lander et al., 1987), so experienced users of genetic analysis software should be already familiar with this scenario.
The input file is a text file, where the first line indicates the
cross-type, and the second line provides information about the number of
individuals, the number of markers, the presence of physical marker
locations, and the presence of phenotypic data. The third line contains
the sample IDs. Then, the genotype information is included separately
for each marker. The character *
indicates the beginning of
information input for a new marker, followed by the marker name. Next,
there is a code indicating the marker type, according to Wu’s et
al. (2002a) notation. It is recommended to check Wu’s et al. (2002a)
paper before using OneMap
.
Marker types must be one of the following: A.1
,
A.2
, A.3
, A.4
, B1.5
,
B2.6
, B3.7
, C.8
,
D1.9
, D1.10
, D1.11
,
D1.12
, D1.13
, D2.14
,
D2.15
, D2.16
, D2.17
or
D2.18
, each one corresponding to a row of the following
table:
Parent | Offspring | |||||
---|---|---|---|---|---|---|
Crosstype | Cross | Observed bands | Observed bands | Segregation | ||
\(A\) | 1 | \(ab \times cd\) | \(ab \times cd\) | \(ac, ad, bc, bd\) | \(1:1:1:1\) | |
2 | \(ab \times ac\) | \(ab \times ac\) | \(a, ac, ba, bc\) | \(1:1:1:1\) | ||
3 | \(ab \times co\) | \(ab \times c\) | \(ac, a, bc,b\) | \(1:1:1:1\) | ||
4 | \(ao \times bo\) | \(a \times b\) | \(ab, a, b, o\) | \(1:1:1:1\) | ||
\(B\) | \(B_1\) | 5 | \(ab \times ao\) | \(ab \times a\) | \(ab, 2a,b\) | \(1:2:1\) |
\(B_2\) | 6 | \(ao \times ab\) | \(a \times ab\) | \(ab,2a,b\) | \(1:2:1\) | |
\(B_3\) | 7 | \(ab \times ab\) | \(ab \times ab\) | \(a, 2ab,b\) | \(1:2:1\) | |
\(C\) | 8 | \(ao \times ao\) | \(a \times a\) | \(3a, o\) | \(3:1\) | |
\(D\) | \(D_1\) | 9 | \(ab \times cc\) | \(ab \times c\) | \(ac, bc\) | \(1:1\) |
10 | \(ab \times aa\) | \(ab \times a\) | \(a, ab\) | \(1:1\) | ||
11 | \(ab \times oo\) | \(ab \times o\) | \(a, b\) | \(1:1\) | ||
12 | \(bo \times aa\) | \(b \times a\) | \(ab, a\) | \(1:1\) | ||
13 | \(ao \times oo\) | \(a \times o\) | \(a, o\) | \(1:1\) | ||
\(D_2\) | 14 | \(cc \times ab\) | \(c \times ab\) | \(ac, bc\) | \(1:1\) | |
15 | \(aa \times ab\) | \(a \times ab\) | \(a, ab\) | \(1:1\) | ||
16 | \(oo \times ab\) | \(o \times ab\) | \(a, b\) | \(1:1\) | ||
17 | \(aa \times bo\) | \(a \times b\) | \(ab, a\) | \(1:1\) | ||
18 | \(oo \times ao\) | \(o \times a\) | \(a, o\) | \(1:1\) |
Letters A
, B
, C
and
D
indicate the segregation type (i.e.,
1:1:1:1
, 1:2:1
, 3:1
or
1:1
, respectively), while the number after the dot
(e.g., A.1
) indicates the observed bands in the
offspring. The paper cited above gives details with respect to marker
types; we will not discuss them here, but it is easy to see that each
marker is classified based on the band patterns of parents and
progeny.
Finally, after each marker name, comes the genotype data for the
segregating population. The coding for marker genotypes used by
OneMap
is also the same one proposed by Wu et al. (2002a),
and the possible values vary according to the specific marker type.
Missing data are indicated with the character -
(minus
sign), and an empty space separates the information for each individual.
Phenotype information, if present, follows genotypic data with a similar
structure. Details are found with the help of function
read_onemap
.
Here is an example of such a file for 10 individuals and 5 markers (the three zeros in the second line indicate that there is no chromosome information, physical position information, or phenotypic data, respectively). It is very similar to a MAPMAKER/EXP file, but has additional information about the cross-type.
data type outcross
10 5 0 0 0
I1 I2 I3 I4 I5 I6 I7 I8 I9 I10
*M1 B3.7 ab ab - ab b ab ab - ab b
*M2 D2.18 o - a a - o a - o o
*M3 D1.13 o a a o o - a o a o
*M4 A.4 ab b - ab a b ab b - a
*M5 D2.18 a a o - o o a o o o
In case you have physical chromosome and position information:
data type outcross
10 5 1 1 0
I1 I2 I3 I4 I5 I6 I7 I8 I9 I10
*M1 B3.7 ab ab - ab b ab ab - ab b
*M2 D2.18 o - a a - o a - o o
*M3 D1.13 o a a o o - a o a o
*M4 A.4 ab b - ab a b ab b - a
*M5 D2.18 a a o - o o a o o o
*CHROM 1 1 1 2 2
*POS 2391 3812 5281 1823 3848
Notice that once the marker type is identified, no variations of
symbols presented on the table for the observed bands
are allowed. For example, for A.1
, only ac
,
ad
, bc
, and bd
genotypes are
expected (plus missing values). We notice in FAQs that this is a
common mistake made by users, so please be careful.
The input file must be saved in text format, with extensions like
.raw
. It is a good idea to open the text file called
onemap_example_out.raw
(available with OneMap
and saved in the directory you installed it) to see how this file should
be. You can see where OneMap
is installed using the
command:
system.file(package = "onemap")
In the session Importing data from VCF file
below, you will see how to import VCF files as OneMap
objects.
Once the input file is created, the data can be loaded and saved into
an R onemap
object. The function used to import data is
named read_onemap
. Its usage is quite simple:
<- read_onemap(dir = "C:/workingdirectory", inputfile = "onemap_example_out.raw") onemap_example_out
The first argument is the directory where the input file is located,
so modify it accordingly. The second one is the data file name. In this
example, an object named onemap_example_out
was created. If
you leave the argument dir
blank, the file will be loaded
from your working directory
.
You can change the working directory in R using function
setwd()
or in the toolbar clicking
File -> Change dir
. If you set your working directory to
the one containing the input file, you can just type:
<- read_onemap(inputfile = "onemap_example_out.raw") onemap_example_out
If no error has occurred, a message will display some basic information about the data, such as number of individuals and number of markers:
#> Working...
#>
#> --Read the following data:
#> Type of cross: outcross
#> Number of individuals: 100
#> Number of markers: 30
#> Chromosome information: no
#> Position information: no
#> Number of traits: 3
#> Missing trait values:
#> Pheno1: 0
#> Pheno2: 3
#> Pheno3: 0
Because this particular dataset is distributed along with the package, as an alternative you can load it by typing:
data("onemap_example_out")
Loading the data creates an object of class onemap
,
which will further be used in the analysis. R command print
recognizes objects of this class. Thus, if you type:
onemap_example_out
you will see some information about the object:
#> This is an object of class 'onemap'
#> Type of cross: outcross
#> No. individuals: 100
#> No. markers: 30
#> CHROM information: no
#> POS information: no
#> Percent genotyped: 100
#>
#> Segregation types:
#> A.1 --> 3
#> A.2 --> 1
#> A.4 --> 4
#> B1.5 --> 1
#> B2.6 --> 2
#> B3.7 --> 5
#> C.8 --> 2
#> D1.10 --> 2
#> D1.12 --> 1
#> D1.13 --> 2
#> D2.15 --> 1
#> D2.16 --> 2
#> D2.17 --> 2
#> D2.18 --> 2
#>
#> No. traits: 3
#> Missing trait values:
#> Pheno1: 0
#> Pheno2: 3
#> Pheno3: 0
Also, you can use the plot.onemap
function to see
graphically markers genotypes:
plot(onemap_example_out)
Changing the argument all
to FALSE
, the
markers will be separated by their type. In this case, you can note that
the graphic cell size will adapt to the number of markers of the same
type. In other words, the higher is the number of markers with the same
type, the lower is the cell for this type.
plot(onemap_example_out, all = FALSE)
This function can take quite some time, depending on the number of
markers involved. More information about this plot function can be found
using ?plot.onemap
.
Also, you can see the number of markers by segregation pattern with
the plot_by_segreg_type
function:
You can import information from VCF
to
OneMap
using onemap_read_vcfR
function.
With the onemap_read_vcfR
you can convert the object
from vcfR
package directly to onemap
. The
onemap_read_vcfR
function keeps chromosome and position
information for each marker in the onemap object generated.
We will use the example file vcf_example_out.vcf.gz
to
show how it works, which contains markers from the same population of
onemap_example_out.raw
.
Here we use the the vcfR
package internally to help this
conversion. The vcfR
authors mentioned in their tutorials
that RAM memory use is an important consideration when using the
package. Depending of your dataset, the object created can be huge and
occupy a lot of memory.
You can use onemap_read_vcfR
function to convert the VCF
file to onemap
object. The parameters used are the
vcf
with the VCF file path, the identification of each
parent (here, you must define only one sample for each parent) and the
cross type.
<- onemap_read_vcfR(vcf = system.file("extdata/vcf_example_out.vcf.gz", package = "onemap"),
vcf_example_out parent1 = "P1",
parent2 = "P2",
cross = "outcross")
Depending on your dataset, this function can take some time to run.
Note that the conversion filter out markers which are not informative
for the informed cross type. For example, in outcrossing species,
markers that have both parents homozygous (aa x bb) do not inform
recombinations and are removed of the data set. Only markers types
contained the table at Creating the
data file are kept in the onemap object. Function
onemap_read_vcfR
print at the screen the reason why markers
were filtered.
You can also have more missing data in the returned object compared
with the VCF because the onemap_read_vcfR
replace by
missing data the genotypes that are not expected for that marker type.
For example, for a marker type D1.10 (ab
x
aa
), we only expect aa
and ab
genotypes, if there are bb
genotypes they will be replaced
by missing data. You can see the percentage of missing data at the
resulted onemap object with:
vcf_example_out#> This is an object of class 'onemap'
#> Type of cross: outcross
#> No. individuals: 92
#> No. markers: 24
#> CHROM information: yes
#> POS information: yes
#> Percent genotyped: 99
#>
#> Segregation types:
#> B3.7 --> 18
#> D1.10 --> 6
#>
#> No. traits: 0
NOTE:From version 2.0.6 to 2.1.1005,
OneMap
had the vcf2raw
function to convert
vcf
to .raw
. Now, this function is defunct,
but it can be replaced by a combination of onemap_read_vcfR
and write_onemap_raw
functions. See Exporting .raw file from
onemap object session to further information about
write_onemap_raw
.
If your onemap object has too many missing genotypes you can face problems during the analysis. Check the percentage of missing genotypes in our data set printing the onemap object:
vcf_example_out#> This is an object of class 'onemap'
#> Type of cross: outcross
#> No. individuals: 92
#> No. markers: 24
#> CHROM information: yes
#> POS information: yes
#> Percent genotyped: 99
#>
#> Segregation types:
#> B3.7 --> 18
#> D1.10 --> 6
#>
#> No. traits: 0
Our example has 1% of missing genotypes (99% are genotyped). If you
want to filter markers according to their percentage of missing data,
you can use the function filter_missing
:
<- filter_missing(vcf_example_out, threshold = 0.25)
vcf_filtered #> Number of markers removed from the onemap object: 0
Any of our markers were filtered, because, in this example, we do not have much missing data.
Function create_depth_profile
generates dispersion
graphics with x and y-axis representing, respectively, the reference and
alternative allele depths. The function is only available for biallelic
markers in VCF files with allele counts information. Each dot represents
a genotype for mks
markers and inds
individuals. If both arguments receive NULL
, all markers
and individuals are considered. Dots are colored according to the
genotypes present in the onemap object (GTfrom = onemap
) or
VCF file (GTfrom = vcf
). A rds file is generated with the
data in the graphic (rds.file
). The alpha
argument controls the transparency of the color of each dot. Control
this parameter is a good idea when having a large number of markers and
individuals. The x_lim
and y_lim
control the
axis scale limits; by default, it uses the maximum value of the
counts.
Here is an example for the vcf_example_out
dataset.
# For outcrossing population
create_depths_profile(onemap.obj = vcf_example_out,
vcf = system.file("extdata/vcf_example_out.vcf.gz", package = "onemap"),
parent1 = "P1",
parent2 = "P2",
vcf.par = "AD",
recovering = FALSE,
mks = NULL,
inds = NULL,
GTfrom = "vcf",
alpha = 0.1,
rds.file = "depths_out.rds")
Because the genotypes are from VCF file, the legend points the VCF
codification? ./.
represent missing data; 0/0
homozygotes for reference alleles; 0/1
heterozygotes;
1/1
homozygotes for alternative alleles. You can also have
phased genotypes represented which have pipe |
instead of
bar /
.
OneMap
datasetsIf you have more than one dataset of markers, all from the same
mapping population, you can use the function combine_onemap
to merge them into only one onemap
object.
In our example, we have two datasets:
onemap_example_out
with 30 markers and 100
individualsvcf_example_out
with 24 biallelic markers and
92 individuals.The combine_function
recognizes the correspondent
individuals by the ID, thus, it is important to define the same IDs to
respective individuals in both raw
files. Compared with the
first file, the second file does not have markers information for 8
individuals. The combine_onemap
will complete this
information with NA.
In our examples, we have only genotypic information, but the function can also merge the phenotypic information if it exists.
<- combine_onemap(onemap_example_out, vcf_example_out)
comb_example
comb_example#> This is an object of class 'onemap'
#> Type of cross: outcross
#> No. individuals: 100
#> No. markers: 54
#> CHROM information: yes
#> POS information: yes
#> Percent genotyped: 96
#>
#> Segregation types:
#> A.1 --> 3
#> A.2 --> 1
#> A.4 --> 4
#> B1.5 --> 1
#> B2.6 --> 2
#> B3.7 --> 23
#> C.8 --> 2
#> D1.10 --> 8
#> D1.12 --> 1
#> D1.13 --> 2
#> D2.15 --> 1
#> D2.16 --> 2
#> D2.17 --> 2
#> D2.18 --> 2
#>
#> No. traits: 3
#> Missing trait values:
#> Pheno1: 0
#> Pheno2: 3
#> Pheno3: 0
The function arguments are the names of the onemap
objects you want to combine.
Plotting markers genotypes from the outputted onemap
object, we can see that there are more missing data -
(black vertical lines) for some individuals because they were missing in
the second file.
plot(comb_example)
It is possible that there are redundant markers in your dataset, especially when dealing with too many markers. Redundant markers have the same genotypic information that other markers because they didn’t happen recombination events between each other. They will not increase information on the map but will increase computational effort during the map building. Therefore, it is a good practice to remove them to build the map and, once the map is already built, they can be added again.
First, we use the function find_bins
to group the
markers into bins according to their genotypic information. In other
words, markers with the same genotypic information will be in the same
bin.
<- find_bins(comb_example, exact = FALSE)
bins
bins#> This is an object of class 'onemap_bin'
#> No. individuals: 100
#> No. markers in original dataset: 54
#> No. of bins found: 52
#> Average of markers per bin: 1.038462
#> Type of search performed: non exact
The first argument is the onemap
object and the
exact
argument specifies if only markers with the same
information will be at the same bin. Using FALSE
at this
second argument, missing data will not be considered, and the marker
with the lowest amount of missing data will be the representative marker
on the bin.
Our example dataset has only two redundant markers. We can create a
new onemap
object without them, using the
create_data_bins
function. This function keeps only the
most representative marker of each bin from the bins
object.
<- create_data_bins(comb_example, bins)
bins_example
bins_example#> This is an object of class 'onemap'
#> Type of cross: outcross
#> No. individuals: 100
#> No. markers: 52
#> CHROM information: yes
#> POS information: yes
#> Percent genotyped: 96
#>
#> Segregation types:
#> A.1 --> 3
#> A.2 --> 1
#> A.4 --> 4
#> B1.5 --> 1
#> B2.6 --> 2
#> B3.7 --> 22
#> C.8 --> 2
#> D1.10 --> 7
#> D1.12 --> 1
#> D1.13 --> 2
#> D2.15 --> 1
#> D2.16 --> 2
#> D2.17 --> 2
#> D2.18 --> 2
#>
#> No. traits: 3
#> Missing trait values:
#> Pheno1: 0
#> Pheno2: 3
#> Pheno3: 0
The arguments for the create_data_bins
function are the
onemap
object and the object created by the
find_bins
function.
The functions onemap_read_vcfR
generates new onemap
objects without use a input .raw
file. Also, the functions
combine_onemap
and create_data_bins
manipulate
the information of the original .raw
file and creates a new
dataset. In both cases, you do not have an input file .raw
that contains the same information as the analyzed data. If you want to
create a new input file with the dataset you are working on after using
these functions, you can use the function
write_onemap_raw
.
write_onemap_raw(bins_example, file.name = "new_dataset.raw")
The file new_dataset.raw
will be generated in your
working directory. In our example, it contains only non-redundant
markers from onemap_example_out
and
vcf_example_out
datasets.
For the map building process, it is also important to know which markers have deviations in the expected segregation pattern. It can be a good practice to remove them from the map building process, because they can adversely affect the map building, and, once the map is built, they can be inserted.
The function test_segregation_of_a_marker
performs a
chi-square test according to Mendelian segregation to check if a
specific marker is following the expected segregation pattern.
test_segregation_of_a_marker(bins_example, 4)
#> $Hypothesis
#> [1] "1:1:1:1"
#>
#> $qui.quad
#> X-squared
#> 2.64
#>
#> $p.val
#> [1] 0.4505201
#>
#> $perc.genot
#> [1] 100
The arguments are the onemap
object and the number of
the marker you want to test.
You can also test all the markers in your onemap
object
using the test_segregation
function. The results can be
viewed by printing the output object of class
onemap_segreg_test
.
<- test_segregation(bins_example)
segreg_test print(segreg_test)
#> Marker H0 Chi-square p-value % genot.
#> 1 M1 1:2:1 1.76000000 4.147829e-01 100
#> 2 M2 1:1 0.04000000 8.414806e-01 100
#> 3 M3 1:1 0.36000000 5.485062e-01 100
#> 4 M4 1:1:1:1 2.64000000 4.505201e-01 100
#> 5 M5 1:1 1.96000000 1.615133e-01 100
#> 6 M6 1:2:1 1.52000000 4.676664e-01 100
#> 7 M7 1:1 0.16000000 6.891565e-01 100
#> 8 M8 1:2:1 0.86000000 6.505091e-01 100
#> 9 M9 1:1 0.04000000 8.414806e-01 100
#> 10 M10 1:1 0.36000000 5.485062e-01 100
#> 11 M11 1:1 0.16000000 6.891565e-01 100
#> 12 M12 1:1:1:1 6.48000000 9.045460e-02 100
#> 13 M13 3:1 0.00000000 1.000000e+00 100
#> 14 M14 1:1:1:1 0.40000000 9.402425e-01 100
#> 15 M15 1:1:1:1 2.24000000 5.241127e-01 100
#> 16 M16 1:1 1.44000000 2.301393e-01 100
#> 17 M17 1:2:1 35.58000000 1.878889e-08 100
#> 18 M18 1:1:1:1 1.44000000 6.961859e-01 100
#> 19 M19 1:2:1 37.98000000 5.659105e-09 100
#> 20 M20 1:1:1:1 4.88000000 1.807980e-01 100
#> 21 M21 1:1 1.44000000 2.301393e-01 100
#> 22 M22 1:1 1.00000000 3.173105e-01 100
#> 23 M23 3:1 0.48000000 4.884223e-01 100
#> 24 M24 1:2:1 1.50000000 4.723666e-01 100
#> 25 M25 1:2:1 35.04000000 2.461278e-08 100
#> 26 M26 1:1:1:1 1.52000000 6.776621e-01 100
#> 27 M27 1:1 1.00000000 3.173105e-01 100
#> 28 M28 1:1:1:1 1.20000000 7.530043e-01 100
#> 29 M29 1:1 0.00000000 1.000000e+00 100
#> 30 M30 1:2:1 3.42000000 1.808658e-01 100
#> 31 SNP1 1:2:1 3.73913043 1.541907e-01 92
#> 32 SNP2 1:2:1 4.76086957 9.251035e-02 92
#> 33 SNP3 1:2:1 4.76086957 9.251035e-02 92
#> 34 SNP5 1:2:1 4.95652174 8.388899e-02 92
#> 35 SNP6 1:2:1 1.80434783 4.056868e-01 92
#> 36 SNP7 1:2:1 1.45652174 4.827478e-01 92
#> 37 SNP8 1:2:1 0.26086957 8.777137e-01 92
#> 38 SNP9 1:2:1 1.39130435 4.987491e-01 92
#> 39 SNP10 1:2:1 0.06521739 9.679172e-01 92
#> 40 SNP11 1:2:1 0.52173913 7.703814e-01 92
#> 41 SNP12 1:2:1 0.17391304 9.167170e-01 92
#> 42 SNP13 1:2:1 0.08695652 9.574534e-01 92
#> 43 SNP14 1:1 0.53846154 4.630710e-01 91
#> 44 SNP16 1:2:1 0.00000000 1.000000e+00 92
#> 45 SNP17 1:1 0.00000000 1.000000e+00 90
#> 46 SNP18 1:1 0.27472527 6.001795e-01 91
#> 47 SNP20 1:1 0.17777778 6.732900e-01 90
#> 48 SNP21 1:1 0.01123596 9.155825e-01 89
#> 49 SNP22 1:2:1 1.10869565 5.744468e-01 92
#> 50 SNP23 1:2:1 2.15217391 3.409270e-01 92
#> 51 SNP24 1:2:1 2.86956522 2.381671e-01 92
#> 52 SNP25 1:2:1 2.15217391 3.409270e-01 92
The only argument of the function is a onemap
object.
Once we have the onemap_segreg_test
object, the function
select_segreg
can be used to show only the markers
considered with/without segregation distortion. By default, it uses as a
threshold for the test a global \(\alpha=0.05\), corrected for multiple tests
with Bonferroni correction.
select_segreg(segreg_test, distorted = TRUE) #to show the markers names with segregation distortion
#> [1] "M17" "M19" "M25"
select_segreg(segreg_test, distorted = FALSE) #to show the markers names without segregation distortion
#> [1] "M1" "M2" "M3" "M4" "M5" "M6" "M7" "M8" "M9"
#> [10] "M10" "M11" "M12" "M13" "M14" "M15" "M16" "M18" "M20"
#> [19] "M21" "M22" "M23" "M24" "M26" "M27" "M28" "M29" "M30"
#> [28] "SNP1" "SNP2" "SNP3" "SNP5" "SNP6" "SNP7" "SNP8" "SNP9" "SNP10"
#> [37] "SNP11" "SNP12" "SNP13" "SNP14" "SNP16" "SNP17" "SNP18" "SNP20" "SNP21"
#> [46] "SNP22" "SNP23" "SNP24" "SNP25"
It is not recommended, but you can define a different threshold value
by changing the threshold
argument of the function
select_segreg
.
For the next steps, it will be useful to know the numbers of each marker with segregation distortion, so then you can keep those out of your map building analysis. These numbers refer to the lines where markers are located on the data file.
To access the corresponding number for these markers you can change
the numbers
argument:
<- select_segreg(segreg_test, distorted = TRUE, numbers = TRUE) #to show the markers numbers with segregation distortion
dist
dist#> [1] 17 19 25
<- select_segreg(segreg_test, distorted = FALSE, numbers = TRUE) #to show the markers numbers without segregation distortion
no_dist
no_dist#> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 20 21 22 23 24 26 27 28
#> [26] 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
You can also see the results graphically by:
plot(segreg_test)
Now, we start the map building analysis. In this example, we follow two different strategies:
Using only recombinations information.
Using the recombinations and also the reference genome
information, once our example has CHROM
and
POS
information for some of the markers.
First, we will apply the strategy using only recombinations information. In the second part of this tutorial, we show a way to use also reference genome information. At the end of our analysis, we will be able to compare these two strategies by drawing the resulted genetic maps.
The first step is estimating the recombination fraction between all pairs of markers, using two-point tests.
<- rf_2pts(bins_example) twopts
Although two-point tests were implemented in C language, which is usually much faster than R, this step can take quite some time, depending on the number of markers involved and their segregation type, because all combinations will be estimated and tested. Besides, the results use a a lot of memory and a rather powerful computer is needed.
When the two-point analysis is finished, an object of class
rf_2pts
is created. Typing
twopts
will show a message with the criteria used in the analysis and some other information:
#> This is an object of class 'rf_2pts'
#>
#> Criteria: LOD = 3 , Maximum recombination fraction = 0.5
#>
#> This object is too complex to print
#> Type 'print(object, c(mrk1=marker, mrk2=marker))' to see
#> the analysis for two markers
#> mrk1 and mrk2 can be the names or numbers of both markers
If you want to see the results for given markers, say M1
and M3
, the command is:
print(twopts, c("M1", "M3"))
#> Results of the 2-point analysis for markers: M1 and M3
#> Criteria: LOD = 3 , Maximum recombination fraction = 0.5
#>
#> rf LOD
#> CC 0.2954514 1.646878
#> CR 0.2954514 1.646878
#> RC 0.7045486 1.646878
#> RR 0.7045486 1.646878
Each line corresponds to a possible linkage phase. CC
denotes the coupling phase in both parents, CR
and
RC
denote coupling phase in parent 1 and 2, respectively,
and repulsion in the other, and RR
denotes the repulsion
phase in both parents. Value rf
is the maximum likelihood
estimate of the recombination fraction, with its corresponding LOD
Score.
Once the recombination fractions and linkage phases for all pairs of
markers have been estimated and tested, markers can be assigned to
linkage groups. To do this, first, use the function
make_seq
to create a sequence with the markers you want to
assign.
The function make_seq
is used to create sequences from
objects of several kinds, as will be seen along with this tutorial.
Here, the object is of class rf_2pts
and the second
argument specifies which markers one wants to use. If one wants to use
only a subset of markers, say M1
and M2
, the
option will be a vector with the corresponding numbers of the markers,
as c(1,2)
, you can also use a string "all"
to
specify that you want to analyze all markers. In our example, we will
use the vector with the numbers of the markers with no segregation
distortion.
<- make_seq(twopts, c(no_dist)) mark_no_dist
Because the identification of the markers can be cumbersome, one
should use the function marker type
to see their numbers,
names, and types:
marker_type(mark_no_dist)
#> Marker Marker.name Type
#> 1 1 M1 B3.7
#> 2 2 M2 D2.18
#> 3 3 M3 D1.13
#> 4 4 M4 A.4
#> 5 5 M5 D2.18
#> 6 6 M6 B3.7
#> 7 7 M7 D2.15
#> 8 8 M8 B3.7
#> 9 9 M9 D1.10
#> 10 10 M10 D2.17
#> 11 11 M11 D2.16
#> 12 12 M12 A.2
#> 13 13 M13 C.8
#> 14 14 M14 A.4
#> 15 15 M15 A.4
#> 16 16 M16 D2.17
#> 17 18 M18 A.1
#> 18 20 M20 A.1
#> 19 21 M21 D2.16
#> 20 22 M22 D1.10
#> 21 23 M23 C.8
#> 22 24 M24 B3.7
#> 23 26 M26 A.1
#> 24 27 M27 D1.12
#> 25 28 M28 A.4
#> 26 29 M29 D1.13
#> 27 30 M30 B3.7
#> 28 31 SNP1 B3.7
#> 29 32 SNP2 B3.7
#> 30 33 SNP3 B3.7
#> 31 34 SNP5 B3.7
#> 32 35 SNP6 B3.7
#> 33 36 SNP7 B3.7
#> 34 37 SNP8 B3.7
#> 35 38 SNP9 B3.7
#> 36 39 SNP10 B3.7
#> 37 40 SNP11 B3.7
#> 38 41 SNP12 B3.7
#> 39 42 SNP13 B3.7
#> 40 43 SNP14 D1.10
#> 41 44 SNP16 B3.7
#> 42 45 SNP17 D1.10
#> 43 46 SNP18 D1.10
#> 44 47 SNP20 D1.10
#> 45 48 SNP21 D1.10
#> 46 49 SNP22 B3.7
#> 47 50 SNP23 B3.7
#> 48 51 SNP24 B3.7
#> 49 52 SNP25 B3.7
OneMap has two different functions for grouping markers. The
group
function:
<- group(mark_no_dist)
LGs #> Selecting markers:
#> group 1
#> ........................
#> group 2
#> ................
#> group 3
#> ......
For this function, optional arguments are LOD
and
max.rf
, which define thresholds to be used when assigning
markers to linkage groups. If none is provided (default), it uses as
default values of LOD Score 3
and maximum recombination
fraction 0.50
.
Also, you can use the function suggest_lod
to calculate
a suggested LOD score considering that multiple tests are being
performed.
<- suggest_lod(mark_no_dist)
LOD_sug
LOD_sug#> [1] 3.639312
And apply this suggested value to the two-point tests:
<- group(mark_no_dist, LOD=LOD_sug)
LGs #> Selecting markers:
#> group 1
#> ........................
#> group 2
#> ................
#> group 3
#> ......
The previous command generates an object of class group
and the command print
for such object has two options. If
you type:
LGs
you will get detailed information about the groups, that is, all linkage groups will be printed, displaying the names of markers in each one of them.
#> This is an object of class 'group'
#> It was generated from the object "mark_no_dist"
#>
#> Criteria used to assign markers to groups:
#> LOD = 3.639312 , Maximum recombination fraction = 0.5
#>
#> No. markers: 49
#> No. groups: 3
#> No. linked markers: 49
#> No. unlinked markers: 0
#>
#> Printing groups:
#> Group 1 : 25 markers
#> M1 M2 M3 M5 M6 M10 M11 M12 M14 M15 M26 M28 M30 SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13
#>
#> Group 2 : 17 markers
#> M4 M9 M16 M20 M21 M23 M24 M27 M29 SNP17 SNP18 SNP20 SNP21 SNP22 SNP23 SNP24 SNP25
#>
#> Group 3 : 7 markers
#> M7 M8 M13 M18 M22 SNP14 SNP16
However, in case you just want to see some basic information (such as the number of groups, number of linked markers, etc), use:
print(LGs, detailed = FALSE)
#> This is an object of class 'group'
#> It was generated from the object "mark_no_dist"
#>
#> Criteria used to assign markers to groups:
#> LOD = 3.639312 , Maximum recombination fraction = 0.5
#>
#> No. markers: 49
#> No. groups: 3
#> No. linked markers: 49
#> No. unlinked markers: 0
You can notice that all markers are linked to some linkage group. If the LOD Score threshold is changed to a higher value, some markers are kept unassigned:
<- group(mark_no_dist, LOD = 6)
LGs #> Selecting markers:
#> group 1
#> ........................
#> group 2
#> ..........
#> group 3
#> .....
#> group 4
#> ....
LGs#> This is an object of class 'group'
#> It was generated from the object "mark_no_dist"
#>
#> Criteria used to assign markers to groups:
#> LOD = 6 , Maximum recombination fraction = 0.5
#>
#> No. markers: 49
#> No. groups: 4
#> No. linked markers: 47
#> No. unlinked markers: 2
#>
#> Printing groups:
#> Group 1 : 25 markers
#> M1 M2 M3 M5 M6 M10 M11 M12 M14 M15 M26 M28 M30 SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13
#>
#> Group 2 : 11 markers
#> M4 M9 M16 M20 M21 M23 M27 SNP17 SNP18 SNP20 SNP21
#>
#> Group 3 : 6 markers
#> M8 M13 M18 M22 SNP14 SNP16
#>
#> Group 4 : 5 markers
#> M24 SNP22 SNP23 SNP24 SNP25
#>
#> Unlinked markers: 2 markers
#> M7 M29
Changing back to the previous criteria, now setting the maximum recombination fraction to 0.40:
<- group(mark_no_dist, LOD = LOD_sug, max.rf = 0.4)
LGs #> Selecting markers:
#> group 1
#> ........................
#> group 2
#> ................
#> group 3
#> ......
LGs#> This is an object of class 'group'
#> It was generated from the object "mark_no_dist"
#>
#> Criteria used to assign markers to groups:
#> LOD = 3.639312 , Maximum recombination fraction = 0.4
#>
#> No. markers: 49
#> No. groups: 3
#> No. linked markers: 49
#> No. unlinked markers: 0
#>
#> Printing groups:
#> Group 1 : 25 markers
#> M1 M2 M3 M5 M6 M10 M11 M12 M14 M15 M26 M28 M30 SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13
#>
#> Group 2 : 17 markers
#> M4 M9 M16 M20 M21 M23 M24 M27 M29 SNP17 SNP18 SNP20 SNP21 SNP22 SNP23 SNP24 SNP25
#>
#> Group 3 : 7 markers
#> M7 M8 M13 M18 M22 SNP14 SNP16
(new!)
The other function for grouping is called group_upgma
.
It is an adapted version of MAPpoly grouping
function.
<- group_upgma(mark_no_dist, expected.groups = 5, inter = F)
LGs_upgma plot(LGs_upgma)
You can define the expected number of groups in the
expected.groups
argument and check how the markers are
split in the plotted dendrogram. Using argument inter=TRUE
you can change interactively the number of groups defined by the red
squares in the graphic.
Once marker assignment to linkage groups is finished, the mapping
step can take place. First of all, you must set the mapping function
that should be used to display the genetic map throughout the analysis.
You can choose between Kosambi
or Haldane
mapping functions. To use Haldane, type:
set_map_fun(type = "haldane")
To use Kosambi, type:
set_map_fun(type = "kosambi")
If you do not set one of these functions, the kosambi is used as default.
Now, you must define which linkage group will be mapped. In other
words, a linkage group must be extracted
from the object of
class group
or group.upgma
, in order to be
mapped. For simplicity, we will start here with the smallest one, which
is linkage group 3 (considering the group
function). This
can be easily done using the following code:
<- make_seq(LGs, 3)
LG3 # or
# LG3 <- make_seq(LGs_upgma,3)
The first argument (LGs
) is an object of class
group
or group.upgma
and the second is a
number indicating which linkage group will be extracted, according to
the results stored in object LGs
. The object
LG3
, generated by function make_seq
, is of
class sequence
, showing that this function can be used with
several types of objects.
If you type
LG3
you will see which markers are comprised in the sequence, and also that no parameters have been estimated so far.
#>
#> Number of markers: 7
#> Markers in the sequence:
#> M7 M8 M13 M18 M22 SNP14 SNP16
#>
#> Parameters not estimated.
To order these markers, one can use a two-point based algorithm such as Seriation (Buetow and Chakravarti, 1987), Rapid Chain Delineation (Doerge, 1996), Recombination Counting and Ordering (Van Os et al., 2005) and Unidirectional Growth (Tan and Fu, 2006):
<- seriation(LG3, hmm = FALSE))
(LG3_ser <- rcd(LG3, hmm = FALSE))
(LG3_rcd <- record(LG3, hmm = FALSE))
(LG3_rec <- ug(LG3, hmm = FALSE)) (LG3_ug
Argument hmm
defines if the function should run the HMM
chain multipoint approach to estimate the genetic distances given the
marker order provided by the two-points ordering algorithm. We set here
the argument hmm=FALSE
because we just want to obtain the
marker order. We are not yet estimating the genetic distances. We
suggest to use hmm=TRUE
only when you already decided which
order is the best because the HMM chain is the most computationally
intensive step in the map building. You can use
rf_graph_table
to check the ordering quality (see details
below) and make editions in the marker order using
drop_marker
. After, you can use map
or
map_avoid_unlinked
functions to estimate the genetic
distances (check session #Map-estimation-for-an-arbitrary-order).
Here we can check, there are some differences between each ordering algorithm results (results not shown).
Alternatively, you can also use the mds_onemap
function
to obtain a first draft for the order of the markers. The
mds_onemap
is a wrapper function that makes an interface
between OneMap
and MDSMap
package. The
ordering approach presented in MDSMap
provides a faster and
efficient way of ordering markers using multi-dimensional scaling. The
method also provides diagnostics graphics and parameters to find
outliers to help users to filter the dataset. You can find more
information in MDSMap
vignette. Here we
will show a simple example of how it can be used for ordering our
example markers from an outcrossing population.
<- mds_onemap(LG3, hmm = FALSE)
LG3_mds #> Stress: 0.225277102704001
#> Mean Nearest Neighbour Fit: 22.0492750538469
If you only specify the input sequence, mds_onemap will use the
default parameters. It will also generate an MDSMap input file in the
out.file
file. You can use out.file
in the
MDSMap package to try other parameters too. The default method used is
the principal curves, know more about using ?mds_onemap
and
reading the MDSMap vignette.
Besides these algorithms use a two-point approach to order the
markers, if you set hmm=TRUE
a multipoint approach is
applied to estimate the genetic distances after the order is estimated.
Thus, it can happen that some markers are not considered linked when
evaluated by multipoint information, and the function will return an
error like this:
ERROR: The linkage between markers 1 and 2 did not reach the OneMap default criteria. They are probably segregating independently
You can automatically remove these markers setting argument
rm_unlinked = TRUE
. The marker will be removed, and the
ordering algorithms will be restarted. Warning messages will inform
which markers were removed. If you don’t get warning messages, it means
that any marker needed to be removed. This is our case in this example,
but if you obtain an error or warning running your dataset, you already
know what happened.
NOTE: (new!) If
your sequence has many markers (more than 60), we suggest to first use
hmm=FALSE to check the ordering and after speed up mds
,
seriation
, rcd
, record
and
ug
using BatchMap parallelization approach. See section Speed up analysis with parallelization
for more information.
To order by comparing all possible orders (exhaustive search), the
function compare
can be used:
<- compare(LG3) LG3_comp
WARNING: This algorithm can take some time to run, depending on marker types in the linkage group. If you are working on a personal computer, without high capacity, we recommend using a maximum of ten markers.
If you have more markers in your group, we suggest using the
following explained approaches order_seq
.
In the example, LG3
contains only seven markers. Two of
them are of type D1, and one is segregating in 3:1 fashion (type C).
Thus, although the number of possible orders is relatively small (360),
for each order, there are various possible combinations of linkage
phases. Also, the convergence of the EM algorithm takes considerably
more time, because markers of type C and D are not very informative.
The first argument to the compare
function is an object
of class sequence
(the extracted group LG3
),
and the object generated by this function is of class
compare
.
To see the results of the previous step, type:
LG3_comp
NOTE: Check the GitHub vignette version to visualize the results.
Remember that for outcrossing populations, one needs to estimate
marker order and also linkage phases between markers for a given order.
However, because two-point analysis provides information about linkage
phases, this information is taken into consideration in the
compare
function, reducing the number of combinations to be
evaluated. If a given linkage phase has LOD greater than 0.005 in the
two-point analysis, we assume that this phase is very unlikely and so
does not need to be evaluated in the multipoint procedure used by
compare
. We did extensive simulations, which showed that
this is a good procedure.
By default, OneMap
stores 50 orders, which may or may
not be unique. The value of LOD
refers to the overall LOD
Score, considering all orders tested. Nested LOD
refers to
LOD Scores within a given order, that is, scores for different
combinations of linkage phases for the same marker order.
For example, order 1 has the largest value of log-likelihood and,
therefore, its LOD Score is zero for a given combination of linkage
phases (CC, CC, RR, RR). For this same order and other linkage phases,
LOD Score is -5.20. Analyzing the results for order 2, notice that its
highest LOD Score is very close to zero, indicating that this order is
also quite plausible. Notice also that Nested LOD
will
always contain at least one zero value, corresponding to the best
combination of phases for markers in a given order. Due to the
information provided by a two-point analysis, not all combinations are
tested, and that is the reason why the number of Nested LOD values is
different for each order.
Unless one has some biological information, it is a good idea to choose the order with the highest likelihood. The final map can then be obtained with the command.
<- make_seq(LG3_comp, 1, 1) LG3_final
The first argument is the object of class compare
. The
second argument indicates which order is chosen: 1 is for the order with
the highest likelihood, 2 is for the second-best, and so on. The third
argument indicates which combination of phases is chosen for a given
order: 1 also means the combination with the highest likelihood among
all combinations of phases (based on Nested LOD).
For simplicity, these values are defaults, so typing
<- make_seq(LG3_comp) LG3_final
have the same effect.
To see the final map, type:
LG3_final#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 43 SNP14 0.00 a | | b a | | a
#> 22 M22 13.51 a | | b a | | a
#> 7 M7 23.08 a | | a a | | b
#> 18 M18 64.51 a | | b c | | d
#> 8 M8 70.20 b | | a b | | a
#> 13 M13 72.86 a | | o a | | o
#> 44 SNP16 78.74 b | | a b | | a
#>
#> 7 markers log-likelihood: -398.6832
At the leftmost position, marker names are displayed.
Position
shows the cumulative distance using the Kosambi
mapping function. Finally, Parent 1
and
Parent 2
show the diplotypes of both parents, that is, the
combination in which alleles are arranged in the chromosomes, given the
estimated linkage phase. The notation is the same as that used by Wu et
al. (2002a). Details about how ordering algorithms can be chosen and
used are presented by Mollinari et al. (2009).
A careful examination of the results can be done using the function
rf_graph_table
to provide graphical view:
rf_graph_table(LG3_final)
With the default arguments, this function plots the recombination
fractions between the markers pointed in the axes. You can change the
number of colors from the rainbow
palette with the argument
n.colors
. Hot colors (more close to red) represent lower
values of recombination fractions, as shown in the scale at the right
side of the graphic. White cells indicate combinations of markers for
which the recombination fractions cannot be estimated (D1 and D2). If
you want to analyze the LOD values between the markers, use
graph.LOD = TRUE
.
rf_graph_table(LG3_final, graph.LOD = TRUE)
If you change the inter
argument to TRUE
,
you should also specify an output HTML file name in
html.file
. This HTML contains an iterative plot graphic. If
you hover the mouse cursor over the cells it shows some extra
information about cells, as a percentage of missing data, marker name,
and type. The output HTML file is generated in your work directory and
opens automatically in your internet browser.
rf_graph_table(LG3_final, inter = TRUE, mrk.axis= "names", html.file = "LG3.html")
For example, passing on the cell corresponding to markers
8
and 13
, you can see their names
(M8
and M13
), types (B3.7
and
C.8
), recombination fraction (rf = 0.03
) and
LOD Scores for each possible linkage phase. This is quite useful in
helping to interpret the results.
If you want to see corresponding marker numbers (not the names) in
the axis, just change the argument mrk.axis
to
numbers
. It can make the next steps easier.
rf_graph_table(LG3_final, inter = FALSE, mrk.axis = "numbers")
The rf_graph_table
can also be used to check the order
of markers based on the monotonicity of the matrix: as we get away from
the secondary diagonal, the recombination fraction values should
increase.
It is possible to see a gap between markers M7
and
M18
(numbers 7 and 18). In some cases, gaps could indicate
that the group must be divided at this position, but here
SNP18
(number 43) also shows linkage with M8
,
which points that probably it is only a gap. Adding more markers to
these groups could fill this gap.
Changing other arguments of the function, you can add/remove labels of the axes (‘lab.xy’) and add a title to the graph (‘main’).
rf_graph_table(LG3_final, main = "LG3", inter = FALSE, n.colors = 7, lab.xy = c("markers", "markers"))
Now, let us map the markers in linkage group number 2.
Again, extract
that group from the object
LGs
:
<- make_seq(LGs, 2)
LG2
LG2#>
#> Number of markers: 17
#> Markers in the sequence:
#> M4 M9 M16 M20 M21 M23 M24 M27 M29 SNP17 SNP18 SNP20 SNP21 SNP22 SNP23 SNP24
#> SNP25
#>
#> Parameters not estimated.
Note that there are more than 10 markers in this group, so it is
infeasible to use the compare
function with all of them
because it will take a very long time to proceed.
First, use rcd
to get a preliminary order estimate:
<- rcd(LG2, hmm = F)
LG2_rcd
LG2_rcd#>
#> Number of markers: 17
#> Markers in the sequence:
#> M27 SNP20 M16 M20 M4 M21 M23 SNP17 SNP18 M9 SNP21 SNP24 SNP23 M24 SNP22 SNP25
#> M29
#>
#> Parameters not estimated.
rf_graph_table(LG2_rcd)
Use the marker_type
function to check the segregation
types of all markers in this group:
marker_type(LG2)
#> Marker Marker.name Type
#> 1 4 M4 A.4
#> 2 9 M9 D1.10
#> 3 16 M16 D2.17
#> 4 20 M20 A.1
#> 5 21 M21 D2.16
#> 6 23 M23 C.8
#> 7 24 M24 B3.7
#> 8 27 M27 D1.12
#> 9 29 M29 D1.13
#> 10 45 SNP17 D1.10
#> 11 46 SNP18 D1.10
#> 12 47 SNP20 D1.10
#> 13 48 SNP21 D1.10
#> 14 49 SNP22 B3.7
#> 15 50 SNP23 B3.7
#> 16 51 SNP24 B3.7
#> 17 52 SNP25 B3.7
Based on their segregation types and distribution on the preliminary
map, markers M4
, M20
, M24
,
SNP22
, SNP23
, SNP24
and
SNP25
are the most informative ones (type A
is
better, followed by type B
). So, let us create a framework
of ordered markers using compare
for the most informative
ones:
<- make_seq(twopts, c(4, 20, 24, 49,50,51, 52)) LG2_init
Here there is a automatic way of obtain a new sequence only with markers selected by type: (new!)
<- seq_by_type(sequence = LG2, mk_type = c("A", "B"))
LG2_init marker_type(LG2_init)
#> Marker Marker.name Type
#> 1 4 M4 A.4
#> 2 20 M20 A.1
#> 3 24 M24 B3.7
#> 4 49 SNP22 B3.7
#> 5 50 SNP23 B3.7
#> 6 51 SNP24 B3.7
#> 7 52 SNP25 B3.7
# If I want to reduce even more the number of markers
# I can use drop_marker function
<- drop_marker(LG2_init, 52)
LG2_init marker_type(LG2_init)
#> Marker Marker.name Type
#> 1 4 M4 A.4
#> 2 20 M20 A.1
#> 3 24 M24 B3.7
#> 4 49 SNP22 B3.7
#> 5 50 SNP23 B3.7
#> 6 51 SNP24 B3.7
Now, the first argument to make_seq
is an object of
class rf_2pts
, and the second argument is a vector of
integers, specifying which molecular markers comprise the sequence.
<- compare(LG2_init)
LG2_comp #> Warning in compare_outcross(input.seq = input.seq, n.best = n.best, tol = tol, : This operation may take a VERY long time
Select the best order:
<- make_seq(LG2_comp) LG2_frame
Also, we can obtain a useful diagnostic graphic using the function
rf_graph_table
.
rf_graph_table(LG2_frame, mrk.axis = "numbers")
The graphic shows that there are two groups of markers, once
M20
and M4
are far from the other markers.
These markers could be in other linkage groups, or they are distant in
the same group. Adding more markers will give more information to solve
this issue.
Next, let us try to map the remaining markers, one at a time. First,
we will try to add the remaining most informative markers. Starting with
SNP25
:
<- try_seq(LG2_frame, 52) LG2_extend
LG2_extend#>
#> LOD scores correspond to the best linkage phase combination
#> for each position
#>
#> The symbol "*" outside the box indicates that more than one
#> linkage phase is possible for the corresponding position
#>
#>
#> Marker tested: 52
#>
#> Markers LOD
#> =====================
#> | |
#> | -26.44 | 1
#> | 20 |
#> | -50.47 | 2
#> | 4 |
#> | -8.34 | 3
#> | 51 |
#> | -8.79 | 4
#> | 50 |
#> | -8.43 | 5
#> | 24 |
#> | -2.46 | 6
#> | 49 |
#> | 0.00 | 7
#> | |
#> =====================
Based on the LOD Scores, marker SNP25 is probably better located
after SNP22
(number 49). Detailed results can be seen
with:
print(LG2_extend, 7)
#>
#> LOD is the overall LOD score (among all orders)
#>
#> NEST.LOD is the LOD score within the order
#>
#> Marker tested: 52
#> --------------
#> | | |
#> | 20 | |
#> | | CR |
#> | 4 | |
#> | | CC |
#> | 51 | |
#> | | CC |
#> | 50 | |
#> | | CC |
#> | 24 | |
#> | | CC |
#> | 49 | |
#> | | CC |
#> | 52 | |
#> | | |
#> |------------|
#> | LOD | 0.0|
#> |------------|
#> |NEST.| |
#> | LOD | 0.0|
#> --------------
The second argument indicates the position where to place the marker. Note that the first allele arrangement is the most likely one.
It should be pointed out that the framework created by the function
compare
(with M20
, M4
,
SNP24
, SNP23
, M24
and
SNP22
, or numbers 20
, 4
,
51
,50
, 24
and 49
)
could be in reverse order (SNP22
, M24
,
SNP23
, SNP24
, M4
and
M20
, or numbers 49
, 24
,
50
, 51
, 4
, 20
) and
still represent the same map. Thus, the positioning of markers with the
try_seq
command can be different on your computer. For
example, here marker SNP25
(number 52
) was
better placed at position 7; however, if you obtain a reversed order,
marker SNP25
would be better placed in position 1. In both
cases, the best position for this marker is after
SNP22
.
We can better evaluate the order with rf_graph_table
. It
requires an object of the sequence
class with mapping
information.
<- make_seq(LG2_extend, 7, 1) LG2_test
When using make_seq
with an object of class
try
, the second argument is the position on the map
(according to the scale on the right of the output) and the last
argument indicates linkage phases (defaults to 1, higher nested
LOD).
rf_graph_table(LG2_test, mrk.axis = "numbers")
We can see that SNP25
(or marker 52) was positioned at
the end of the sequence and the color pattern shows that it is strongly
linked with its neighbors, indicating that it is well-positioned. We
will maintain this marker at this position:
<- LG2_test LG2_frame
Adding other markers, one by one (output not shown):
<- try_seq(LG2_frame, 9)
LG2_extend <- make_seq(LG2_extend, 3)
LG2_frame <- try_seq(LG2_frame, 16)
LG2_extend <- make_seq(LG2_extend, 1)
LG2_frame <- try_seq(LG2_frame, 21)
LG2_extend <- make_seq(LG2_extend, 4)
LG2_frame <- try_seq(LG2_frame, 23)
LG2_extend <- make_seq(LG2_extend, 5)
LG2_frame <- try_seq(LG2_frame, 27)
LG2_extend <- make_seq(LG2_extend, 1)
LG2_frame <- try_seq(LG2_frame, 29)
LG2_extend <- make_seq(LG2_extend, 12)
LG2_frame <- try_seq(LG2_frame, 45)
LG2_extend <- make_seq(LG2_extend, 8)
LG2_frame <- try_seq(LG2_frame, 46)
LG2_extend <- make_seq(LG2_extend, 7)
LG2_frame <- try_seq(LG2_frame, 47)
LG2_extend <- make_seq(LG2_extend, 7)
LG2_frame <- try_seq(LG2_frame, 48)
LG2_extend <- make_seq(LG2_extend, 10) LG2_final
Checking graphically:
rf_graph_table(LG2_final)
The process of adding markers can be automated with the use of
function order_seq
.
<- order_seq(LG2, n.init = 5, THRES = 3)
LG2_ord #>
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|==== | 5%
|
|===== | 7%
|
|====== | 8%
|
|======= | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|============ | 17%
|
|============= | 18%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 23%
|
|================== | 25%
|
|=================== | 27%
|
|==================== | 28%
|
|===================== | 30%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 35%
|
|========================== | 37%
|
|=========================== | 38%
|
|============================ | 40%
|
|============================= | 42%
|
|============================== | 43%
|
|================================ | 45%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================== | 60%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================== | 65%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================= | 70%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|============================================================ | 85%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|=============================================================== | 90%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|======================================================================| 100%
This function automates what the try_seq
function does,
using some predefined rules. In the function, n.init = 5
means that five markers (the most informative ones) will be used in the
compare
step; THRES = 3
indicates that the
try_seq
step will only add markers to the sequence which
can be mapped with LOD Score greater than 3.
NOTE: Although very useful, this function can be misleading, especially if there are not many fully informative markers, so use it carefully. Results can vary between multiple runs on the same markers, of course.
Check the final order:
LG2_ord#>
#> Best sequence found.
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.76 a | | a b | | o
#> 20 M20 22.94 a | | b c | | d
#> 4 M4 35.71 a | | o o | | b
#> 23 M23 53.69 a | | o a | | o
#> 45 SNP17 72.86 a | | b a | | a
#> 50 SNP23 103.35 a | | b b | | a
#> 24 M24 108.37 a | | b b | | a
#> 49 SNP22 112.25 a | | b b | | a
#>
#> 9 markers log-likelihood: -551.4083
#>
#>
#>
#> The following markers could not be uniquely positioned.
#> Printing most likely positions for each unpositioned marker:
#>
#> ------------------------------------------------------
#> | | 9 | 21 | 29 | 46 | 47 | 48 | 51 | 52 |
#> |----|-----|-----|-----|-----|-----|-----|-----|-----|
#> | | | | | | | | | |
#> | 27 | | | | | | | | |
#> | | | | | | | | | |
#> | 16 | | | | | | | | |
#> | | | | | | | | | |
#> | 20 | | | | | | | | |
#> | | | | | | | | | |
#> | 4 | | | | | | | | |
#> | | | *** | | | | | | |
#> | 23 | | | | | | | | |
#> | | *** | * | | * | | *** | | |
#> | 45 | | | | | | | | |
#> | | * | | ** | *** | *** | | *** | |
#> | 50 | | | | | | | | |
#> | | | | | | | | | |
#> | 24 | | | | | | | | |
#> | | | | | | | | | |
#> | 49 | | | | | | | | |
#> | | | | *** | | | | ** | *** |
#> ------------------------------------------------------
#>
#> '***' indicates the most likely position(s) (LOD = 0.0)
#>
#> '**' indicates very likely positions (LOD > -1.0)
#>
#> '*' indicates likely positions (LOD > -2.0)
Note that markers 9
, 21
, 29
,
46
, 47
, 48
, 51
and
52
could not be safely mapped to a single position
(LOD Score > THRES
in absolute value). The output
displays the safe
order and the most likely positions for
markers not mapped, where ***
indicates the most likely
position and *
corresponds to other plausible
positions.
To get the safe order (i.e., without markers 9
,
21
, 29
, 46
, 47
,
48
, 51
and 52
), use
<- make_seq(LG2_ord, "safe") LG2_safe
and to get the order with all markers, use
<- make_seq(LG2_ord, "force")
LG2_all
LG2_all#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.76 a | | a b | | o
#> 20 M20 22.94 a | | b c | | d
#> 4 M4 35.71 a | | o o | | b
#> 21 M21 49.40 o | | o a | | b
#> 23 M23 55.18 a | | o a | | o
#> 48 SNP21 67.73 a | | b a | | a
#> 9 M9 74.32 a | | b a | | a
#> 46 SNP18 80.85 a | | b a | | a
#> 45 SNP17 96.14 a | | b a | | a
#> 47 SNP20 121.33 a | | b a | | a
#> 29 M29 170.70 o | | a o | | o
#> 50 SNP23 192.38 a | | b b | | a
#> 24 M24 197.52 a | | b b | | a
#> 49 SNP22 200.80 a | | b b | | a
#> 52 SNP25 205.70 a | | b b | | a
#> 51 SNP24 216.64 a | | b b | | a
#>
#> 17 markers log-likelihood: -871.9124
Notice that, for this linkage group, the forced
map
obtained with order_seq
is different from that obtained
with compare
plus try_seq
. It depends on which
markers we choose to try to add first when doing manually.
The order_seq
function can also perform two rounds of
the try_seq
algorithms, first using THRES
and
then THRES - 1
as a threshold. This generally results in
safe orders with more markers mapped but may take longer to run. To do
this, use the touchdown
option:
<- order_seq(LG2, n.init = 5, THRES = 3, touchdown = TRUE) LG2_ord
LG2_ord#>
#> Best sequence found.
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.76 a | | a b | | o
#> 20 M20 22.94 a | | b c | | d
#> 4 M4 35.71 a | | o o | | b
#> 23 M23 52.88 a | | o a | | o
#> 48 SNP21 65.21 a | | b a | | a
#> 45 SNP17 79.22 a | | b a | | a
#> 47 SNP20 103.39 a | | b a | | a
#> 50 SNP23 143.48 a | | b b | | a
#> 24 M24 148.50 a | | b b | | a
#> 49 SNP22 151.80 a | | b b | | a
#> 52 SNP25 156.80 a | | b b | | a
#>
#> 12 markers log-likelihood: -670.5
#>
#>
#>
#> The following markers could not be uniquely positioned.
#> Printing most likely positions for each unpositioned marker:
#>
#> ------------------------------------
#> | | 9 | 21 | 29 | 46 | 51 |
#> |----|-----|-----|-----|-----|-----|
#> | | | | | | |
#> | 27 | | | | | |
#> | | | | | | |
#> | 16 | | | | | |
#> | | | | | | |
#> | 20 | | | | | |
#> | | | | | | |
#> | 4 | | | | | |
#> | | | *** | | | |
#> | 23 | | | | | |
#> | | | | | | |
#> | 48 | | | | | |
#> | | *** | | | *** | |
#> | 45 | | | | | |
#> | | | | | * | |
#> | 47 | | | | | |
#> | | | | *** | | ** |
#> | 50 | | | | | |
#> | | | | | | |
#> | 24 | | | | | |
#> | | | | | | |
#> | 49 | | | | | |
#> | | | | | | |
#> | 52 | | | | | |
#> | | | | ** | | *** |
#> ------------------------------------
#>
#> '***' indicates the most likely position(s) (LOD = 0.0)
#>
#> '**' indicates very likely positions (LOD > -1.0)
#>
#> '*' indicates likely positions (LOD > -2.0)
For this particular sequence, the touchdown
step could
map safely markers 46
and 52
, but this depends
on the specific dataset.
rf_graph_table(LG2_all, mrk.axis = "numbers")
Finally, to check for alternative orders (because we did not use
exhaustive search), use the ripple_seq
function:
ripple_seq(LG2_all, ws = 4, LOD = LOD_sug)
#> 27-16-20-4-|-21-... OK
#>
#> ...-27-|-16-20-4-21-|-23-... OK
#>
#> ...-16-|-20-4-21-23-|-48-...
#> Alternative orders:
#> ... 16 20 4 21 23 48 ... : 0.00 ( linkage phases: ... 1 2 2 1 1 ... )
#> ... 16 20 4 23 21 48 ... : -2.21 ( linkage phases: ... 1 2 2 1 1 ... )
#>
#> ...-20-|-4-21-23-48-|-9-...
#> Alternative orders:
#> ... 20 4 21 23 48 9 ... : 0.00 ( linkage phases: ... 2 2 1 1 1 ... )
#> ... 20 4 23 21 48 9 ... : -2.21 ( linkage phases: ... 2 2 1 1 1 ... )
#>
#> ...-4-|-21-23-48-9-|-46-...
#> Alternative orders:
#> ... 4 21 23 48 9 46 ... : 0.00 ( linkage phases: ... 2 1 1 1 1 ... )
#> ... 4 23 21 48 9 46 ... : -2.21 ( linkage phases: ... 2 1 1 1 1 ... )
#>
#> ...-21-|-23-48-9-46-|-45-...
#> Alternative orders:
#> ... 21 23 48 9 46 45 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 21 23 48 46 9 45 ... : -1.05 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-23-|-48-9-46-45-|-47-...
#> Alternative orders:
#> ... 23 48 9 46 45 47 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 48 46 9 45 47 ... : -1.05 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 48 45 9 46 47 ... : -1.80 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 48 9 45 46 47 ... : -3.21 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 45 48 9 46 47 ... : -3.58 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-48-|-9-46-45-47-|-29-...
#> Alternative orders:
#> ... 48 9 46 45 47 29 ... : 0.00 ( linkage phases: ... 1 1 1 1 3 ... )
#> ... 48 46 9 45 47 29 ... : -1.05 ( linkage phases: ... 1 1 1 1 3 ... )
#> ... 48 45 9 46 47 29 ... : -1.80 ( linkage phases: ... 1 1 1 1 3 ... )
#> ... 48 9 45 46 47 29 ... : -3.21 ( linkage phases: ... 1 1 1 1 3 ... )
#>
#> ...-9-|-46-45-47-29-|-50-...
#> Alternative orders:
#> ... 9 46 45 47 29 50 ... : 0.00 ( linkage phases: ... 1 1 1 3 4 ... )
#> ... 9 45 46 47 29 50 ... : -3.21 ( linkage phases: ... 1 1 1 3 4 ... )
#>
#> ...-46-|-45-47-29-50-|-24-... OK
#>
#> ...-45-|-47-29-50-24-|-49-... OK
#>
#> ...-47-|-29-50-24-49-|-52-... OK
#>
#> ...-29-|-50-24-49-52-|-51-...
#> Alternative orders:
#> ... 29 50 24 49 52 51 : 0.00 ( linkage phases: ... 4 1 1 1 1 )
#> ... 29 52 49 24 50 51 : -0.90 ( linkage phases: ... 4 1 1 1 1 )
#> ... 29 50 24 52 49 51 : -2.76 ( linkage phases: ... 4 1 1 1 1 )
#> ... 29 50 52 49 24 51 : -2.79 ( linkage phases: ... 4 1 1 1 1 )
#> ... 29 50 49 52 24 51 : -3.31 ( linkage phases: ... 4 1 1 1 1 )
#> ... 29 49 52 24 50 51 : -3.31 ( linkage phases: ... 4 1 1 1 1 )
#>
#> 50-|-24-49-52-51
#> Alternative orders:
#> ... 50 24 49 52 51 : 0.00 ( linkage phases: ... 1 1 1 1 )
#> ... 50 24 52 49 51 : -2.76 ( linkage phases: ... 1 1 1 1 )
#> ... 50 52 49 24 51 : -2.79 ( linkage phases: ... 1 1 1 1 )
#> ... 50 51 24 49 52 : -2.85 ( linkage phases: ... 1 1 1 1 )
#> ... 50 49 52 24 51 : -3.31 ( linkage phases: ... 1 1 1 1 )
We should do this to any of the orders we found, either using
try_seq
or order_seq
. Here, we choose
LG2_all
for didactic purposes only. The second argument,
ws = 4
, means that subsets (windows) of four markers will
be permuted sequentially (4!
orders for each window), to
search for other plausible orders. The LOD
argument means
that only orders with LOD Score smaller than 3.68 will be printed.
The output shows sequences of four numbers, because
ws = 4
. They are followed by an OK
if there is
no alternative order with LOD Score smaller than
LOD = LOD_sug
in absolute value, or by a list of
alternative orders. In the example, some sequences showed alternative
orders with LOD smaller than LOD = LOD_sug
. However, the
best order was the original one (LOD = 0.00
).
If there were an alternative order more likely than the original, one should check the difference between these orders (and linkage phases).
In some cases, even if there are no better alternative orders
suggested by ripple_seq
, the graphic showed a color pattern
different from the expected. Then, we can remove doubtful markers (for
this groups markers 23
and 29
) and try to
position them again. First, we use the function drop_marker
to remove the selected marker of our sequence.
<- drop_marker(LG2_all, c(23,29)) LG2_test_seq
The function will provide a sequence with the same order as the
estimated map (LG2_all
). After, we should estimate the map
again using this predefined order (see section Map estimation for an arbitrary order
for further information). For this we use the map
function:
<- map(LG2_test_seq)) (LG2_test_map
Warning: If you find an error message like:
Error in as_mapper(.f, ...) : argument ".f" is missing, with no default
It’s because the map
function has a very common name,
and you can have in your environment other functions with the same name.
In the case of the pointed error, R is using the map
function from purrr
package instead of OneMap
,
to solve this, simply specify that you want the OneMap
function with ::
command from stringr
package:
library(stringr)
<- onemap::map(LG2_test_seq))
(LG2_test_map #>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.75 a | | a o | | b
#> 20 M20 22.93 a | | b d | | c
#> 4 M4 35.70 a | | o b | | o
#> 21 M21 52.67 o | | o b | | a
#> 48 SNP21 53.59 a | | b a | | a
#> 9 M9 60.40 a | | b a | | a
#> 46 SNP18 66.83 a | | b a | | a
#> 45 SNP17 81.83 a | | b a | | a
#> 47 SNP20 106.46 a | | b a | | a
#> 50 SNP23 145.68 a | | b a | | b
#> 24 M24 150.69 a | | b a | | b
#> 49 SNP22 153.97 a | | b a | | b
#> 52 SNP25 158.87 a | | b a | | b
#> 51 SNP24 169.81 a | | b a | | b
#>
#> 15 markers log-likelihood: -780.5746
NOTE: (new!) If
your sequence has many markers (more than 60), we suggest to speed up
map
using BatchMap parallelization approach. See section Speed up analysis with parallelization
for more information.
Now, we have the map without markers 23
and
51
.
rf_graph_table(LG2_test_map, mrk.axis = "numbers")
We use the try_seq
function to positioned them
again:
<- try_seq(LG2_test_map, 23)
LG2_test_seq <- make_seq(LG2_test_seq, 6)
LG2_test_23 <- try_seq(LG2_test_23, 29)
LG2_test_seq <- make_seq(LG2_test_seq, 12)
LG2_test_23_29 rf_graph_table(LG2_test_23_29, mrk.axis = "numbers")
Marker 23
kept its previous position, but marker
29
was re-positioned, configuring now a gap between markers
47
and 50
. We removed marker 29
from our map because its color pattern is too different from expected.
Then, our final map is:
<- LG2_test_23
LG2_final
LG2_final#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.75 a | | a o | | b
#> 20 M20 22.93 a | | b d | | c
#> 4 M4 35.70 a | | o b | | o
#> 21 M21 50.34 o | | o b | | a
#> 23 M23 54.87 a | | o o | | a
#> 48 SNP21 70.11 a | | b a | | a
#> 9 M9 76.77 a | | b a | | a
#> 46 SNP18 83.19 a | | b a | | a
#> 45 SNP17 98.20 a | | b a | | a
#> 47 SNP20 122.83 a | | b a | | a
#> 50 SNP23 161.89 a | | b a | | b
#> 24 M24 166.89 a | | b a | | b
#> 49 SNP22 170.18 a | | b a | | b
#> 52 SNP25 175.08 a | | b a | | b
#> 51 SNP24 186.02 a | | b a | | b
#>
#> 16 markers log-likelihood: -811.7886
rf_graph_table(LG2_final, mrk.axis = "numbers")
Finally, linkage group 1 (the largest one) will be analyzed. Extract markers:
<- make_seq(LGs, 1) LG1
Construct the linkage map, by automatically using the try algorithm:
<- order_seq(LG1, n.init = 6, touchdown = TRUE)
LG1_ord #> Warning in compare_outcross(input.seq = input.seq, n.best = n.best, tol = tol, : This operation may take a VERY long time
Notice that the second round of try_seq
added markers
10
, 31
, 32
, 35
,
36
, and 40
.
LG1_ord#>
#> Best sequence found.
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 34 SNP5 0.00 a | | b a | | b
#> 32 SNP2 13.48 a | | b a | | b
#> 12 M12 17.26 b | | a c | | a
#> 31 SNP1 22.81 a | | b a | | b
#> 3 M3 46.09 o | | a o | | o
#> 14 M14 59.22 a | | o b | | o
#> 2 M2 66.93 o | | o o | | a
#> 36 SNP7 92.49 b | | a b | | a
#> 35 SNP6 105.45 b | | a b | | a
#> 1 M1 109.86 b | | a b | | a
#> 10 M10 115.37 a | | a o | | b
#> 28 M28 132.05 a | | o b | | o
#> 26 M26 164.42 a | | b d | | c
#> 42 SNP13 179.54 b | | a b | | a
#> 6 M6 185.89 b | | a b | | a
#> 39 SNP10 192.06 b | | a b | | a
#> 41 SNP12 208.53 b | | a b | | a
#> 40 SNP11 225.65 b | | a b | | a
#> 15 M15 267.89 o | | a o | | b
#>
#> 19 markers log-likelihood: -1313.595
#>
#>
#>
#> The following markers could not be uniquely positioned.
#> Printing most likely positions for each unpositioned marker:
#>
#> ------------------------------------------
#> | | 5 | 11 | 30 | 33 | 37 | 38 |
#> |----|-----|-----|-----|-----|-----|-----|
#> | | | | | | | |
#> | 34 | | | | | | |
#> | | | | | *** | | |
#> | 32 | | | | | | |
#> | | | | *** | * | | |
#> | 12 | | | | | | |
#> | | | | | | | |
#> | 31 | | | | | | |
#> | | | | | | | |
#> | 3 | | | | | | |
#> | | | | | | | |
#> | 14 | | | | | | |
#> | | | | | | | |
#> | 2 | | | | | | |
#> | | | | | | * | |
#> | 36 | | | | | | |
#> | | | | | | ** | * |
#> | 35 | | | | | | |
#> | | | | | | | |
#> | 1 | | | | | | |
#> | | | | | | *** | *** |
#> | 10 | | | | | | |
#> | | | | | | | |
#> | 28 | | | | | | |
#> | | | | | | | |
#> | 26 | | | | | | |
#> | | | | | | | |
#> | 42 | | | | | | |
#> | | | | | | | |
#> | 6 | | | | | | |
#> | | | | | | | |
#> | 39 | | | | | | |
#> | | | | | | | |
#> | 41 | | | | | | |
#> | | | | | | | |
#> | 40 | | | | | | |
#> | | ** | ** | | | | |
#> | 15 | | | | | | |
#> | | *** | *** | | | | |
#> ------------------------------------------
#>
#> '***' indicates the most likely position(s) (LOD = 0.0)
#>
#> '**' indicates very likely positions (LOD > -1.0)
#>
#> '*' indicates likely positions (LOD > -2.0)
Now, get the order with all markers:
<- make_seq(LG1_ord, "force"))
(LG1_frame #>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 34 SNP5 0.00 a | | b a | | b
#> 33 SNP3 12.27 a | | b a | | b
#> 32 SNP2 18.45 a | | b a | | b
#> 30 M30 21.21 a | | b a | | b
#> 12 M12 22.22 b | | a c | | a
#> 31 SNP1 27.77 a | | b a | | b
#> 3 M3 51.05 o | | a o | | o
#> 14 M14 64.18 a | | o b | | o
#> 2 M2 72.05 o | | o o | | a
#> 36 SNP7 95.86 b | | a b | | a
#> 37 SNP8 112.04 b | | a b | | a
#> 35 SNP6 126.75 b | | a b | | a
#> 1 M1 131.15 b | | a b | | a
#> 38 SNP9 137.28 b | | a b | | a
#> 10 M10 148.79 a | | a o | | b
#> 28 M28 165.35 a | | o b | | o
#> 26 M26 197.72 a | | b d | | c
#> 42 SNP13 212.84 b | | a b | | a
#> 6 M6 219.19 b | | a b | | a
#> 39 SNP10 225.36 b | | a b | | a
#> 41 SNP12 241.82 b | | a b | | a
#> 40 SNP11 258.95 b | | a b | | a
#> 15 M15 301.19 o | | a o | | b
#> 5 M5 306.20 o | | o o | | a
#> 11 M11 326.21 o | | o b | | a
#>
#> 25 markers log-likelihood: -1554.628
Check the map graphically:
rf_graph_table(LG1_frame, mrk.axis = "numbers")
Check for alternative orders:
ripple_seq(LG1_frame)
#> 34-33-32-30-|-12-...
#> Alternative orders:
#> 34 33 32 30 12 ... : 0.00 ( linkage phases: 1 1 1 4 ... )
#> 34 32 33 30 12 ... : -1.29 ( linkage phases: 1 1 1 4 ... )
#> 34 33 30 32 12 ... : -1.55 ( linkage phases: 1 1 1 4 ... )
#>
#> ...-34-|-33-32-30-12-|-31-...
#> Alternative orders:
#> 34 33 12 30 32 31 ... : 0.00 ( linkage phases: 1 4 4 1 1 ... )
#> 34 33 30 12 32 31 ... : -0.25 ( linkage phases: 1 1 4 4 1 ... )
#> 34 33 32 30 12 31 ... : -0.59 ( linkage phases: 1 1 1 4 4 ... )
#> 34 32 33 30 12 31 ... : -1.88 ( linkage phases: 1 1 1 4 4 ... )
#> 34 33 30 32 12 31 ... : -2.14 ( linkage phases: 1 1 1 4 4 ... )
#> 34 33 32 12 30 31 ... : -2.74 ( linkage phases: 1 1 4 4 1 ... )
#>
#> ...-33-|-32-30-12-31-|-3-...
#> Alternative orders:
#> ... 33 12 30 32 31 3 ... : 0.00 ( linkage phases: ... 4 4 1 1 4 ... )
#> ... 33 30 12 32 31 3 ... : -0.25 ( linkage phases: ... 1 4 4 1 4 ... )
#> ... 33 32 30 12 31 3 ... : -0.59 ( linkage phases: ... 1 1 4 4 4 ... )
#> ... 33 30 32 12 31 3 ... : -2.14 ( linkage phases: ... 1 1 4 4 4 ... )
#> ... 33 32 12 30 31 3 ... : -2.74 ( linkage phases: ... 1 4 4 1 4 ... )
#>
#> ...-32-|-30-12-31-3-|-14-...
#> Alternative orders:
#> ... 32 30 12 31 3 14 ... : 0.00 ( linkage phases: ... 1 4 4 4 4 ... )
#> ... 32 12 30 31 3 14 ... : -2.15 ( linkage phases: ... 4 4 1 4 4 ... )
#>
#> ...-30-|-12-31-3-14-|-2-... OK
#>
#> ...-12-|-31-3-14-2-|-36-... OK
#>
#> ...-31-|-3-14-2-36-|-37-... OK
#>
#> ...-3-|-14-2-36-37-|-35-...
#> Alternative orders:
#> ... 3 14 2 36 37 35 ... : 0.00 ( linkage phases: ... 4 2 3 1 1 ... )
#> ... 3 14 2 37 36 35 ... : -0.89 ( linkage phases: ... 4 2 3 1 1 ... )
#>
#> ...-14-|-2-36-37-35-|-1-...
#> Alternative orders:
#> ... 14 2 36 37 35 1 ... : 0.00 ( linkage phases: ... 2 3 1 1 1 ... )
#> ... 14 2 37 36 35 1 ... : -0.89 ( linkage phases: ... 2 3 1 1 1 ... )
#>
#> ...-2-|-36-37-35-1-|-38-...
#> Alternative orders:
#> ... 2 36 37 1 35 38 ... : 0.00 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 37 35 1 38 ... : -0.25 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 37 36 35 1 38 ... : -1.14 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 37 36 1 35 38 ... : -2.34 ( linkage phases: ... 3 1 1 1 1 ... )
#>
#> ...-36-|-37-35-1-38-|-10-...
#> Alternative orders:
#> ... 36 37 1 35 38 10 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 35 1 38 10 ... : -0.25 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 35 1 38 37 10 ... : -1.23 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 38 35 1 10 ... : -1.35 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 38 35 1 37 10 ... : -1.42 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 38 1 35 10 ... : -1.54 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 38 1 35 37 10 ... : -1.69 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 35 38 1 37 10 ... : -2.45 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 1 35 38 37 10 ... : -2.50 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 1 38 35 10 ... : -2.78 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 35 38 1 10 ... : -2.85 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-37-|-35-1-38-10-|-28-...
#> Alternative orders:
#> ... 37 1 35 38 10 28 ... : 0.00 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 37 35 1 38 10 28 ... : -0.25 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 37 38 35 1 10 28 ... : -1.35 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 37 38 1 35 10 28 ... : -1.54 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 37 1 38 35 10 28 ... : -2.78 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 37 35 38 1 10 28 ... : -2.85 ( linkage phases: ... 1 1 1 1 4 ... )
#>
#> ...-35-|-1-38-10-28-|-26-...
#> Alternative orders:
#> ... 35 1 38 10 28 26 ... : 0.00 ( linkage phases: ... 1 1 1 4 2 ... )
#> ... 35 38 1 10 28 26 ... : -2.60 ( linkage phases: ... 1 1 1 4 2 ... )
#>
#> ...-1-|-38-10-28-26-|-42-... OK
#>
#> ...-38-|-10-28-26-42-|-6-... OK
#>
#> ...-10-|-28-26-42-6-|-39-...
#> Alternative orders:
#> ... 10 28 26 42 6 39 ... : 0.00 ( linkage phases: ... 4 2 3 1 1 ... )
#> ... 10 28 26 6 42 39 ... : -2.30 ( linkage phases: ... 4 2 3 1 1 ... )
#>
#> ...-28-|-26-42-6-39-|-41-...
#> Alternative orders:
#> ... 28 26 42 6 39 41 ... : 0.00 ( linkage phases: ... 2 3 1 1 1 ... )
#> ... 28 26 6 42 39 41 ... : -2.30 ( linkage phases: ... 2 3 1 1 1 ... )
#>
#> ...-26-|-42-6-39-41-|-40-...
#> Alternative orders:
#> ... 26 42 6 39 41 40 ... : 0.00 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 26 6 42 39 41 40 ... : -2.30 ( linkage phases: ... 3 1 1 1 1 ... )
#>
#> ...-42-|-6-39-41-40-|-15-... OK
#>
#> ...-6-|-39-41-40-15-|-5-... OK
#>
#> ...-39-|-41-40-15-5-|-11-...
#> Alternative orders:
#> ... 39 41 40 15 5 11 : 0.00 ( linkage phases: ... 1 1 1 1 1 )
#> ... 39 41 40 5 15 11 : -2.57 ( linkage phases: ... 1 1 1 1 1 )
#>
#> 41-|-40-15-5-11
#> Alternative orders:
#> ... 41 40 15 5 11 : 0.00 ( linkage phases: ... 1 1 1 1 )
#> ... 41 40 11 5 15 : -0.82 ( linkage phases: ... 1 1 1 1 )
#> ... 41 40 5 15 11 : -2.57 ( linkage phases: ... 1 1 1 1 )
#> ... 41 40 11 15 5 : -2.58 ( linkage phases: ... 1 1 1 1 )
No better order was observed.
Let’s check how it behaves with the MDS approach:
<- mds_onemap(LG1, rm_unlinked = TRUE, hmm = F)
LG1_mds #> Stress: 0.303935806009324
#> Mean Nearest Neighbour Fit: 30.2669402204845
rf_graph_table(LG1_mds)
Based on the drafts from order_seq
or/and
mds_onemap
, we can remove some doubtful markers accordingly
with the graphic, try to position them again, and decide if and where we
will maintain them.
<- drop_marker(LG1_frame, c(10,11,28,42))
LG1_test_seq <- onemap::map(LG1_test_seq) LG1_test_map
rf_graph_table(LG1_test_map, mrk.axis = "numbers")
<- try_seq(LG1_test_map,10)
LG1_extend <- make_seq(LG1_extend,15)
LG1_test_map <- try_seq(LG1_test_map,11)
LG1_extend <- make_seq(LG1_extend,23) # We choose to remove this marker
LG1_test <- try_seq(LG1_test_map,28)
LG1_extend <- make_seq(LG1_extend,16)
LG1_test_map <- try_seq(LG1_test_map,42)
LG1_extend <- make_seq(LG1_extend,17) LG1_final
Print it:
LG1_final#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 34 SNP5 0.00 a | | b a | | b
#> 33 SNP3 12.27 a | | b a | | b
#> 32 SNP2 18.45 a | | b a | | b
#> 30 M30 21.21 a | | b a | | b
#> 12 M12 22.22 b | | a c | | a
#> 31 SNP1 27.77 a | | b a | | b
#> 3 M3 51.05 o | | a o | | o
#> 14 M14 64.18 a | | o b | | o
#> 2 M2 72.10 o | | o o | | a
#> 36 SNP7 95.93 b | | a b | | a
#> 37 SNP8 112.12 b | | a b | | a
#> 35 SNP6 126.83 b | | a b | | a
#> 1 M1 131.23 b | | a b | | a
#> 38 SNP9 137.34 b | | a b | | a
#> 10 M10 150.29 a | | a o | | b
#> 28 M28 166.58 a | | o b | | o
#> 42 SNP13 203.65 b | | a b | | a
#> 26 M26 219.73 a | | b d | | c
#> 6 M6 228.71 b | | a b | | a
#> 39 SNP10 234.88 b | | a b | | a
#> 41 SNP12 251.25 b | | a b | | a
#> 40 SNP11 268.27 b | | a b | | a
#> 15 M15 309.06 o | | a o | | b
#> 5 M5 314.07 o | | o o | | a
#>
#> 24 markers log-likelihood: -1518.748
rf_graph_table(LG1_final)
As an option, different algorithms to order markers could be applied:
<- seriation(LG1, hmm = F)
LG1_ser rf_graph_table(LG1_ser)
<- rcd(LG1, hmm = F)
LG1_rcd rf_graph_table(LG1_rcd)
<- record(LG1, hmm = F)
LG1_rec rf_graph_table(LG1_rec)
<- ug(LG1, hmm = F)
LG1_ug rf_graph_table(LG1_ug)
There are some differences between the results. Seriation did not provide good results in this case. See Mollinari et al. (2009) for an evaluation of these methods.
In our example, we have reference genome chromosome and position information for some of the markers; here, we will exemplify one method of using this information to help build the genetic map.
With the CHROM
information in the input file, you can
identify markers belonging to some chromosome using the function
make_seq
with the rf_2pts
object. For example,
assign the string "1"
for the second argument to get
chromosome 1 makers. The output sequence will be automatically ordered
by POS
information.
<- make_seq(twopts, "1")
CHR1
CHR1#>
#> Number of markers: 12
#> Markers in the sequence:
#> SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13
#>
#> Parameters not estimated.
<- make_seq(twopts, "2")
CHR2 <- make_seq(twopts, "3") CHR3
Here we use the string "1"
because it is our chromosome
ID, you can have a different string as ID, check this with:
unique(bins_example$CHROM)
#> [1] NA "1" "2" "3"
We can see that we have markers without chromosome information
(NA
) and markers with chromosome ID "1"
,
"2"
and "3"
.
According to CHROM
information we have three defined
linkage groups, now we can try to group the markers without chromosome
information to them using recombination information. For this, we can
use the function group_seq
:
<- group_seq(input.2pts = twopts, seqs = "CHROM", unlink.mks = mark_no_dist,
CHR_mks repeated = FALSE)
#> Selecting markers:
#> group 1
#> ........................
#> group 2
#> ........
#> group 3
#> ....
#> Selecting markers:
#> group 1
#> ......
#> group 2
#> ............
#> group 3
#> ........
#> Selecting markers:
#> group 1
#> ................
#> group 2
#> ............
#> group 3
#> ....
The function works as the function group
but considering
pre-existing sequences. Setting seqs
argument with the
string "CHROM"
, it will consider the pre-existing sequences
according to CHROM
information. You can also indicate other
pre-existing sequences if they make sense for your study. For that, you
should inform a list of objects of class sequences
, as the
example:
<- group_seq(input.2pts = twopts, seqs = list(CHR1=CHR1, CHR2=CHR2, CHR3=CHR3),
CHR_mks unlink.mks = mark_no_dist, repeated = FALSE)
In this case, the command had the same effect as the previous because we indicate chromosome sequences, but other sequences can be used.
The unlink.mks
argument receives an object of class
sequence
; this defines which markers will be tested to
group with the sequences in seqs
. In our example, we will
indicate only the markers with no segregation distortion, using the
sequence mark_no_dist
. It is also possible to use the
string "all"
to test all the remaining markers at the
rf_2pts
object.
In some cases, the same marker can group into more than one sequence;
those markers will be considered repeated
. We can choose if
we want to remove or not (FALSE/TRUE
) them of the output
sequences, with the argument repeated
. Anyway, their
numbers will be informed at the list repeated
in the output
object. In the example case, there are no repeated markers. However, if
they exist, it could indicate that their groups actually constitute the
same group. Also, genotyping errors can generate repeated markers.
Anyway, they deserve better investigations.
We can access detailed information about the results, just printing:
CHR_mks#> This is an object of class 'group_seq'
#> Criteria used to assign markers to groups:
#> LOD = 3 , Maximum recombination fraction = 0.5
#>
#> No. markers in input sequences:
#> CHR1 : 12 markers
#> CHR2 : 2 markers
#> CHR3 : 8 markers
#>
#> No. unlinked input markers: 27 markers
#>
#> No. markers in output sequences:
#> CHR1 : 25 markers
#> CHR2 : 7 markers
#> CHR3 : 17 markers
#> No. unlinked: 0 markers
#> No. repeated: 0 markers
#>
#> Printing output sequences:
#> Group CHR1 : 25 markers
#> SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13 M1 M2 M3 M5 M6 M10 M11 M12 M14 M15 M26 M28 M30
#>
#> Group CHR2 : 7 markers
#> SNP14 SNP16 M7 M8 M13 M18 M22
#>
#> Group CHR3 : 17 markers
#> SNP17 SNP18 SNP20 SNP21 SNP22 SNP23 SNP24 SNP25 M4 M9 M16 M20 M21 M23 M24 M27 M29
#>
#> Unlinked markers: 0 markers
#>
#>
#> Repeated markers: 0 markers
#>
Also, we can access the numbers of repeated markers with:
$repeated
CHR_mks#> [1] NA
We have no repeated markers.
In the same way, we can access the output sequences:
$sequences$CHR1
CHR_mks#>
#> Number of markers: 25
#> Markers in the sequence:
#> SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13 M1 M2 M3 M5 M6
#> M10 M11 M12 M14 M15 M26 M28 M30
#>
#> Parameters not estimated.
# or
$sequences[[1]]
CHR_mks#>
#> Number of markers: 25
#> Markers in the sequence:
#> SNP1 SNP2 SNP3 SNP5 SNP6 SNP7 SNP8 SNP9 SNP10 SNP11 SNP12 SNP13 M1 M2 M3 M5 M6
#> M10 M11 M12 M14 M15 M26 M28 M30
#>
#> Parameters not estimated.
For this function, optional arguments are LOD
and
max.rf
, which define thresholds to be used when assigning
markers to linkage groups. If none is provided (default), criteria
previously defined for the object rf_2pts
are used.
Now we can order the markers in each group as we made before in (Genetic mapping of linkage group 1,2 and 3). As shown, we can choose different approaches to order the markers.
To order those groups, first, we will use the order_seq
function to access a preliminary order, and after, we will edit some
marker’s positions or remove some of them according to their color
pattern in the rf_graph_table
graphic, and other parameters
as likelihood and map size.
# CHR1_frame <- mds_onemap(CHR_mks$sequences$CHR1, hmm = F)
# or
<- order_seq(CHR_mks$sequences$CHR1)
CHR1_ord <- make_seq(CHR1_ord, "force") CHR1_frame
rf_graph_table(CHR1_frame) # graphic not shown
The group is similar to that built before with only recombinations
information. We will better explore differences in the later step. Only
marker 11
does not follow the expected color pattern; then,
we will try to reposition it.
<- drop_marker(CHR1_frame, 11)
CHR1_test_seq <- onemap::map(CHR1_test_seq)
CHR1_test_map <- try_seq(CHR1_test_map, 11)
CHR1_add11_seq <- make_seq(CHR1_add11_seq, 25) # marker 11 was placed at the same position as before CHR1_add11
Based in those results, we decide not to include marker 11 in our map.
CHR1_test_map#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 34 SNP5 0.00 a | | b a | | b
#> 33 SNP3 12.27 a | | b a | | b
#> 30 M30 15.58 a | | b a | | b
#> 12 M12 16.59 b | | a c | | a
#> 32 SNP2 20.15 a | | b a | | b
#> 31 SNP1 28.31 a | | b a | | b
#> 3 M3 51.17 o | | a o | | o
#> 14 M14 64.54 a | | o b | | o
#> 2 M2 72.44 o | | o o | | a
#> 36 SNP7 96.32 b | | a b | | a
#> 35 SNP6 109.13 b | | a b | | a
#> 1 M1 113.53 b | | a b | | a
#> 38 SNP9 119.59 b | | a b | | a
#> 37 SNP8 136.22 b | | a b | | a
#> 10 M10 149.87 a | | a o | | b
#> 28 M28 167.09 a | | o b | | o
#> 26 M26 199.46 a | | b d | | c
#> 42 SNP13 214.58 b | | a b | | a
#> 6 M6 220.93 b | | a b | | a
#> 39 SNP10 227.10 b | | a b | | a
#> 41 SNP12 243.56 b | | a b | | a
#> 40 SNP11 260.69 b | | a b | | a
#> 15 M15 302.93 o | | a o | | b
#> 5 M5 307.95 o | | o o | | a
#>
#> 24 markers log-likelihood: -1507.316
rf_graph_table(CHR1_test_map)
<- CHR1_test_map CHR1_final
Checking for better orders:
ripple_seq(CHR1_final)
#> 34-33-30-12-|-32-...
#> Alternative orders:
#> 34 33 12 30 32 ... : 0.00 ( linkage phases: 1 4 4 1 ... )
#> 34 33 30 12 32 ... : -0.25 ( linkage phases: 1 1 4 4 ... )
#>
#> ...-34-|-33-30-12-32-|-31-...
#> Alternative orders:
#> 34 33 12 30 32 31 ... : 0.00 ( linkage phases: 1 4 4 1 1 ... )
#> 34 33 30 12 32 31 ... : -0.25 ( linkage phases: 1 1 4 4 1 ... )
#> 34 33 32 30 12 31 ... : -0.59 ( linkage phases: 1 1 1 4 4 ... )
#> 34 32 33 30 12 31 ... : -1.88 ( linkage phases: 1 1 1 4 4 ... )
#> 34 33 30 32 12 31 ... : -2.14 ( linkage phases: 1 1 1 4 4 ... )
#> 34 33 32 12 30 31 ... : -2.74 ( linkage phases: 1 1 4 4 1 ... )
#>
#> ...-33-|-30-12-32-31-|-3-...
#> Alternative orders:
#> ... 33 12 30 32 31 3 ... : 0.00 ( linkage phases: ... 4 4 1 1 3 ... )
#> ... 33 30 12 32 31 3 ... : -0.25 ( linkage phases: ... 1 4 4 1 3 ... )
#> ... 33 32 30 12 31 3 ... : -0.59 ( linkage phases: ... 1 1 4 4 3 ... )
#> ... 33 30 32 12 31 3 ... : -2.14 ( linkage phases: ... 1 1 4 4 3 ... )
#> ... 33 32 12 30 31 3 ... : -2.74 ( linkage phases: ... 1 4 4 1 3 ... )
#>
#> ...-30-|-12-32-31-3-|-14-...
#> Alternative orders:
#> ... 30 12 32 31 3 14 ... : 0.00 ( linkage phases: ... 4 4 1 3 3 ... )
#> ... 30 32 12 31 3 14 ... : -1.89 ( linkage phases: ... 1 4 4 3 3 ... )
#>
#> ...-12-|-32-31-3-14-|-2-... OK
#>
#> ...-32-|-31-3-14-2-|-36-... OK
#>
#> ...-31-|-3-14-2-36-|-35-... OK
#>
#> ...-3-|-14-2-36-35-|-1-... OK
#>
#> ...-14-|-2-36-35-1-|-38-...
#> Alternative orders:
#> ... 14 2 36 35 1 38 ... : 0.00 ( linkage phases: ... 2 3 1 1 1 ... )
#> ... 14 2 36 1 35 38 ... : -1.27 ( linkage phases: ... 2 3 1 1 1 ... )
#>
#> ...-2-|-36-35-1-38-|-37-...
#> Alternative orders:
#> ... 2 38 1 35 36 37 ... : 0.00 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 35 1 38 37 ... : -0.99 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 38 35 1 36 37 ... : -1.15 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 35 1 38 36 37 ... : -1.16 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 38 35 1 37 ... : -1.18 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 38 1 35 37 ... : -1.45 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 35 38 1 37 ... : -2.21 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 2 36 1 35 38 37 ... : -2.25 ( linkage phases: ... 3 1 1 1 1 ... )
#>
#> ...-36-|-35-1-38-37-|-10-...
#> Alternative orders:
#> ... 36 37 1 35 38 10 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 35 1 38 10 ... : -0.25 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 35 1 38 37 10 ... : -1.23 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 38 35 1 10 ... : -1.35 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 38 35 1 37 10 ... : -1.42 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 38 1 35 10 ... : -1.54 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 38 1 35 37 10 ... : -1.69 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 35 38 1 37 10 ... : -2.45 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 1 35 38 37 10 ... : -2.50 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 1 38 35 10 ... : -2.78 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 36 37 35 38 1 10 ... : -2.85 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-35-|-1-38-37-10-|-28-...
#> Alternative orders:
#> ... 35 1 38 37 10 28 ... : 0.00 ( linkage phases: ... 1 1 1 1 4 ... )
#> ... 35 38 1 37 10 28 ... : -1.22 ( linkage phases: ... 1 1 1 1 4 ... )
#>
#> ...-1-|-38-37-10-28-|-26-... OK
#>
#> ...-38-|-37-10-28-26-|-42-... OK
#>
#> ...-37-|-10-28-26-42-|-6-... OK
#>
#> ...-10-|-28-26-42-6-|-39-...
#> Alternative orders:
#> ... 10 28 26 42 6 39 ... : 0.00 ( linkage phases: ... 4 2 3 1 1 ... )
#> ... 10 28 26 6 42 39 ... : -2.30 ( linkage phases: ... 4 2 3 1 1 ... )
#>
#> ...-28-|-26-42-6-39-|-41-...
#> Alternative orders:
#> ... 28 26 42 6 39 41 ... : 0.00 ( linkage phases: ... 2 3 1 1 1 ... )
#> ... 28 26 6 42 39 41 ... : -2.30 ( linkage phases: ... 2 3 1 1 1 ... )
#>
#> ...-26-|-42-6-39-41-|-40-...
#> Alternative orders:
#> ... 26 42 6 39 41 40 ... : 0.00 ( linkage phases: ... 3 1 1 1 1 ... )
#> ... 26 6 42 39 41 40 ... : -2.30 ( linkage phases: ... 3 1 1 1 1 ... )
#>
#> ...-42-|-6-39-41-40-|-15-... OK
#>
#> ...-6-|-39-41-40-15-|-5-... OK
#>
#> 39-|-41-40-15-5
#> Alternative orders:
#> ... 39 41 40 15 5 : 0.00 ( linkage phases: ... 1 1 1 1 )
#> ... 39 41 40 5 15 : -0.81 ( linkage phases: ... 1 1 1 1 )
# CHR2_frame <- mds_onemap(CHR_mks$sequences$CHR2)
# or
<- order_seq(CHR_mks$sequences$CHR2)
CHR2_ord <- make_seq(CHR2_ord, "force") CHR2_frame
rf_graph_table(CHR2_frame) # graphic not shown
As did before, we will not change the positions of the markers of this group.
<- CHR2_frame CHR2_final
# CHR2_frame <- mds_onemap(CHR_mks$sequences$CHR2)
# or
<- order_seq(CHR_mks$sequences$CHR2)
CHR2_ord <- make_seq(CHR2_ord, "force") CHR2_frame
rf_graph_table(CHR2_frame) # graphic not shown
As did before, we will not change the positions of the markers of this group.
<- CHR2_frame CHR2_final
# CHR3_frame <- mds_onemap(CHR_mks$sequences$CHR3)
# or
<- order_seq(CHR_mks$sequences$CHR3)
CHR3_ord <- make_seq(CHR3_ord, "force") CHR3_frame
rf_graph_table(CHR3_frame, mrk.axis = "numbers") # graphic not shown
Here, marker 29
has a color pattern too different from
the expected, removing it could be influential in other markers
ordering. Then we will remove them and search for a new order.
<- drop_marker(CHR3_frame, c(29))
CHR3_test_seq <- order_seq(CHR3_test_seq)
CHR3_test_ord <- make_seq(CHR3_test_ord, "force") CHR3_test_map
rf_graph_table(CHR3_test_map, mrk.axis = "numbers") #graphic not shown
Trying to add marker 29
again.
<- try_seq(CHR3_test_map, 29)
CHR3_add29_seq <- make_seq(CHR3_add29_seq, 12) # Marker 29 increase the map size disproportionately, it was removed from the map CHR3_add29
<- CHR3_test_map
CHR3_final rf_graph_table(CHR3_final, inter = FALSE)
Checking for better orders:
ripple_seq(CHR3_final)
#> 27-16-20-4-|-21-... OK
#>
#> ...-27-|-16-20-4-21-|-23-... OK
#>
#> ...-16-|-20-4-21-23-|-48-...
#> Alternative orders:
#> ... 16 20 4 21 23 48 ... : 0.00 ( linkage phases: ... 1 2 2 1 1 ... )
#> ... 16 20 4 23 21 48 ... : -2.27 ( linkage phases: ... 1 2 2 1 1 ... )
#>
#> ...-20-|-4-21-23-48-|-9-...
#> Alternative orders:
#> ... 20 4 21 23 48 9 ... : 0.00 ( linkage phases: ... 2 2 1 1 1 ... )
#> ... 20 4 23 21 48 9 ... : -2.27 ( linkage phases: ... 2 2 1 1 1 ... )
#>
#> ...-4-|-21-23-48-9-|-46-...
#> Alternative orders:
#> ... 4 21 23 48 9 46 ... : 0.00 ( linkage phases: ... 2 1 1 1 1 ... )
#> ... 4 23 21 48 9 46 ... : -2.27 ( linkage phases: ... 2 1 1 1 1 ... )
#>
#> ...-21-|-23-48-9-46-|-45-...
#> Alternative orders:
#> ... 21 23 48 9 46 45 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 21 23 48 46 9 45 ... : -1.00 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-23-|-48-9-46-45-|-47-...
#> Alternative orders:
#> ... 23 48 9 46 45 47 ... : 0.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 48 46 9 45 47 ... : -1.00 ( linkage phases: ... 1 1 1 1 1 ... )
#> ... 23 48 45 9 46 47 ... : -1.67 ( linkage phases: ... 1 1 1 1 1 ... )
#>
#> ...-48-|-9-46-45-47-|-50-...
#> Alternative orders:
#> ... 48 9 46 45 47 50 ... : 0.00 ( linkage phases: ... 1 1 1 1 2 ... )
#> ... 48 46 9 45 47 50 ... : -1.00 ( linkage phases: ... 1 1 1 1 2 ... )
#> ... 48 45 9 46 47 50 ... : -1.67 ( linkage phases: ... 1 1 1 1 2 ... )
#>
#> ...-9-|-46-45-47-50-|-24-... OK
#>
#> ...-46-|-45-47-50-24-|-49-... OK
#>
#> ...-45-|-47-50-24-49-|-52-... OK
#>
#> ...-47-|-50-24-49-52-|-51-...
#> Alternative orders:
#> ... 47 52 49 24 50 51 : 0.00 ( linkage phases: ... 2 1 1 1 1 )
#> ... 47 50 24 49 52 51 : -0.34 ( linkage phases: ... 2 1 1 1 1 )
#> ... 47 49 52 24 50 51 : -2.59 ( linkage phases: ... 2 1 1 1 1 )
#>
#> 50-|-24-49-52-51
#> Alternative orders:
#> ... 50 24 49 52 51 : 0.00 ( linkage phases: ... 1 1 1 1 )
#> ... 50 24 52 49 51 : -2.76 ( linkage phases: ... 1 1 1 1 )
Once all linkage groups were obtained using both strategies, we can
draw a map for each approach using the function draw_map
.
Since version 2.1.1007, OneMap
has a new version of
draw_map
, called draw_map2
. The new function
draws elegant linkage groups and presents new arguments to personalize
your draw.
If you prefer the old function, we also keep it. Follow examples on how to use both of them.
Drawing the map, which was built with only recombinations information.
<- list(LG1_final, LG2_final, LG3_final)
map1 draw_map(map1, names = TRUE, grid = TRUE, cex.mrk = 0.7)
Drawing the map, which built with reference genome and recombinations information
<- list(CHR1_final, CHR2_final, CHR3_final)
map2 draw_map(map2, names = TRUE, grid = TRUE, cex.mrk = 0.7)
We also can draw maps comparing corresponding linkage groups in each strategy:
<- list(LG1_final, CHR1_final)
CHR1_comp draw_map(CHR1_comp, names = TRUE, grid = TRUE, cex.mrk = 0.7)
<- list(LG3_final, CHR2_final)
CHR2_comp draw_map(CHR2_comp, names = TRUE, grid = TRUE, cex.mrk = 0.7)
Both strategies produced the same result for CHR2 (the map is only inverted).
<- list(LG2_final, CHR3_final)
CHR3_comp draw_map(CHR3_comp, names = TRUE, grid = TRUE, cex.mrk = 0.7)
Or groups alone:
draw_map(CHR1_final, names = TRUE, grid = TRUE, cex.mrk = 0.7)
Function draw_map
draws a straightforward graphic
representation of the genetic map. More recently, we developed a new
version called draw_map2
that brings a more sophisticated
figure. Furthermore, once the distances and the linkage phases are
estimated, other map figures can be drawn by the user with any
appropriate software. There are several free software that can be used,
such as MapChart
(Voorrips, 2002).
The same figures did with draw_map
can be done with the
draw_map2
function. But it has different capacities and
arguments. Here are some examples, but you can find more options on the
help page ?write_map2
.
Drawing the map, which was built with only recombinations information:
draw_map2(LG1_final, LG2_final, LG3_final, main = "Only with linkage information",
group.names = c("LG1", "LG2", "LG3"), output = "map.png")
NOTE: Check the GitHub vignette version to visualize the graphic.
The figure will be saved in your work directory with the default name
map.eps
. You can change the file name and extension
specifying them in the argument output
.
Drawing the map, which was built with reference genome and recombinations information
draw_map2(CHR1_final, CHR2_final, CHR3_final, output= "map_ref.pdf",
col.group = "#58A4B0",
col.mark= "#335C81")
NOTE: Check the GitHub vignette version to visualize the graphic.
With the argument tag
, we can highlight some markers
with other colors. The arguments col.group
,
col.mark
and col.tag
can be changed to
personalize the color of the groups, the markers, and the highlighted
markers, respectively.
We also can draw maps comparing corresponding linkage groups in each strategy:
draw_map2(LG1_final, CHR1_final, output = "map_comp.pdf", tag = c("M1","SNP2"))
NOTE: Check the GitHub vignette version to visualize the graphic.
When defining marker names in the tag
argument, all
markers with these names will be highlighted no matter in which group
it/they is/are.
draw_map2(LG2_final, CHR3_final, tag= c("SNP17", "SNP18", "M29"), main = "Chromosome 3",
group.names = c("Only linkage", "With genome"), centered = TRUE, output = "map_comp2.pdf")
NOTE: Check the GitHub vignette version to visualize the graphic.
If for any reason, one wants to estimate parameters for a given
linkage map (e.g., for other orders on published papers), it is
possible to define a sequence and use the map
function. For
example, for markers M30
, M12
,
M3
, M14
and M2
, in this order,
use:
<- make_seq(twopts, c(30, 12, 3, 14, 2))
any_seq <- map(any_seq))
(any_seq_map #>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 30 M30 0.00 a | | b a | | b
#> 12 M12 1.00 b | | a c | | a
#> 3 M3 20.57 o | | a o | | o
#> 14 M14 33.63 a | | o b | | o
#> 2 M2 41.70 o | | o o | | a
#>
#> 5 markers log-likelihood: -320.9012
NOTE: (new!) If
your sequence has many markers (more than 60), we suggest to speed up
map
using BatchMap parallelization approach. See section Speed up analysis with parallelization
for more information.
Warning: If you find an error message like:
Error in as_mapper(.f, ...) : argument ".f" is missing, with no default
It’s because the map
function has a very common name,
and you can have in your environment other functions with the same name.
In the case of the pointed error, R is using map
function
from purrr
package instead of OneMap
, to solve
this, simply specify that you want the OneMap
function with
::
command from stringr
package::
library(stringr)
<- onemap::map(any_seq)) (any_seq_map
This is a subset of the first linkage group. When used this way, the
map
function searches for the best combination of phases
between markers and prints the results.
(new!)
Warning: It is not our case in this example, but,
sometimes, it can happen that some markers in your sequence don’t reach
the OneMap
linkage criteria when linkage are estimated by
HMM multipoint approach using map
, it will produce an error
like this:
ERROR: The linkage between markers 1 and 2 did not reach the OneMap default criteria. They are probably segregating independently
You can evaluate the marker manually, or you can remove them
automatically using the map
argument
rm_unlinked = TRUE
. The map
function will
return a vector with marker numbers excluding the problematic marker;
then, you can repeat the process without the marker, using
make_seq
to create a new sequence and repeat the
map
. You can also do it automatically using the
map_avoid_unlinked
function:
<- map_avoid_unlinked(any_seq) LG2_test_map
Using this, if map
finds a problematic marker, it will
print a warning pointing to the marker number, which was removed, and
will automatically repeat the analysis without it.
Furthermore, a sequence can also have user-defined linkage phases. The next example shows (incorrect) phases used for the same order of markers:
<- make_seq(twopts, c(30, 12, 3, 14, 2), phase = c(4, 1, 4, 3))
any_seq <- map(any_seq)) (any_seq_map
If one needs to add or drop markers from a predefined sequence,
functions add_marker
and drop_marker
can be
used. For example, to add markers 4 to 8 to any_seq
.
<- add_marker(any_seq, 4:8))
(any_seq #>
#> Number of markers: 10
#> Markers in the sequence:
#> M30 M12 M3 M14 M2 M4 M5 M6 M7 M8
#>
#> Parameters not estimated.
Removing markers 3, 4, 5, 12, and 30 from any_seq
:
<- drop_marker(any_seq, c(3, 4, 5, 12, 30)))
(any_seq #>
#> Number of markers: 5
#> Markers in the sequence:
#> M14 M2 M6 M7 M8
#>
#> Parameters not estimated.
After that, the map needs to be re-estimated.
Warning: Only available for outcrossing and f2 intercross populations.
As already mentioned, OneMap
uses HMM multipoint
approach to estimate genetic distances, a very robust method, but it can
take time to run if you have many markers. In 2017, Schiffthaler et
al. released an OneMap
fork with modifications in CRAN and
in GitHub with the possibility of parallelizing the HMM chain dividing
markers in batches and use different cores for each phase. Their
approach speeds up our HMM and keeps the genetic distance estimation
quality. It allows us to divide the job into a maximum of four cores
according to the four possible phases for outcrossing mapping
populations. We add this parallelized approach to the functions:
map
, mds_onemap
, seriation
,
rcd
, record
and ug
. For better
efficiency, batches must be composed of 50 markers or more; therefore,
this approach is only recommended for linkage groups with many
markers.
The parallelization is here available for all types of operational
systems, however, we suggest setting argument
parallelization.type
to FORK
if you are not
using Windows system. It will improve the procedure speed.
Here we will show an example of how to use the BatchMap approach in some functions that requires HMM. For this, we simulated a group with 294 markers (we don’t want this vignette to take too much time to run, but usually maps with markers from high-throughput technologies result in larger groups). Before start, you can see the time spent on each approach (see also Session Info) in this example:
Without parallelization (h) | With parallelization (h) | |
---|---|---|
rcd | 0.6700558 | 0.1612458 |
record_map | 1.4368436 | 0.2907308 |
ug_map | 0.7145778 | 0.1884214 |
mds_onemap | 1.0643083 | 0.2827314 |
map | 2.0994486 | 0.6107456 |
Reading the simulated dataset:
<- read_onemap(system.file(package = "onemap", "extdata/simParall_out.raw")) # dataset available only in onemap github version
simParallel plot(simParallel, all=FALSE)
# Calculates two-points recombination fractions
<- rf_2pts(simParallel)
twopts <- make_seq(twopts, "all")
seq_all # There are no redundant markers
find_bins(simParallel)
# There are no distorted markers
<- plot(test_segregation(simParallel)) p
To prepare the data with a defined bach size, we use the function
pick_batch_sizes
. It selects a batch size that splits the
data into even groups. Argument size
defines the batch size
next to which an optimum size will be searched. overlap
defines the number of markers that overlap between the present batch and
next. This is used because pre-defined phases at these overlap markers
in the present batch are used to start the HMM in the next batch. The
around
argument defines how much the function can vary
around the defined number in size
to search for the optimum
batch size.
Some aspects should be considered to define these arguments because if the batch size were set too high, there would be less gain in execution time. If the overlap size were too small, phases would be incorrectly estimated, and large gaps would appear in the map, inflating its size. In practice, these values will depend on many factors such as population size, marker quality, and species. BatchMap authors recommended trying several configurations on a subset of data and select the best performing one.
<- pick_batch_sizes(input.seq = seq_all,
batch_size size = 80,
overlap = 30,
around = 10)
batch_size
To use parallelized approach you just need to include the arguments when using the functions:
# Without parallelization
<- rcd(input.seq = seq_all)
rcd_map # With parallelization
<- rcd(input.seq = seq_all,
rcd_map_par phase_cores = 4,
size = batch_size,
overlap = 30)
# Without parallelization
<- record(input.seq = seq_all)
record_map # With parallelization
<- record(input.seq = seq_all,
record_map_par phase_cores = 4,
size = batch_size,
overlap = 30)
# Without parallelization
<- ug(input.seq = seq_all)
ug_map # With parallelization
<- ug(input.seq = seq_all,
ug_map_par phase_cores = 4,
size = batch_size,
overlap = 30)
# Without parallelization ok
<- mds_onemap(input.seq = seq_all)
map_mds # With parallelization
<- mds_onemap(input.seq = seq_all,
map_mds_par phase_cores = 4,
size = batch_size,
overlap = 30)
Because we simulate this dataset, we know the correct order. We can
use map_overlapping_batches
to estimate genetic distance in
this case. This is equivalent to map
, but with a
parallelized process.
Similarly, with map
, using argument
rm_unlinked = TRUE
the function will return a vector with
marker numbers without the problematic marker. To repeat the analysis
removing automatically all problematic markers use
map_avoid_unlinked
:
# Without parallelization
<- map_avoid_unlinked(input.seq = seq_all)
batch_map # With parallelization
<- map_avoid_unlinked(input.seq = seq_all,
batch_map_par size = batch_size,
phase_cores = 4,
overlap = 30)
As you can see in the above maps, heuristic ordering algorithms do not return an optimal order result, mostly if you don’t have many individuals in your population. Because of the erroneous order, generated map size is not close to the simulated size (100 cM) and their heatmaps don’t present the expected color pattern. Two of them get close to the color pattern, they are the ug and the MDS method. They present the right global ordering but not local. If you have a reference genome, you can use its position information to rearrange the local order.
In the older version, users could only access the estimated linkage phase by observing the print in the console:
CHR3_final#>
#> Printing map:
#>
#> Markers Position Parent 1 Parent 2
#>
#> 27 M27 0.00 b | | o a | | a
#> 16 M16 11.76 a | | a b | | o
#> 20 M20 22.94 a | | b c | | d
#> 4 M4 35.71 a | | o o | | b
#> 21 M21 49.41 o | | o a | | b
#> 23 M23 55.24 a | | o a | | o
#> 48 SNP21 67.50 a | | b a | | a
#> 9 M9 74.02 a | | b a | | a
#> 46 SNP18 80.43 a | | b a | | a
#> 45 SNP17 95.40 a | | b a | | a
#> 47 SNP20 119.95 a | | b a | | a
#> 50 SNP23 158.43 a | | b b | | a
#> 24 M24 163.43 a | | b b | | a
#> 49 SNP22 166.72 a | | b b | | a
#> 52 SNP25 171.62 a | | b b | | a
#> 51 SNP24 182.56 a | | b b | | a
#>
#> 16 markers log-likelihood: -811.5279
Now, you can export this information into a data.frame using:
<- parents_haplotypes(CHR3_final))
(parents_haplot #> group mk.number mk.names dist P1_1 P1_2 P2_1 P2_2
#> 1 Group - 1 27 M27 0.00000 b o a a
#> 2 Group - 1 16 M16 11.75801 a a b o
#> 3 Group - 1 20 M20 22.93778 a b c d
#> 4 Group - 1 4 M4 35.70756 a o o b
#> 5 Group - 1 21 M21 49.40658 o o a b
#> 6 Group - 1 23 M23 55.24441 a o a o
#> 7 Group - 1 48 SNP21 67.49802 a b a a
#> 8 Group - 1 9 M9 74.02260 a b a a
#> 9 Group - 1 46 SNP18 80.43325 a b a a
#> 10 Group - 1 45 SNP17 95.40132 a b a a
#> 11 Group - 1 47 SNP20 119.95322 a b a a
#> 12 Group - 1 50 SNP23 158.42672 a b b a
#> 13 Group - 1 24 M24 163.43246 a b b a
#> 14 Group - 1 49 SNP22 166.72021 a b b a
#> 15 Group - 1 52 SNP25 171.61828 a b b a
#> 16 Group - 1 51 SNP24 182.55630 a b b a
write.table(parents_haplot, "parents_haplot.txt")
The data.frame contains: group ID (group), marker number (mk.number) and names (mk.names), position in centimorgan (dist) and parents haplotypes (P1_1, P1_2, P2_1, P2_2).
You can also obtain a data.frame with a list of sequences and personalize the group names:
parents_haplotypes(CHR2_final,CHR3_final, group_names=c("CHR2","CHR3"))
#> group mk.number mk.names dist P1_1 P1_2 P2_1 P2_2
#> 1 CHR2 43 SNP14 0.00000 a b a a
#> 2 CHR2 22 M22 13.50513 a b a a
#> 3 CHR2 7 M7 23.08429 a a a b
#> 4 CHR2 18 M18 64.51343 a b c d
#> 5 CHR2 8 M8 70.19895 b a b a
#> 6 CHR2 13 M13 72.86392 a o a o
#> 7 CHR2 44 SNP16 78.74250 b a b a
#> 8 CHR3 27 M27 0.00000 b o a a
#> 9 CHR3 16 M16 11.75801 a a b o
#> 10 CHR3 20 M20 22.93778 a b c d
#> 11 CHR3 4 M4 35.70756 a o o b
#> 12 CHR3 21 M21 49.40658 o o a b
#> 13 CHR3 23 M23 55.24441 a o a o
#> 14 CHR3 48 SNP21 67.49802 a b a a
#> 15 CHR3 9 M9 74.02260 a b a a
#> 16 CHR3 46 SNP18 80.43325 a b a a
#> 17 CHR3 45 SNP17 95.40132 a b a a
#> 18 CHR3 47 SNP20 119.95322 a b a a
#> 19 CHR3 50 SNP23 158.42672 a b b a
#> 20 CHR3 24 M24 163.43246 a b b a
#> 21 CHR3 49 SNP22 166.72021 a b b a
#> 22 CHR3 52 SNP25 171.61828 a b b a
#> 23 CHR3 51 SNP24 182.55630 a b b a
Function progeny_haplotypes
generates a data.frame with
progeny phased haplotypes estimated by OneMap
HMM. For
progeny, the HMM results in probabilities for each possible genotype,
then the generated data.frame contains all possible genotypes. If
most_likely = TRUE
, the most likely genotype receives 1 and
the rest 0 (if there are two most likely both receive 0.5), if
most_likely = FALSE
genotypes probabilities will be
according to the HMM results. You can choose which individual to be
evaluated in ind
. The data.frame is composed by the
information: individual (ind) and group (grp) ID, position in
centimorgan (pos), progeny homologs (homologs), and from each parent the
allele came (parents).
<- progeny_haplotypes(CHR2_final, most_likely = TRUE, ind = c(1,2), group_names = "CHR2_final"))
(progeny_haplot #> ind marker grp pos prob parents parents.homologs allele
#> 1 IND1 SNP14 CHR2_final 0.00000 0 P1 H1 P1_H1
#> 2 IND1 M22 CHR2_final 13.50513 0 P1 H1 P1_H1
#> 3 IND1 M7 CHR2_final 23.08429 0 P1 H1 P1_H1
#> 4 IND1 M18 CHR2_final 64.51343 1 P1 H1 P1_H1
#> 5 IND1 M8 CHR2_final 70.19895 1 P1 H1 P1_H1
#> 6 IND1 M13 CHR2_final 72.86392 1 P1 H1 P1_H1
#> 7 IND1 SNP16 CHR2_final 78.74250 1 P1 H1 P1_H1
#> 8 IND2 SNP14 CHR2_final 0.00000 1 P1 H1 P1_H1
#> 9 IND2 M22 CHR2_final 13.50513 1 P1 H1 P1_H1
#> 10 IND2 M7 CHR2_final 23.08429 1 P1 H1 P1_H1
#> 11 IND2 M18 CHR2_final 64.51343 1 P1 H1 P1_H1
#> 12 IND2 M8 CHR2_final 70.19895 1 P1 H1 P1_H1
#> 13 IND2 M13 CHR2_final 72.86392 1 P1 H1 P1_H1
#> 14 IND2 SNP16 CHR2_final 78.74250 1 P1 H1 P1_H1
#> 15 IND1 SNP14 CHR2_final 0.00000 1 P1 H2 P1_H2
#> 16 IND1 M22 CHR2_final 13.50513 1 P1 H2 P1_H2
#> 17 IND1 M7 CHR2_final 23.08429 1 P1 H2 P1_H2
#> 18 IND1 M18 CHR2_final 64.51343 0 P1 H2 P1_H2
#> 19 IND1 M8 CHR2_final 70.19895 0 P1 H2 P1_H2
#> 20 IND1 M13 CHR2_final 72.86392 0 P1 H2 P1_H2
#> 21 IND1 SNP16 CHR2_final 78.74250 0 P1 H2 P1_H2
#> 22 IND2 SNP14 CHR2_final 0.00000 0 P1 H2 P1_H2
#> 23 IND2 M22 CHR2_final 13.50513 0 P1 H2 P1_H2
#> 24 IND2 M7 CHR2_final 23.08429 0 P1 H2 P1_H2
#> 25 IND2 M18 CHR2_final 64.51343 0 P1 H2 P1_H2
#> 26 IND2 M8 CHR2_final 70.19895 0 P1 H2 P1_H2
#> 27 IND2 M13 CHR2_final 72.86392 0 P1 H2 P1_H2
#> 28 IND2 SNP16 CHR2_final 78.74250 0 P1 H2 P1_H2
#> 29 IND1 SNP14 CHR2_final 0.00000 0 P2 H1 P2_H1
#> 30 IND1 M22 CHR2_final 13.50513 0 P2 H1 P2_H1
#> 31 IND1 M7 CHR2_final 23.08429 0 P2 H1 P2_H1
#> 32 IND1 M18 CHR2_final 64.51343 0 P2 H1 P2_H1
#> 33 IND1 M8 CHR2_final 70.19895 0 P2 H1 P2_H1
#> 34 IND1 M13 CHR2_final 72.86392 0 P2 H1 P2_H1
#> 35 IND1 SNP16 CHR2_final 78.74250 0 P2 H1 P2_H1
#> 36 IND2 SNP14 CHR2_final 0.00000 0 P2 H1 P2_H1
#> 37 IND2 M22 CHR2_final 13.50513 0 P2 H1 P2_H1
#> 38 IND2 M7 CHR2_final 23.08429 0 P2 H1 P2_H1
#> 39 IND2 M18 CHR2_final 64.51343 0 P2 H1 P2_H1
#> 40 IND2 M8 CHR2_final 70.19895 0 P2 H1 P2_H1
#> 41 IND2 M13 CHR2_final 72.86392 0 P2 H1 P2_H1
#> 42 IND2 SNP16 CHR2_final 78.74250 0 P2 H1 P2_H1
#> 43 IND1 SNP14 CHR2_final 0.00000 1 P2 H2 P2_H2
#> 44 IND1 M22 CHR2_final 13.50513 1 P2 H2 P2_H2
#> 45 IND1 M7 CHR2_final 23.08429 1 P2 H2 P2_H2
#> 46 IND1 M18 CHR2_final 64.51343 1 P2 H2 P2_H2
#> 47 IND1 M8 CHR2_final 70.19895 1 P2 H2 P2_H2
#> 48 IND1 M13 CHR2_final 72.86392 1 P2 H2 P2_H2
#> 49 IND1 SNP16 CHR2_final 78.74250 1 P2 H2 P2_H2
#> 50 IND2 SNP14 CHR2_final 0.00000 1 P2 H2 P2_H2
#> 51 IND2 M22 CHR2_final 13.50513 1 P2 H2 P2_H2
#> 52 IND2 M7 CHR2_final 23.08429 1 P2 H2 P2_H2
#> 53 IND2 M18 CHR2_final 64.51343 1 P2 H2 P2_H2
#> 54 IND2 M8 CHR2_final 70.19895 1 P2 H2 P2_H2
#> 55 IND2 M13 CHR2_final 72.86392 1 P2 H2 P2_H2
#> 56 IND2 SNP16 CHR2_final 78.74250 1 P2 H2 P2_H2
You can also have a view of progeny estimated haplotypes using
plot
. It shows which markers came from each parent’s
homologs. position
argument defines if haplotypes will be
plotted by homologs (stack
) or alleles
(split
). split
option is a good way to view
the likelihoods of each allele.
plot(progeny_haplot, position = "stack")
plot(progeny_haplot, position = "split")
sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United States.1252
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] rmarkdown_2.18 knitr_1.40 stringr_1.4.1 onemap_3.0.0
#>
#> loaded via a namespace (and not attached):
#> [1] minqa_1.2.4 colorspace_2.0-3 ggsignif_0.6.4
#> [4] deldir_1.0-6 ellipsis_0.3.2 class_7.3-20
#> [7] htmlTable_2.4.1 base64enc_0.1-3 rstudioapi_0.14
#> [10] proxy_0.4-27 mice_3.14.0 farver_2.1.1
#> [13] ggpubr_0.4.0 fansi_1.0.3 codetools_0.2-18
#> [16] splines_4.1.0 doParallel_1.0.17 cachem_1.0.6
#> [19] Formula_1.2-4 polynom_1.4-1 jsonlite_1.8.0
#> [22] nloptr_2.0.3 broom_1.0.1 cluster_2.1.3
#> [25] png_0.1-7 compiler_4.1.0 httr_1.4.4
#> [28] backports_1.4.1 assertthat_0.2.1 Matrix_1.4-1
#> [31] fastmap_1.1.0 lazyeval_0.2.2 cli_3.4.1
#> [34] htmltools_0.5.2 tools_4.1.0 gtable_0.3.1
#> [37] glue_1.4.2 rebus.base_0.0-3 reshape2_1.4.4
#> [40] dplyr_1.0.9 Rcpp_1.0.8.3 carData_3.0-5
#> [43] jquerylib_0.1.4 vctrs_0.5.1 ape_5.6-2
#> [46] gdata_2.18.0.1 nlme_3.1-157 iterators_1.0.14
#> [49] pinfsc50_1.2.0 xfun_0.31 rebus.datetimes_0.0-2
#> [52] lme4_1.1-29 lifecycle_1.0.3 rebus.numbers_0.0-1
#> [55] weights_1.0.4 gtools_3.9.2.2 rstatix_0.7.1
#> [58] dendextend_1.16.0 princurve_2.1.6 candisc_0.8-6
#> [61] MASS_7.3-57 scales_1.2.1 heplots_1.4-2
#> [64] parallel_4.1.0 smacof_2.1-5 RColorBrewer_1.1-3
#> [67] yaml_2.3.5 gridExtra_2.3 ggplot2_3.4.0
#> [70] sass_0.4.0 rpart_4.1.16 latticeExtra_0.6-30
#> [73] stringi_1.7.6 highr_0.9 foreach_1.5.2
#> [76] plotrix_3.8-2 permute_0.9-7 e1071_1.7-11
#> [79] checkmate_2.1.0 boot_1.3-28 rlang_1.0.6
#> [82] pkgconfig_2.0.3 evaluate_0.18 lattice_0.20-45
#> [85] purrr_0.3.4 labeling_0.4.2 htmlwidgets_1.5.4
#> [88] tidyselect_1.2.0 plyr_1.8.7 magrittr_2.0.1
#> [91] R6_2.5.1 generics_0.1.3 nnls_1.4
#> [94] Hmisc_4.7-0 DBI_1.1.3 mgcv_1.8-40
#> [97] pillar_1.8.1 foreign_0.8-82 withr_2.5.0
#> [100] rebus_0.1-3 survival_3.3-1 abind_1.4-5
#> [103] nnet_7.3-17 rebus.unicode_0.0-2 tibble_3.1.7
#> [106] car_3.1-1 interp_1.1-3 wordcloud_2.6
#> [109] utf8_1.2.2 ellipse_0.4.3 plotly_4.10.1
#> [112] vcfR_1.12.0 viridis_0.6.2 jpeg_0.1-9
#> [115] grid_4.1.0 data.table_1.14.2 vegan_2.6-2
#> [118] digest_0.6.29 tidyr_1.2.0 munsell_0.5.0
#> [121] viridisLite_0.4.1 bslib_0.4.1
Buetow, K. H., Chakravarti, A. Multipoint gene mapping using seriation. I. General methods. American Journal of Human Genetics 41, 180-188, 1987.
Doerge, R.W. Constructing genetic maps by rapid chain delineation. Journal of Agricultural Genomics 2, 1996.
Mollinari, M., Margarido, G. R. A., Vencovsky, R. and Garcia, A. A. F. Evaluation of algorithms used to order markers on genetics maps. Heredity 103, 494-502, 2009.
Schiffthaler, B., Bernhardsson, C., Ingvarsson, P. K., & Street, N. R. BatchMap: A parallel implementation of the OneMap R package for fast computation of F1 linkage maps in outcrossing species. PLoS ONE, 12(12), 1–12, 2017.
Tan, Y., Fu, Y. A novel method for estimating linkage maps. Genetics 173, 2383-2390, 2006.
Van Os H, Stam P, Visser R.G.F., Van Eck H.J. RECORD: a novel method for ordering loci on a genetic linkage map. Theor Appl Genet 112, 30-40, 2005.
Voorrips, R.E. MapChart: software for the graphical presentation of linkage maps and QTLs. Journal of Heredity 93, 77-78, 2002.
Voorrips, R. E., Maliepaard, C. A. The simulation of meiosis in diploid and tetraploid organisms using various genetic models. BMC Bioinformatics, 13(1), 248, 2012.
Wu, R., Ma, C.X., Painter, I. and Zeng, Z.-B. Simultaneous maximum likelihood estimation of linkage and linkage phases in outcrossing species. Theoretical Population Biology 61, 349-363, 2002a.
Wu, R., Ma, C.-X., Wu, S. S. and Zeng, Z.-B. Linkage mapping of sex-specific differences. Genetical Research 79, 85-96, 2002b.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.