Mapping and QTL Example Tools for Polyploid RGC11
Jeekin Lau
2023-03-17
Prerequisites
For this tutorial, you will need to install mappoly
polymapr qtlpoly viewpoly. If you
want to download the stable version and the versions please run the
following lines.
install.packages(c("mappoly","polymapR","qtlpoly","viewpoly"))If you want to run the more experimental versions of the software, run following code:
install.packages("devtools")
devtools::install_github("mmollina/MAPpoly")
devtools::install_github("gabrielgesteira/QTLpoly")
devtools::install_github("mmollina/viewpoly")Introduction
This is a tutorial using a very clean dataset from a SNP array on a well characterized autotetraploid garden rose mapping population that has been published already Rose Rosette Disease Resistance Loci Detected in Two Interconnected Tetraploid Garden Rose Populations
To download the dataset please go to: https://github.com/jeekinlau/mapping_example/tree/main/data
Tutorial link:
More in depth Tutorials from the developers
mappoly
qtlpoly
viewpoly
polymapR
Import the data
Download the dataset and import the dataset into R. The data is
organized in a way that mappoly recognizes it. must be in
this format if you want to use mappoly::read_geno_csv() to
read in your data. If you have a vcf file you can also use
mappoly::read_vcf().
Normally in the
mappoly::read_geno_csv(filter.non.conforming=T) but in this
case below we set toFalse in order to do some preliminary
QC done in polymapR.
CSV input file
setwd("~/GitHub/mapping_example/tutorial")
library(mappoly)## =====================================================
## MAPpoly Package [Version 0.3.3]
## More information: https://github.com/mmollina/MAPpoly
## =====================================================dat = read_geno_csv(file="../data/BExMG_subset_with_contaminants_2075mrks.csv", ploidy = 4, filter.non.conforming = F)## Reading the following data:
## Ploidy level: 4
## No. individuals: 169
## No. markers: 2075
## No. informative markers: 2075 (100%)
## ...
## Done with reading.Importing VCFs
This is not part of the tutorial. mappoly is able to take VCF files and this is just showing how to import data.
dat_vcf = read_vcf("../data/sweetpotato_chr1.vcf.gz", parent.1 = "PARENT1", parent.2 = "PARENT2",ploidy = 6)QC
Here we will look for contaminants and duplicates and remove them
from the analysis. Good idea because you have no idea what kind of
errors you have in genotyping. We will use the functions the developers
have implemented to help going between different software. What we will
do here is use polymapR to find the contaminants and
duplicate individuals and remove them from our analysis. We will produce
a new genetic input file.
library(polymapR)##
## Attaching package: 'polymapR'## The following object is masked from 'package:mappoly':
##
## compare_mapspolymapR_dosage_matrix=export_data_to_polymapR(dat)
PCA_progeny(polymapR_dosage_matrix)
contaminants = c("X16009_N001", "X16009_N010","SW","X16009_N005","X16009_N006","X16009_N009")
screened_data <- screen_for_duplicate_individuals(dosage_matrix = ALL_dosages,
cutoff = 0.95,
plot_cor = T)##
## Combining F1_060 & F1_081 into F1_060
##
## Combining F1_046 & F1_085 into F1_046## Warning in screen_for_duplicate_individuals(dosage_matrix = ALL_dosages, : Multiple duplicates of
## single genotype identified at this threshold. Attempting to merge...
##
## Combining F1_032 & F1_100 into F1_032
##
## Combining F1_032 & F1_114 into F1_032
##
## Combining F1_079 & F1_101 into F1_079
##
## Combining F1_106 & F1_113 into F1_106
##
## Combining F1_112 & F1_199 into F1_112
##
## #### 7 individuals removed:
##
## _ _ _ _
## [1,] "F1_081" "F1_085" "F1_100" "F1_114"
## [2,] "F1_101" "F1_113" "F1_199" ""
## |_ |_ |_ |_ |
## |:------|:------|:------|:------|
## |F1_081 |F1_085 |F1_100 |F1_114 |
## |F1_101 |F1_113 |F1_199 | |duplicates = c("X16035_N028.1", "X16035_N029.1","X16400_N023","X16400_N038","X16405_N113","X16400_N006")
remove = c(contaminants,duplicates)
genotype_data = read.csv("../data/BExMG_subset_with_contaminants_2075mrks.csv")
genotype_data_after_QC=genotype_data[,which(!colnames(genotype_data)%in%remove)]
write.csv(genotype_data_after_QC,"../data/BExMG_subset_with_contaminants_2075mrks_afterQC.csv", row.names = F)Read back in the new fixed genotype file
setwd("~/GitHub/mapping_example/tutorial")
library(mappoly)
dat = read_geno_csv(file="../data/BExMG_subset_with_contaminants_2075mrks_afterQC.csv", ploidy = 4, filter.non.conforming = T)## Reading the following data:
## Ploidy level: 4
## No. individuals: 157
## No. markers: 2075
## No. informative markers: 2075 (100%)
## ...
## Done with reading.
## Filtering non-conforming markers.
## ...
## Done with filtering.Filtering Data
Filter by missing data
Typically, there will be many markers filtered out when filtered by marker, or by individual. Because this is a very well curated dataset, There are no markers filtered out. Start with 0.05 or 0.10 but there are no “rules” for what level to use.
## Filtering dataset by marker
dat = filter_missing(input.data = dat, type = "marker",
filter.thres = 0.10, inter = F)
## Filtering dataset by individual
dat = filter_missing(input.data = dat, type = "individual",
filter.thres = 0.10, inter = F)
dat## This is an object of class 'mappoly.data'
## Ploidy level: 4
## No. individuals: 157
## No. markers: 2074
## Missing data: 1.35%
## Redundant markers: 0.05%
##
## This dataset contains chromosome information.
## ----------
## No. of markers per dosage combination in both parents:
## P1 P2 freq
## 0 1 182
## 0 2 52
## 0 3 15
## 1 0 168
## 1 1 215
## 1 2 102
## 1 3 52
## 1 4 12
## 2 0 48
## 2 1 113
## 2 2 116
## 2 3 108
## 2 4 52
## 3 0 11
## 3 1 50
## 3 2 115
## 3 3 223
## 3 4 203
## 4 1 5
## 4 2 52
## 4 3 180Filter by segregation distortion
Also very few markers removed due to good dataset. The figure below shows the marker types, the dosage of all the marker types, and the segregation distortion of the markers. Here also we apply an alpha=0.05 with Bonfferoni correction.
## Filter according to segregation test
pval.bonf = 0.05/dat$n.mrk
mrks.chi.filt = filter_segregation(dat, chisq.pval.thres = pval.bonf, inter = F)
seq.init = make_seq_mappoly(mrks.chi.filt)
seq.init## This is an object of class 'mappoly.sequence'
## ------------------------
## Parameters not estimated
## ------------------------
## Ploidy level: 4
## No. individuals: 157
## No. markers: 2062
##
## ----------
## No. markers per sequence:
## chrom No.mrk
## 1 267
## 2 432
## 3 234
## 4 263
## 5 234
## 6 325
## 7 307
##
## ----------
## No. of markers per dosage in both parents:
## dP1 dP2 freq
## 0 1 182
## 0 2 52
## 0 3 14
## 1 0 168
## 1 1 213
## 1 2 102
## 1 3 51
## 1 4 12
## 2 0 48
## 2 1 113
## 2 2 116
## 2 3 108
## 2 4 52
## 3 0 11
## 3 1 46
## 3 2 115
## 3 3 221
## 3 4 202
## 4 1 4
## 4 2 52
## 4 3 180plot(seq.init)
Calculating recombination fractions
This step allows for the pairwise recombination fraction calculation
between all markers. If you computer has multiple cores use the
parallel::detectCores() function to find out how many cores
your local computer has. This step takes a while especially with a
laptop with limited cores and RAM memory. (there are tricks
around this email me for alternative strategies around
this)
## Accounting for recombination fractions
all.rf.pairwise = est_pairwise_rf(input.seq = seq.init, ncpus = 6)## INFO: Using 6 CPUs for calculation.
## INFO: Done with 2124891 pairs of markers
## INFO: Calculation took: 226.69 secondsPlotting RF matrix
Plotting the RF matrix shows a heat map of markers that are closely linked versus those that are far away. We see that those markers group into the 7 chromosomes of the Rose genome.
mat = rf_list_to_matrix(all.rf.pairwise)## INFO: Going singlemode. Using one CPU.plot(mat)
Grouping
Next step is to group the markers into a list containing the necessary objects needed for mapping (ie markers and their two point rf estimate and a sequence of markers). There are several ways to approach this. Below we compare the grouping using UPGMA method, genomic (physical position according to genome position), and an intersect between the two.
## Getting group information (UPGMA)
grs = group_mappoly(input.mat = mat,
expected.groups = 7,
comp.mat = TRUE,
inter = F)
grs## This is an object of class 'mappoly.group'
## ------------------------------------------
## Criteria used to assign markers to groups:
##
## - Number of markers: 2062
## - Number of linkage groups: 7
## - Number of markers per linkage groups:
## group n.mrk
## 1 275
## 2 448
## 3 229
## 4 292
## 5 263
## 6 236
## 7 319
## ------------------------------------------
## 1 2 3 7 4 5 6 NoChr
## 1 266 0 6 0 1 0 2 0
## 2 1 431 1 6 1 1 7 0
## 3 0 1 226 2 0 0 0 0
## 4 0 0 1 291 0 0 0 0
## 5 0 0 0 1 261 0 1 0
## 6 0 0 0 4 0 232 0 0
## 7 0 0 0 3 0 1 315 0
## ------------------------------------------plot(grs)
## Making groups with intersection between UPGMA (rf) + genomic information (discards scaffolds)
LGS.inter=vector("list", 7)
for(j in 1:7){
temp1 = make_seq_mappoly(grs, j, genomic.info=1)
tpt = make_pairs_mappoly(all.rf.pairwise, input.seq = temp1)
temp2 = rf_snp_filter(input.twopt = tpt, diagnostic.plot = FALSE)
tpt2 = make_pairs_mappoly(tpt, input.seq = temp2)
LGS.inter[[as.numeric(names(table(temp2$chrom))[which.max(table(temp2$chrom))])]] = list(seq = temp2, tpt = tpt2)
}
## Making groups only with genomic information (tends to bring more markers, but also discards scaffolds)
LGS.genomic = vector("list", 7)
for (i in 1:7){
tempseq1 = make_seq_mappoly(dat, arg = paste0("seq",i), genomic.info = 1)
mrks = intersect(tempseq1$seq.mrk.names, seq.init$seq.mrk.names)
tempseq = make_seq_mappoly(dat, arg = mrks)
temptpt = make_pairs_mappoly(all.rf.pairwise, input.seq = tempseq)
rffilt = rf_snp_filter(input.twopt = temptpt, diagnostic.plot = FALSE)
temptpt2 = make_pairs_mappoly(temptpt, input.seq = rffilt)
LGS.genomic[[as.numeric(unique(rffilt$chrom))]] = list(seq = rffilt, tpt = temptpt2)
}
## Making groups with UPGMA (rf) information (includes scaffolds)
LGS.upgma=vector("list", 7)
for(j in 1:7){
temp1 = make_seq_mappoly(grs, j)
tpt = make_pairs_mappoly(all.rf.pairwise, input.seq = temp1)
temp2 = rf_snp_filter(input.twopt = tpt, diagnostic.plot = FALSE)
tpt2 = make_pairs_mappoly(tpt, input.seq = temp2)
LGS.upgma[[as.numeric(names(table(temp2$chrom))[which.max(table(temp2$chrom))])]] = list(seq = temp2, tpt = tpt2)
}
## Comparing number of markers in each group
comp = data.frame(UPGMA_Genomic = unlist(lapply(LGS.inter, function(x) length(x$seq$seq.num))),
Genomic = unlist(lapply(LGS.genomic, function(x) length(x$seq$seq.num))),
UPGMA = unlist(lapply(LGS.upgma, function(x) length(x$seq$seq.num))))
comp## UPGMA_Genomic Genomic UPGMA
## 1 253 253 264
## 2 410 412 427
## 3 214 222 218
## 4 250 252 250
## 5 223 222 224
## 6 299 308 304
## 7 279 294 279We see that all three are very similar and for the purpose of this example lets just take the upgma grouping and run with that.
LGS<-LGS.interHMM chain
The HMM chain is used to order the markers. Below is a example of just ordering one chromosome. Later we will show how to parallelize this task by chromosome. (NOTE the map is only running quickly here because of the small number of markers and the good quality of the data)
single_chrom <- est_rf_hmm_sequential(input.seq = LGS[[1]]$seq,
start.set = 3,
thres.twopt = 10,
thres.hmm = 50,
extend.tail = 30,
twopt = LGS[[1]]$tpt,
verbose = TRUE,
phase.number.limit = 20,
sub.map.size.diff.limit = 5)## Number of markers: 253## ════════════════════════════════════════════════════════════════════════════════ Initial sequence ══## 3 markers...
## • Trying sequence: 1 2 3 :
## 4 phase(s): . . . .
## • Trying sequence: 2 3 4 :
## 3 phase(s): . . .## ══════════════════════════════════════════════════════════════════════ Done with initial sequence ══## Making 'reestimate.single.ph.configuration = TRUE' to use map expansion## 4 /5 :(2%) 6 : 2 ph (1/2) -- tail: 3 •||| •••• ||•| ••••
## 5 /6 :(2.4%) 7 : 1 ph (1/1) -- tail: 4 ••|• •|||
## 6 /7 :(2.8%) 9 : 1 ph (1/1) -- tail: 5 |||| |||•
## 7 /8 :(3.2%) 10 : 1 ph (1/1) -- tail: 6 |||| |||•
## 8 /9 :(3.6%) 11 : 1 ph (1/1) -- tail: 7 |||| ||••
## 9 /10 :(4%) 12 : 1 ph (1/1) -- tail: 8 |||| |||•
## 10 /11 :(4.3%) 14 : 4 ph (1/4) -- tail: 9 ••|| ••|| ... ••|| •|•|
## 11 /12 :(4.7%) 16 : 1 ph (1/1) -- tail: 10 •••• •••|
## 12 /13 :(5.1%) 17 : 1 ph (1/1) -- tail: 11 •••• ••||
## 13 /14 :(5.5%) 18 : 1 ph (1/1) -- tail: 12 ••|• •••|
## 14 /15 :(5.9%) 20 : 4 ph (1/4) -- tail: 13 •••| •••• ... ••|• ••••
## 15 /16 :(6.3%) 21 : 1 ph (1/1) -- tail: 14 ••|• •••|
## 16 /17 :(6.7%) 22 : 1 ph (1/1) -- tail: 15 ••|• •••|
## 17 /18 :(7.1%) 23 : 2 ph (1/2) -- tail: 16 |||| ||•| |||| |||•
## 18 /19 :(7.5%) 24 : 1 ph (1/1) -- tail: 17 ••|• •••|
## 19 /20 :(7.9%) 26 : 1 ph (1/1) -- tail: 18 •••• •••|
## 20 /21 :(8.3%) 28 : 2 ph (1/2) -- tail: 19 •••• •|•| •••• |••|
## 21 /22 :(8.7%) 29 : 1 ph (1/1) -- tail: 20 •••• |•••
## 22 /23 :(9.1%) 30 : 4 ph (1/4) -- tail: 21 •••| |••| ... ••|• |••|
## 23 /24 :(9.5%) 31 : 2 ph (1/2) -- tail: 22 •••| |•|• |••• |•|•
## 24 /25 :(9.9%) 32 : 1 ph (1/1) -- tail: 23 |||• ||••
## 25 /26 :(10.3%) 34 : 1 ph (1/1) -- tail: 24 |||• |||•
## 26 /27 :(10.7%) 35 : 1 ph (1/1) -- tail: 25 •••• •••|
## 27 /28 :(11.1%) 36 : 1 ph (1/1) -- tail: 26 |||• |||•
## 28 /29 :(11.5%) 37 : 1 ph (1/1) -- tail: 27 •••| •••|
## 29 /30 :(11.9%) 38 : 2 ph (1/2) -- tail: 28 |||| |•|| |||| ||•|
## 30 /31 :(12.3%) 39 : 1 ph (1/1) -- tail: 29 |||• |||•
## 31 /32 :(12.6%) 40 : 1 ph (1/1) -- tail: 30 |||| |••|
## 32 /33 :(13%) 41 : 1 ph (1/1) -- tail: 30 •••• •||•
## 33 /34 :(13.4%) 42 : 1 ph (1/1) -- tail: 30 |••• |•••
## 34 /35 :(13.8%) 44 : 1 ph (1/1) -- tail: 30 •••• •||•
## 35 /36 :(14.2%) 45 : 1 ph (1/1) -- tail: 30 ||•| ||||
## 36 /37 :(14.6%) 46 : 1 ph (1/1) -- tail: 30 ••|| •••|
## 37 /38 :(15%) 48 : 1 ph (1/1) -- tail: 30 ||•| ||||
## 38 /39 :(15.4%) 50 : 1 ph (1/1) -- tail: 30 |||| |••|
## 39 /40 :(15.8%) 51 : 1 ph (1/1) -- tail: 30 •||| •|||
## 40 /41 :(16.2%) 52 : 1 ph (1/1) -- tail: 30 •••| •••|
## 41 /42 :(16.6%) 53 : 1 ph (1/1) -- tail: 30 ||•| ||||
## 42 /43 :(17%) 54 : 1 ph (1/1) -- tail: 30 •••• |••|
## 43 /44 :(17.4%) 55 : 3 ph (1/3) -- tail: 30 ||•• ••|• ... ||•• •|••
## 44 /45 :(17.8%) 56 : 1 ph (1/1) -- tail: 30 ||•| |•||
## 45 /46 :(18.2%) 57 : 1 ph (1/1) -- tail: 30 ••|| •|••
## 46 /47 :(18.6%) 58 : 1 ph (1/1) -- tail: 30 |||• |||•
## 47 /48 :(19%) 59 : 1 ph (1/1) -- tail: 30 |||| |||•
## 48 /49 :(19.4%) 60 : 1 ph (1/1) -- tail: 30 |||| |||•
## 49 /50 :(19.8%) 62 : 1 ph (1/1) -- tail: 30 ••|| •••|
## 50 /51 :(20.2%) 63 : 1 ph (1/1) -- tail: 30 |||• ||||
## 51 /52 :(20.6%) 64 : 1 ph (1/1) -- tail: 30 ||•| ||||
## 52 /53 :(20.9%) 65 : 1 ph (1/1) -- tail: 30 |||| ||•|
## 53 /54 :(21.3%) 66 : 1 ph (1/1) -- tail: 30 •••• ••|•
## 54 /55 :(21.7%) 67 : 9 ph (1/9) -- tail: 30 •||• •||• ... •||• |•|•
## 55 /56 :(22.1%) 68 : 1 ph (1/1) -- tail: 30 •••• •••|
## 56 /57 :(22.5%) 69 : 2 ph (1/2) -- tail: 30 •|•| |||| |••| ||||
## 57 /58 :(22.9%) 70 : 2 ph (1/2) -- tail: 30 |•|| ••|• ||•| ••|•
## 58 /59 :(23.3%) 71 : 1 ph (1/1) -- tail: 30 •••| ••||
## 59 /60 :(23.7%) 72 : 1 ph (1/1) -- tail: 30 •••| ••|•
## 60 /61 :(24.1%) 73 : 1 ph (1/1) -- tail: 30 •••| ••|•
## 61 /62 :(24.5%) 74 : 1 ph (1/1) -- tail: 30 •••• •••|
## 62 /63 :(24.9%) 75 : 1 ph (1/1) -- tail: 30 |||| |||•
## 63 /64 :(25.3%) 76 : 1 ph (1/1) -- tail: 30 •••• ••|•
## 64 /65 :(25.7%) 77 : 1 ph (1/1) -- tail: 30 •••| ••||
## 65 /66 :(26.1%) 78 : 1 ph (1/1) -- tail: 30 •••| ••||
## 66 /67 :(26.5%) 79 : 1 ph (1/1) -- tail: 30 |•|• ||••
## 67 /68 :(26.9%) 80 : 1 ph (1/1) -- tail: 30 •|•| ||||
## 68 /69 :(27.3%) 81 : 1 ph (1/1) -- tail: 30 •||• ••••
## 69 /70 :(27.7%) 84 : 1 ph (1/1) -- tail: 30 ••|• ••••
## 70 /71 :(28.1%) 85 : 1 ph (1/1) -- tail: 30 •••| ••••
## 71 /72 :(28.5%) 86 : 2 ph (1/2) -- tail: 30 ••|• ||•• |••• ||••
## 72 /73 :(28.9%) 87 : 1 ph (1/1) -- tail: 30 •••| ••||
## 73 /74 :(29.2%) 88 : 1 ph (1/1) -- tail: 30 •||| ••||
## 74 /75 :(29.6%) 89 : 1 ph (1/1) -- tail: 30 ••|• ••••
## 75 /76 :(30%) 90 : 1 ph (1/1) -- tail: 30 |||| ||•|
## 76 /77 :(30.4%) 92 : 1 ph (1/1) -- tail: 31 •|•| ••••
## 77 /78 :(30.8%) 93 : 1 ph (1/1) -- tail: 32 ••|• ••••
## 78 /79 :(31.2%) 94 : 1 ph (1/1) -- tail: 33 ||•| ••||
## 79 /80 :(31.6%) 95 : 3 ph (1/3) -- tail: 34 •••| ••|• ... ••|• ••|•
## 80 /81 :(32%) 96 : 4 ph (1/4) -- tail: 35 •||| |||• ... |•|| |||•
## 81 /82 :(32.4%) 97 : 1 ph (1/1) -- tail: 36 |••• •••|
## 82 /83 :(32.8%) 98 : 1 ph (1/1) -- tail: 37 •••• ••|•
## 83 /84 :(33.2%) 99 : 1 ph (1/1) -- tail: 38 |||| ||•|
## 84 /85 :(33.6%) 100: 1 ph (1/1) -- tail: 39 |••• •••|
## 85 /86 :(34%) 101: 1 ph (1/1) -- tail: 40 ••|• ••••
## 86 /87 :(34.4%) 102: 1 ph (1/1) -- tail: 41 |••• •••|
## 87 /88 :(34.8%) 103: 1 ph (1/1) -- tail: 42 ••|• ••••
## 88 /89 :(35.2%) 104: 1 ph (1/1) -- tail: 43 •||| |||•
## 89 /90 :(35.6%) 105: 1 ph (1/1) -- tail: 44 •••• ••|•
## 90 /91 :(36%) 106: 1 ph (1/1) -- tail: 45 ||•| ||||
## 91 /92 :(36.4%) 107: 1 ph (1/1) -- tail: 46 |•|• ||||
## 92 /93 :(36.8%) 108: 1 ph (1/1) -- tail: 47 •••• •••|
## 93 /94 :(37.2%) 109: 1 ph (1/1) -- tail: 48 •••• •••|
## 94 /95 :(37.5%) 110: 1 ph (1/1) -- tail: 49 ||•| ||||
## 95 /96 :(37.9%) 111: 1 ph (1/1) -- tail: 50 |•|• |||•
## 96 /97 :(38.3%) 112: 1 ph (1/1) -- tail: 51 •||| ••||
## 97 /98 :(38.7%) 113: 1 ph (1/1) -- tail: 52 •||| ••||
## 98 /99 :(39.1%) 114: 1 ph (1/1) -- tail: 53 •|•| ••••
## 99 /100:(39.5%) 115: 1 ph (1/1) -- tail: 54 •||| ••||
## 100/101:(39.9%) 116: 1 ph (1/1) -- tail: 55 •||| ••||
## 101/102:(40.3%) 117: 1 ph (1/1) -- tail: 56 |••• ||••
## 102/103:(40.7%) 118: 1 ph (1/1) -- tail: 57 ||•| ||||
## 103/104:(41.1%) 119: 1 ph (1/1) -- tail: 58 |••• ||••
## 104/105:(41.5%) 120: 4 ph (1/4) -- tail: 59 |||| •||| ... |||| |•||
## 105/106:(41.9%) 121: 1 ph (1/1) -- tail: 33 ||•| ||||
## 106/107:(42.3%) 122: 2 ph (1/2) -- tail: 34 |••• •|•• |••• |•••
## 107/108:(42.7%) 123: 1 ph (1/1) -- tail: 35 •||| •|||
## 108/109:(43.1%) 124: 4 ph (1/4) -- tail: 36 •••| •|•• ... ••|• •|••
## 109/110:(43.5%) 125: 1 ph (1/1) -- tail: 30 •||| •|||
## 110/111:(43.9%) 126: 3 ph (1/3) -- tail: 30 |•|| ||•• ... ||•| ||••
## 111/112:(44.3%) 127: 1 ph (1/1) -- tail: 30 •••| •|••
## 112/113:(44.7%) 128: 1 ph (1/1) -- tail: 30 |•|• |•••
## 113/114:(45.1%) 129: 1 ph (1/1) -- tail: 30 •||| •|||
## 114/115:(45.5%) 130: 1 ph (1/1) -- tail: 30 •••| •|••
## 115/116:(45.8%) 131: 1 ph (1/1) -- tail: 30 |••• |•••
## 116/117:(46.2%) 132: 1 ph (1/1) -- tail: 30 |••| |•••
## 117/118:(46.6%) 133: 1 ph (1/1) -- tail: 30 |•|| ||••
## 118/119:(47%) 134: 1 ph (1/1) -- tail: 30 •||• •|||
## 119/120:(47.4%) 135: 1 ph (1/1) -- tail: 30 •||• •|||
## 120/121:(47.8%) 136: 1 ph (1/1) -- tail: 30 •||• •|||
## 121/122:(48.2%) 137: 1 ph (1/1) -- tail: 30 •••• •|••
## 122/123:(48.6%) 138: 1 ph (1/1) -- tail: 30 |••| ||••
## 123/124:(49%) 139: 1 ph (1/1) -- tail: 30 |••| |•••
## 124/125:(49.4%) 140: 1 ph (1/1) -- tail: 30 •||• ||||
## 125/126:(49.8%) 141: 1 ph (1/1) -- tail: 30 •||• ||||
## 126/127:(50.2%) 142: 1 ph (1/1) -- tail: 31 |••| •|••
## 127/128:(50.6%) 143: 1 ph (1/1) -- tail: 32 ||•| ||||
## 128/129:(51%) 144: 1 ph (1/1) -- tail: 33 |••| ||••
## 129/130:(51.4%) 145: 1 ph (1/1) -- tail: 34 •|•• •|||
## 130/131:(51.8%) 146: 1 ph (1/1) -- tail: 35 •|•• ••||
## 131/132:(52.2%) 147: 1 ph (1/1) -- tail: 36 •||• |•||
## 132/133:(52.6%) 148: 1 ph (1/1) -- tail: 37 •||• ••||
## 133/134:(53%) 149: 1 ph (1/1) -- tail: 38 •••• •|••
## 134/135:(53.4%) 150: 6 ph (1/6) -- tail: 39 •||| ••|| ... •||| |••|
## 135/136:(53.8%) 151: 1 ph (1/1) -- tail: 30 •||| |•|•
## 136/137:(54.2%) 152: 1 ph (1/1) -- tail: 30 |••• •|•|
## 137/138:(54.5%) 153: 1 ph (1/1) -- tail: 30 |•|• ||•|
## 138/139:(54.9%) 154: 1 ph (1/1) -- tail: 30 ••|• ||•|
## 139/140:(55.3%) 155: 1 ph (1/1) -- tail: 30 ••|• ||•|
## 140/141:(55.7%) 156: 1 ph (1/1) -- tail: 30 |••• •|•|
## 141/142:(56.1%) 157: 1 ph (1/1) -- tail: 30 •||| |•|•
## 142/143:(56.5%) 158: 1 ph (1/1) -- tail: 30 |••• •|•|
## 143/144:(56.9%) 159: 1 ph (1/1) -- tail: 30 |•|• ||•|
## 144/145:(57.3%) 160: 1 ph (1/1) -- tail: 30 |••• •|•|
## 145/146:(57.7%) 161: 1 ph (1/1) -- tail: 30 •|•• ••••
## 146/147:(58.1%) 162: 2 ph (1/2) -- tail: 30 •••| ••|• •|•• ••|•
## 147/148:(58.5%) 163: 1 ph (1/1) -- tail: 30 |||• ||•|
## 148/149:(58.9%) 164: 1 ph (1/1) -- tail: 30 •||| |•|•
## 149/150:(59.3%) 165: 1 ph (1/1) -- tail: 30 •||| |•|•
## 150/151:(59.7%) 166: 1 ph (1/1) -- tail: 30 |••• •|•|
## 151/152:(60.1%) 167: 1 ph (1/1) -- tail: 30 |•|• •|•|
## 152/153:(60.5%) 168: 1 ph (1/1) -- tail: 30 •|•• ••••
## 153/154:(60.9%) 169: 1 ph (1/1) -- tail: 30 •||| |•|•
## 154/155:(61.3%) 170: 1 ph (1/1) -- tail: 30 |•|• ||•|
## 155/156:(61.7%) 171: 1 ph (1/1) -- tail: 30 |••• •|•|
## 156/157:(62.1%) 172: 1 ph (1/1) -- tail: 30 ••|• ••••
## 157/158:(62.5%) 173: 1 ph (1/1) -- tail: 30 |•|• ||•|
## 158/159:(62.8%) 174: 1 ph (1/1) -- tail: 30 •||| |•|•
## 159/160:(63.2%) 175: 1 ph (1/1) -- tail: 30 •|•| ••|•
## 160/161:(63.6%) 176: 1 ph (1/1) -- tail: 30 •|•| ••|•
## 161/162:(64%) 177: 1 ph (1/1) -- tail: 30 •||| |•|•
## 162/163:(64.4%) 178: 1 ph (1/1) -- tail: 30 •||| |•|•
## 163/164:(64.8%) 179: 1 ph (1/1) -- tail: 30 •••• •|•|
## 164/165:(65.2%) 180: 1 ph (1/1) -- tail: 31 •••• •|•|
## 165/166:(65.6%) 181: 1 ph (1/1) -- tail: 32 |••• •|•|
## 166/167:(66%) 182: 1 ph (1/1) -- tail: 33 |•|• ||•|
## 167/168:(66.4%) 183: 1 ph (1/1) -- tail: 34 |••• •|•|
## 168/169:(66.8%) 184: 1 ph (1/1) -- tail: 35 •||| |•|•
## 169/170:(67.2%) 185: 1 ph (1/1) -- tail: 36 •|•| ••|•
## 170/171:(67.6%) 186: 1 ph (1/1) -- tail: 37 •|•| ••|•
## 171/172:(68%) 187: 1 ph (1/1) -- tail: 38 •||| |•|•
## 172/173:(68.4%) 188: 1 ph (1/1) -- tail: 39 •||| ||||
## 173/174:(68.8%) 189: 1 ph (1/1) -- tail: 40 •••| ••|•
## 174/175:(69.2%) 190: 1 ph (1/1) -- tail: 41 •••| ••|•
## 175/176:(69.6%) 191: 1 ph (1/1) -- tail: 42 |||• ||•|
## 176/177:(70%) 192: 2 ph (1/2) -- tail: 43 •••• ••|• •••• |•••
## 177/178:(70.4%) 193: 1 ph (1/1) -- tail: 44 |||• •|•|
## 178/179:(70.8%) 194: 1 ph (1/1) -- tail: 45 •••• |•••
## 179/180:(71.1%) 195: 2 ph (1/2) -- tail: 46 ••|| |•|• •|•| |•|•
## 180/181:(71.5%) 196: 1 ph (1/1) -- tail: 47 ••|| |•|•
## 181/182:(71.9%) 197: 1 ph (1/1) -- tail: 48 |||| ||•|
## 182/183:(72.3%) 198: 1 ph (1/1) -- tail: 49 |||| ||•|
## 183/184:(72.7%) 199: 1 ph (1/1) -- tail: 50 ••|• |•|•
## 184/185:(73.1%) 200: 1 ph (1/1) -- tail: 51 •••• ••|•
## 185/186:(73.5%) 201: 1 ph (1/1) -- tail: 52 ••|| |•|•
## 186/187:(73.9%) 202: 1 ph (1/1) -- tail: 53 |||| ||•|
## 187/188:(74.3%) 203: 1 ph (1/1) -- tail: 54 •••• ••|•
## 188/189:(74.7%) 204: 1 ph (1/1) -- tail: 55 •••• |•|•
## 189/190:(75.1%) 205: 1 ph (1/1) -- tail: 56 |||• •|•|
## 190/191:(75.5%) 206: 1 ph (1/1) -- tail: 57 ••|• |•|•
## 191/192:(75.9%) 207: 1 ph (1/1) -- tail: 58 ••|| |•|•
## 192/193:(76.3%) 208: 1 ph (1/1) -- tail: 59 •••• ••|•
## 193/194:(76.7%) 209: 1 ph (1/1) -- tail: 60 ••|• ••••
## 194/195:(77.1%) 210: 1 ph (1/1) -- tail: 61 ||•| •|•|
## 195/196:(77.5%) 211: 1 ph (1/1) -- tail: 62 |||• ||||
## 196/197:(77.9%) 212: 1 ph (1/1) -- tail: 63 •••| |•|•
## 197/198:(78.3%) 213: 1 ph (1/1) -- tail: 64 |||• •|||
## 198/199:(78.7%) 214: 1 ph (1/1) -- tail: 65 ••|• |•|•
## 199/200:(79.1%) 215: 1 ph (1/1) -- tail: 66 |||• ||||
## 200/201:(79.4%) 216: 1 ph (1/1) -- tail: 67 |||| ||•|
## 201/202:(79.8%) 217: 1 ph (1/1) -- tail: 68 |||• ||||
## 202/203:(80.2%) 218: 1 ph (1/1) -- tail: 69 ||•• •|•|
## 203/204:(80.6%) 219: 1 ph (1/1) -- tail: 70 ••|• ••••
## 204/205:(81%) 220: 1 ph (1/1) -- tail: 71 ||•• •|•|
## 205/206:(81.4%) 221: 1 ph (1/1) -- tail: 72 ••|| |•|•
## 206/207:(81.8%) 222: 1 ph (1/1) -- tail: 73 ••|| |•|•
## 207/208:(82.2%) 223: 1 ph (1/1) -- tail: 74 •••| ••|•
## 208/209:(82.6%) 224: 1 ph (1/1) -- tail: 75 •••| ••|•
## 209/210:(83%) 225: 1 ph (1/1) -- tail: 76 ••|| |•|•
## 210/211:(83.4%) 226: 1 ph (1/1) -- tail: 77 ••|• |•••
## 211/212:(83.8%) 227: 1 ph (1/1) -- tail: 78 ••|• |•|•
## 212/213:(84.2%) 228: 1 ph (1/1) -- tail: 79 ||•| •|•|
## 213/214:(84.6%) 229: 1 ph (1/1) -- tail: 80 ••|• |•••
## 214/215:(85%) 230: 1 ph (1/1) -- tail: 81 ••|• |•|•
## 215/216:(85.4%) 231: 1 ph (1/1) -- tail: 82 ••|• |•|•
## 216/217:(85.8%) 232: 1 ph (1/1) -- tail: 83 ||•| •|||
## 217/218:(86.2%) 233: 1 ph (1/1) -- tail: 84 ||•| •|•|
## 218/219:(86.6%) 234: 1 ph (1/1) -- tail: 85 ••|• ••••
## 219/220:(87%) 235: 1 ph (1/1) -- tail: 86 |||| •|||
## 220/221:(87.4%) 236: 1 ph (1/1) -- tail: 87 ||•| ||||
## 221/222:(87.7%) 237: 1 ph (1/1) -- tail: 88 |||| ||•|
## 222/223:(88.1%) 238: 1 ph (1/1) -- tail: 89 •••• ••|•
## 223/224:(88.5%) 239: 1 ph (1/1) -- tail: 90 |||| •|||
## 224/225:(88.9%) 240: 1 ph (1/1) -- tail: 91 ||•| ||||
## 225/226:(89.3%) 241: 1 ph (1/1) -- tail: 92 |||| ||•|
## 226/227:(89.7%) 242: 2 ph (1/2) -- tail: 93 ||•| •••| ||•| •|••
## 227/228:(90.1%) 243: 1 ph (1/1) -- tail: 55 |||| ||•|
## 228/229:(90.5%) 244: 1 ph (1/1) -- tail: 56 ||•• •••|
## 229/230:(90.9%) 245: 1 ph (1/1) -- tail: 57 ••|• ••••
## 230/231:(91.3%) 246: 1 ph (1/1) -- tail: 58 ||•• •••|
## 231/232:(91.7%) 248: 2 ph (1/2) -- tail: 59 |||| |•|| |||| |||•
## 232/233:(92.1%) 249: 1 ph (1/1) -- tail: 60 •||| |||•
## 233/234:(92.5%) 250: 1 ph (1/1) -- tail: 30 ••|• ••••
## 234/235:(92.9%) 251: 1 ph (1/1) -- tail: 30 |||| |||•
## 235/236:(93.3%) 252: 3 ph (1/3) -- tail: 30 ••|| ••|• ... •|•| ••|•
## 236/237:(93.7%) 253: 1 ph (1/1) -- tail: 30 ••|• ••••
## 237/238:(94.1%) 254: 1 ph (1/1) -- tail: 30 ••|| |||•
## 238/239:(94.5%) 255: 1 ph (1/1) -- tail: 30 ||•• •••|
## 239/240:(94.9%) 256: 1 ph (1/1) -- tail: 30 •••• |•••
## 240/241:(95.3%) 257: 1 ph (1/1) -- tail: 30 ||•| ||||
## 241/242:(95.7%) 258: 2 ph (1/2) -- tail: 30 |•|| ||•| |||• ||•|
## 242/243:(96%) 259: 2 ph (1/2) -- tail: 30 •||| |||• |•|| |||•
## 243/244:(96.4%) 260: 1 ph (1/1) -- tail: 30 ••|| |||•
## 244/245:(96.8%) 261: 1 ph (1/1) -- tail: 30 |||| |•||
## 245/246:(97.2%) 262: 1 ph (1/1) -- tail: 30 ||•• •••|
## 246/247:(97.6%) 263: 1 ph (1/1) -- tail: 30 ||•• •••|
## 247/248:(98%) 264: 1 ph (1/1) -- tail: 30 |•|| |||•
## 248/249:(98.4%) 265: 2 ph (1/2) -- tail: 30 ||•• •|•| ||•• |••|
## 249/250:(98.8%) 266: 1 ph (1/1) -- tail: 30 |||| |•||
## 250/251:(99.2%) 267: 1 ph (1/1) -- tail: 30 •••| |||•
## 251/252:(99.6%) 268: 1 ph (1/1) -- tail: 30 ||•• •••|
## 252/253:(100%) 269: 1 ph (1/1) -- tail: 30 ||•• •••|## ══════════════════════════════════════════════════════ Reestimating final recombination fractions ══## Markers in the initial sequence: 253
## Mapped markers : 252 (99.6%)## ════════════════════════════════════════════════════════════════════════════════════════════════════plot(single_chrom)
Parallel construction of all 7 LGs
## Performing parallel computation of hidden markov
phasing_and_hmm_rf <- function(X){
fl <- paste0("output_map_ch_", unique(X$seq$chrom), ".txt")
sink(fl)
map <- est_rf_hmm_sequential(input.seq = X$seq,
start.set = 3,
thres.twopt = 10,
thres.hmm = 50,
extend.tail = 30,
twopt = X$tpt,
verbose = TRUE,
phase.number.limit = 20,
sub.map.size.diff.limit = 5)
sink()
return(map)
}
# 5% genotyping error
my.error.func<-function(X){
x<-est_full_hmm_with_global_error(input.map = X,
error = 0.05,
tol = 10e-4,
verbose = FALSE)
return(x)
}
#upgma
#hidden markov based off physical location
ptm <- proc.time()
cl <- parallel::makeCluster(7)
parallel::clusterEvalQ(cl, require(mappoly))## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUEparallel::clusterExport(cl, "dat")
MAPs <- parallel::parLapply(cl,LGS, phasing_and_hmm_rf)
parallel::stopCluster(cl)
# account for error
cl <- parallel::makeCluster(7)
parallel::clusterEvalQ(cl, require(mappoly))## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUEparallel::clusterExport(cl, "dat")
MAPs.err <- parallel::parLapply(cl,MAPs,my.error.func)
parallel::stopCluster(cl)
par(mfrow=c(2,1))
plot_map_list(MAPs, col = "ggstyle")
plot_map_list(MAPs.err, col = "ggstyle")
proc.time() - ptm## user system elapsed
## 1.58 0.97 1404.33After constructing the maps a 5% global error rate is applied to the chain to allow for the chain to make small changes that correct for possible genotyping error. The difference is seen the the picture above where the second map is shorter indicating a better map.
plot_genome_vs_map(MAPs.err, same.ch.lg = T)
summary_maps(MAPs.err)##
## Your dataset contains removed (redundant) markers. Once finished the maps, remember to add them back with the function 'update_map'.## LG Genomic sequence Map length (cM) Markers/cM Simplex Double-simplex Multiplex Total Max gap
## 1 1 1 70.52 3.57 80 44 128 252 2.75
## 2 2 2 111.38 2.83 117 44 154 315 6.92
## 3 3 3 66.91 3.18 75 49 89 213 4.5
## 4 4 4 66.34 3.77 84 40 126 250 3.53
## 5 5 5 72.84 3.06 87 47 89 223 3.12
## 6 6 6 80.2 3.69 131 47 118 296 3.96
## 7 7 7 81.92 3.41 88 104 87 279 2.95
## 8 Total <NA> 550.11 3.36 662 375 791 1828 3.96Genotypic probabilites for use in QTL mapping
Genotypic probabilites will be used in QTL mapping thus we need to save this for use in next session
genoprob.err <- vector("list", 7)
for(i in 1:7)
genoprob.err[[i]] <- calc_genoprob_error(input.map = MAPs.err[[i]], error = 0.05)## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......
## Ploidy level:4
## Number of individuals:157
## ..................................................
## ..................................................
## ..................................................
## .......saveRDS(genoprob.err, file="genoprob.err.RDS")homolog prob
homolog_probs = calc_homologprob(genoprob.err)##
## Linkage group 1 ...
## Linkage group 2 ...
## Linkage group 3 ...
## Linkage group 4 ...
## Linkage group 5 ...
## Linkage group 6 ...
## Linkage group 7 ...plot(homolog_probs, lg = 1, ind = 1)HMM chain changing and imputing genotypes
The HMM chain will make small changes to the progeny genotype calls
(dosage) so it fits the better in the map. As you increase the global
error rate, you allow the chain more freedom to make more changes. This
is represented in the figure below. 
QTLpoly for QTL analysis
MAPpoly’s sister software QTLpoly was made to take the outputs of MAPpoly to do QTL mapping. There are two main ways in QTLpoly REMIM random-effect multiple interval mapping and FEIM fixed effect interval mapping.
Import data
First we need the to read in the genotypic probabilites file
genoprob.err object that was generated in the previous
linkage mapping session.
library(qtlpoly)
genoprob.err = readRDS(file="genoprob.err.RDS")
pheno = read.csv("../data/blackspot_cercospora_defoliation_RRD.csv", row.names = 1)
data = read_data(ploidy = 4, geno.prob = genoprob.err, pheno = pheno, step = 1)## Reading the following data:
## Ploidy level: 4
## No. individuals: 157
## No. linkage groups: 7
## Step size: 1 cM
## Map size: 550.11 cM (448 positions)
## No. phenotypes: 5# this step takes a long time so do not run import it
#data.sim = simulate_qtl(data = data, mu = 0, h2.qtl = NULL, var.error = 1, n.sim = 1000, missing = TRUE, seed = 123)
#score.null = null_model(data = data.sim$results, n.clusters = 7, plot = NULL)Didactic example of one trait for teaching REMIM
The REMIM method consists of a few steps: this is a simplification for more math terms go here: https://guilherme-pereira.github.io/QTLpoly/1-tutorial#Perform_QTL_detection 1. Null model: model starts with no QTL 2. Forward search: QTL are added one at a time conditional to any previous if any QTL already in the model 3. Model optimization: each QTL is then tested again conditional to all other QTL under a more stringent threshold - steps 2 and 3 are repeated until no more QTL are added and dropped from model. 4. QTL profiling: score statistics are updated conditional to QTL that are retained
1. Null model
null.mod <- null_model(data = data, pheno.col = 1 ,n.clusters = 4, plot = "null")## INFO: Using 4 CPUs for calculation
##
## Null model for trait 1 'blackspot'## Calculation took 14.44 secondsprint(null.mod)## This is an object of class 'qtlpoly.null'
##
## * Trait 1 'blackspot'
## There are no QTL in the model2. Forward search
search.mod <- search_qtl(data = data, model = null.mod, w.size = 15, sig.fwd = 0.01, n.clusters = 4, plot = "search")## INFO: Using 4 CPUs for calculation
##
## Forward search for trait 1 'blackspot'; there are no QTL in the model
## QTL was found on LG 3 at 49.13 cM (position number 178)## QTL was found on LG 5 at 36.71 cM (position number 282)## QTL was found on LG 6 at 16.08 cM (position number 321)## QTL was found on LG 1 at 60.06 cM (position number 53)## No more QTL were found. A putative QTL on LG 4 at 55.13 cM (position number 239) did not reach the threshold; its p-value was 0.01944
## Calculation took 249.38 secondsprint(search.mod)## This is an object of class 'qtlpoly.search'
##
## * Trait 1 'blackspot'
## LG Pos Nmrk Mrk Score Pval
## 1 3 49.13 178 Affx-86777792 271.66 5.10e-09
## 2 5 36.71 282 Affx-86838636 207.94 5.58e-07
## 3 6 16.08 321 Affx-86817386 60.56 1.50e-02
## 4 1 60.06 53 Affx-86815556 53.16 2.24e-023. Optimize QTL
optimize.mod <- optimize_qtl(data = data, model = search.mod, sig.bwd = 1e-04, n.clusters = 4, plot = "optimize")## INFO: Using 4 CPUs for calculation
##
## Model optimization for trait 1 'blackspot'; there are 4 QTL in the model already
## Refining QTL positions## ... 178## ... 282## ... 321## ... 53
## Excluding non-significant QTL ... 321 ... 53
## Refining QTL positions## ... 178## ... 282
## Calculation took 86.19 secondsprint(optimize.mod)## This is an object of class 'qtlpoly.optimize'
##
## * Trait 1 'blackspot'
## LG Pos Nmrk Mrk Score Pval
## 1 3 49.13 178 Affx-86777792 286.81 1.02e-09
## 2 5 36.71 282 Affx-86838636 241.71 2.42e-083. Profile QTL
profile.mod <- profile_qtl(data = data, model = optimize.mod, d.sint = 1.5, polygenes = FALSE, n.clusters = 4, plot = "profile")## INFO: Using 4 CPUs for calculation
##
## QTL profile for trait 1 'blackspot'; there are 2 QTL in the model
## Profiling QTL ... 178## ... 282## Calculation took 68.17 secondsprint(profile.mod)## This is an object of class 'qtlpoly.profile'
##
## * Trait 1 'blackspot'
## LG Pos Nmrk Mrk Score Pval
## 1 3 49.13 178 Affx-86777792 286.81 1.02e-09
## 2 5 36.71 282 Affx-86838636 241.71 2.42e-08print(profile.mod, sint="lower")## This is an object of class 'qtlpoly.profile'
##
## * Trait 1 'blackspot'
## LG Pos_lower Nmrk_lower Mrk_lower Score_lower Pval_lower
## 1 3 47.29 176 Affx-86800060 261.97 7.08e-09
## 2 5 34.07 281 Affx-86782784 200.45 4.60e-07print(profile.mod, sint="upper")## This is an object of class 'qtlpoly.profile'
##
## * Trait 1 'blackspot'
## LG Pos_upper Nmrk_upper Mrk_upper Score_upper Pval_upper
## 1 3 53.45 181 Affx-86786863 250.81 1.25e-08
## 2 5 42.41 286 Affx-86839796 218.69 2.07e-07Plotting the QTL profiles
plot_profile(data=data,model=profile.mod, sup.int = T) 
Run all of the traits
Because the score.null object (genome wide significance levels) take so long to calculate, I have provided done it in advance. If you have a new map, or if you have a new experiment, you will need to recalculate.
score.null = readRDS("../data/score.null.RDS")
remim.mod = remim(data, score.null = score.null, w.size = 15, sig.fwd = 0.2, sig.bwd = 0.05, n.clusters = 7)## INFO: Using 7 CPUs for calculation
##
## REMIM for trait 1 'blackspot'
## QTL was found on LG 3 at 49.13 cM (position number 178)
## QTL was found on LG 5 at 36.71 cM (position number 282)
## No more QTL were found. A putative QTL on LG 6 at 16.08 cM (position number 321) did not reach the threshold; its p-value was 0.00385
## Refining QTL positions ... 178 ... 282
## Profiling QTL ... 178 ... 282
## Calculation took 157.11 seconds
##
## REMIM for trait 2 'cercospora'
## QTL was found on LG 1 at 11.57 cM (position number 12)
## QTL was found on LG 5 at 38.12 cM (position number 283)
## QTL was found on LG 6 at 8.34 cM (position number 314)
## No more QTL were found. A putative QTL on LG 2 at 88.17 cM (position number 121) did not reach the threshold; its p-value was 0.02052
## Refining QTL positions ... 12 ... 283 ... 314
## Excluding non-significant QTL ... 314
## Refining QTL positions ... 12 ... 283
## Profiling QTL ... 12 ... 283
## Calculation took 237.94 seconds
##
## REMIM for trait 3 'defoliation'
## QTL was found on LG 3 at 49.13 cM (position number 178)
## QTL was found on LG 5 at 36.71 cM (position number 282)
## QTL was found on LG 1 at 26.3 cM (position number 22)
## No more QTL were found. A putative QTL on LG 6 at 80.19 cM (position number 376) did not reach the threshold; its p-value was 0.00394
## Refining QTL positions ... 178 ... 282 ... 22
## Excluding non-significant QTL ... 178
## Refining QTL positions ... 282 ... 22
## Profiling QTL ... 22 ... 282
## Calculation took 236.12 seconds
##
## REMIM for trait 4 'RRD_severity'
## QTL was found on LG 5 at 4.99 cM (position number 254)
## QTL was found on LG 7 at 25.61 cM (position number 400)
## No more QTL were found. A putative QTL on LG 2 at 83.16 cM (position number 116) did not reach the threshold; its p-value was 0.00664
## Refining QTL positions ... 258 ... 382
## Refining QTL positions ... 258 ... 382
## Profiling QTL ... 258 ... 382
## Calculation took 90 seconds
##
## REMIM for trait 5 'RRV_Ct'
## QTL was found on LG 5 at 17.18 cM (position number 265)
## No more QTL were found. A putative QTL on LG 3 at 2.85 cM (position number 142) did not reach the threshold; its p-value was 0.04121
## Refining QTL positions ... 265
## Profiling QTL ... 265
## Calculation took 39.8 secondsplot_profile(data, remim.mod, sup.int = T)
plot_profile(data, remim.mod, sup.int = T, grid=T)
plot_sint(data,remim.mod)
fit_model() calculates all the QTL’s variance components
and estimates heritability of the QTL (comparable to the R-squared of
the QTL)
fitted.mod <- fit_model(data = data, model = remim.mod, probs = "joint", polygenes = "none")## There are 2 QTL in the model for trait 1 'blackspot'. Fitting model... Done!
##
## There are 2 QTL in the model for trait 2 'cercospora'. Fitting model... Done!
##
## There are 2 QTL in the model for trait 3 'defoliation'. Fitting model... Done!
##
## There are 2 QTL in the model for trait 4 'RRD_severity'. Fitting model... Done!
##
## There is 1 QTL in the model for trait 5 'RRV_Ct'. Fitting model... Done!summary(fitted.mod) #some reason this does not display right please run in the console to see results## This is an object of class 'qtlpoly.fitted'
##
## * Trait 1 'blackspot'
## LG Pos Nmrk Intercept Var(g) Var(e) h2
## 1 3 49.13 178 NA 0.1628624 NA 0.3798968
## 2 5 36.71 282 NA 0.1132065 NA 0.2640683
## 3 NA NA NA 3.514389 NA 0.1526328 0.643965
##
## * Trait 2 'cercospora'
## LG Pos Nmrk Intercept Var(g) Var(e) h2
## 1 1 11.57 12 NA 2.120141 NA 0.4179725
## 2 5 38.12 283 NA 1.364544 NA 0.2690113
## 3 NA NA NA 2.050893 NA 1.587756 0.6869838
##
## * Trait 3 'defoliation'
## LG Pos Nmrk Intercept Var(g) Var(e) h2
## 1 1 26.3 22 NA 0.1589865 NA 0.1873592
## 2 5 36.71 282 NA 0.1756531 NA 0.2070002
## 3 NA NA NA 4.076364 NA 0.5139255 0.3943594
##
## * Trait 4 'RRD_severity'
## LG Pos Nmrk Intercept Var(g) Var(e) h2
## 1 5 8.6 258 NA 0.3318958 NA 0.3478961
## 2 7 5.13 382 NA 0.1464115 NA 0.1534698
## 3 NA NA NA 0.8961831 NA 0.4757011 0.5013659
##
## * Trait 5 'RRV_Ct'
## LG Pos Nmrk Intercept Var(g) Var(e) h2
## 1 5 17.18 265 30.04077 14.86612 50.55313 0.2272438qtl_effects
qtl_effects() shows the user for each QTL what parental
homolog is responsible for increasing or decreasing the phenotypic
average.
est.effects <- qtl_effects(ploidy = 4, fitted = fitted.mod)## There are 2 QTL in the model for trait 1 'blackspot'. Computing effects for QTL ... 178 ... 282. Done!
##
## There are 2 QTL in the model for trait 2 'cercospora'. Computing effects for QTL ... 12 ... 283. Done!
##
## There are 2 QTL in the model for trait 3 'defoliation'. Computing effects for QTL ... 22 ... 282. Done!
##
## There are 2 QTL in the model for trait 4 'RRD_severity'. Computing effects for QTL ... 258 ... 382. Done!
##
## There is 1 QTL in the model for trait 5 'RRV_Ct'. Computing effects for QTL ... 265. Done!plot(est.effects)
Save results for visualization in Viewpoly
What we need to do is now save all the objects so that we can now look at a more beautiful output for mining your results.
# linkage map file
save(MAPs.err,file = "../viewpoly_files/MAPs.err.RData")
save(data, file = "../viewpoly_files/QTLpoly_data.RData")
save(remim.mod, file = "../viewpoly_files/QTLpoly_remim.mod.RData")
save(est.effects, file = "../viewpoly_files/QTLpoly_est.effects.RData")
save(fitted.mod, file = "../viewpoly_files/QTLpoly_fitted.mod.RData")