fitOneMarker: Function to fit multiple mixture models to signal ratios of a...
In fitPoly: Genotype Calling for Bi-Allelic Marker Assays
Description
Usage
Arguments
Details
Value
Examples
Description
This function takes a data frame with allele signal ratios for
multiple bi-allelic markers and samples, and fits multiple mixture models to
a selected marker. It returns a list, reporting on the performance of these
models, selecting the best one based on the BIC criterion, optionally
plotting results.
Usage
1
2
3
4
5
6
fitOneMarker(ploidy, marker, data, diplo=NULL, select=TRUE,
diploselect=TRUE, pop.parents=NULL, population=NULL, parentalPriors=NULL,
samplePriors=NULL, startmeans=NULL, maxiter=40, maxn.bin=200, nbin=200,
sd.threshold=0.1, p.threshold=0.99, call.threshold=0.6, peak.threshold=0.85,
try.HW=TRUE, dip.filter=1, sd.target=NA,
plot="none", plot.type="png", plot.dir, sMMinfo=NULL)
Arguments
ploidy
The ploidy level, 2 or higher: 2 for diploids, 3 for triploids
etc.
marker
A marker name of number. Used to select the data for one marker,
referring to the MarkerName column of parameter data. If a number, the number
of the marker based on alphabetic order of the MarkerNames in data.
data
A data frame with the polyploid samples, with (at least) columns
MarkerName, SampleName and ratio, where ratio is the Y-allele signal
divided by the sum of the X- and Y-allele signals: ratio == Y/(X+Y)
diplo
NULL or a data frame like data, with the diploid samples and (a subset
of) the same markers as in data. Genotypic scores for diploid samples are
calculated according to the best-fitting model calculated for the polyploid
samples and therefore may range from 0 (nulliplex) to <ploidy>, with the
expected dosages 0 and <ploidy> for the homozygotes and <ploidy/2> for the
heterozygotes.
Note that diplo can also be used for any other samples that need to be
scored, but that should not affect the fitted models.
select
A logical vector, recycled if shorter than nrow(data):
indicates which rows of data are to be used (default TRUE, i.e. keep all rows)
diploselect
A logical vector like select, matching diplo instead of data
pop.parents
NULL or a data.frame specifying the population structure. The
data frame has 3 columns: the first containing population IDs, the 2nd and 3rd
with the population IDs of the parents of these populations (if F1's) or NA
(if not). The poopulation IDs should match those in parameter population. If
pop.parents is NULL all samples are considered to be in one population, and
parameter population should also be NULL (default).
population
NULL or a data.frame specifying to which population each
sample belongs. The data frame has two columns, the first containing
the SampleName (containing all SampleNames occurring in data),
the second column containing population IDs that match pop.parents. In both
columns NA values are not allowed. Parameters pop.parents and population
should both be NULL (default) or both be specified.
parentalPriors
NULL or a data frame specifying the prior dosages for
the parental populations. The data frame has one column MarkerName
followed by one column for each F1 parental population. Column names (except
first) are population IDs matching the parental populations in pop.parents.
In case there is just one F1 population in pop.parents, it is possible to
have two columns for both parental populations instead of one (allowing two
specify two different prior dosages); in that case both columns for each
parent have the same caption. Each row specifies the priors for
one marker. The contents of the data frame are dosages, as integers from 0
to <ploidy>; NA values are allowed.
Note: when reading the data frame with read.table or read.csv, set
check.names=FALSE so column names (population IDs) are not changed.
samplePriors
NULL or a data.frame specifying prior dosages for individual
samples. The first column called MarkerName is followed by one column per
sample; not all samples in data need to have a column here, only
those samples for which prior dosages for one or more markers are available.
Each row specifies the priors for one marker. The contents of the data frame
are dosages, as integers from 0 to <ploidy>; NA values are allowed.
Note: when reading the data frame with read.table or read.csv, set
check.names=FALSE so column names (population IDs) are not changed.
startmeans
NULL or a data.frame specifying the prior means of
the mixture distributions. The data frame has one column MarkerName,
followed by <ploidy+1> columns with the prior means on the original
(untransformed) scale. Each row specifies the
means for one marker in strictly ascending order (all means NA is allowed, but
markers without start means can also be omitted).
maxiter
A single integer, passed to CodomMarker, see there for explanation
maxn.bin
A single integer, passed to CodomMarker, see there for explanation
nbin
A single integer, passed to CodomMarker, see there for explanation
sd.threshold
The maximum value allowed for the (constant) standard
deviation of each peak on the arcsine - square root transformed scale,
default 0.1. If the optimal model has a larger standard deviation the marker
is rejected. Set to a large value (e.g. 1) to disable this filter.
p.threshold
The minimum P-value required to assign a genotype (dosage)
to a sample; default 0.99. If the P-value for all possible genotypes is less
than p.threshold the sample is assigned genotype NA. Set to 1 to disable
this filter.
call.threshold
The minimum fraction of samples to have genotypes
assigned ("called"); default 0.6. If under the optimal model the fraction of
"called" samples is less than call.threshold the marker is rejected. Set to 0
to disable this filter.
peak.threshold
The maximum allowed fraction of the scored samples that
are in one peak; default 0.85. If any of the possible genotypes (peaks in the
ratio histogram) contains more than peak.threshold of the samples the marker
is rejected (because the remaining samples offers too little information for
reliable model fitting).
try.HW
Logical: if TRUE (default), try models with and without a
constraint on the mixing proportions according to Hardy-Weinberg equilibrium
ratios. If FALSE, only try models without this constraint. Even when the HW
assumption is not applicable, setting try.HW to TRUE often still leads to
a better model. For more details on how try.HW is used see the Details
section.
dip.filter
if 1 (default), select best model only from models
that do not have a dip (a lower peak surrounded by higher peaks: these are not
expected under Hardy-Weinberg equilibrium or in cross progenies). If all
fitted models have a dip still the best of these is selected. If 2, similar,
but if all fitted models have a dip the marker is rejected. If 0, select best
model among all fitted models, including those with a dip.
sd.target
If the fitted standard deviation of the peaks on the
transformed scale is larger than sd.target a penalty is given (see Details);
default NA i.e. no penalty is given.
plot
String, "none" (default), "fitted" or "all". If "fitted" a plot
of the best fitting model and the assigned genotypes is saved with filename
<marker number><marker name>.<plot.type>, preceded by "rejected_" if the
marker was rejected. If "all", small plots of all models are saved to files
(8 per file) with filename
<"plots"><marker number><A..F><marker name>.<plot.type> in addition to the
plot of the best fitting model.
plot.type
String, "png" (default), "emf", "svg" or "pdf". Indicates
format for saving the plots.
plot.dir
String, the directory where to save the plot files. Must be
specified if plot is not "none". Set this to "" to save plot files
in the current working directory.
sMMinfo
NULL (default), for internal use only. Prevents unneeded checking
and recalculation of input parameters when called from saveMarkerModels.
Details
fitOneMarker fits a series of mixture models for the given marker by
repeatedly calling CodomMarker and selects the optimal one. The initial
models vary according to the values of try.HW, pop.parents,
parentalPriors, samplePriors and startmeans:
no pop.parents, try.HW FALSE: 4 models with different constraints
on the means (different or equal X and Y background signal, ratio a linear or
quadratic function of dosage), no restrictions on the mixing proportions
(the fractions of samples in each dosage peak)
no pop.parents, try.HW TRUE: The previous 4 models are fitted and
also 4 models with the same restrictions on the means and the mixing
proportions restricted to Hardy-Weinberg ratios (assuming polysomic
inheritance)
pop.parents specified, no parentalPriors / samplePriors / startmeans,
try.HW FALSE: 4 models
are fitted with the same restrictions on the means as above, but with
different restrictions on the mixing proportions for each population:
no restriction on parental populations, none on accession panels, polysomic
F1 segregation ratios on F1 populations. Additionally 4 models are fitted
with all samples considered as one population, with the same 4 models for
the means and no restrictions on mixing proportions.
pop.parents specified, no parentalPriors / samplePriors / startmeans,
try.HW TRUE: 4 models
are fitted with the same restrictions on the means as above, but with
different restrictions on the mixing proportions for each population:
no restriction on parental populations, HW-ratios for accession panels,
polysomic F1 segregation ratios on F1 populations. Additionally 4 models
are fitted with all samples considered as one population, with the same 4
models for the means and mixing proportions according to HW ratios.
pop.parents and parentalPriors specified, try.HW FALSE: 4 models
are fitted with the same restrictions on the means as above, but with
different restrictions on the mixing proportions for each population:
no restriction on parental populations and the accession panels,
polysomic F1 segregation ratios on F1 populations ignoring the parental
priors. Additionally 4 models are fitted with the same restrictions on the
means and mixing proportions of the accession panels, but where the
mixing proportions of the parental populations are set to (almost) 1 for
the prior dosage and (almost) 0 for all other dosages, and those for the
F1 populations to the polysomic segregation ratios expected for the
parental priors.
pop.parents and parentalPriors specified, try.HW TRUE: same as with
try.HW FALSE, except that the mixing proportions of accession panels are now
restricted to HW ratios.
if parentalPriors and/or samplePriors are specified, these and the
signal ratios of the corresponding samples are (also) used to estimate starting
values of the mixture component means in the EM algorithm. Alternatively
startmeans can be specified directly.
Because convergence to the optimal solution often fails, the models are fitted
with several start values for the <ploidy+1> means of the mixture
distributions: (1) based on initial clustering of the ratios, (2) based on
a uniform distribition from 0.02 to pi/2-0.02 on the asin(sqrt(x)) scale,
and (3) if startmeans are specified or can be calculated from samplePriors
and/or parentalPriors these are used for a third set of model fits.
The main difference between parentalPriors and samplePriors is that
parentalPriors are treated as fixed (and if both parents of an F1 population
have priors, the F1 segregation is also fixed) while samplePriors are only
used to calculate starting ratio means for each dosage. Depending on the
confidence the user has in the prior dosages of the parents they can be
supplied as parentalPriors or samplePriors.
In some cases an additional fit is performed with a modified set of initial
means.
An optimal model is selected based on the Bayesian Information Criterium
(BIC), which takes into account the Log-Likelihood and the number of fitted
parameters of the models. If sd.target is specified and the standard
deviation of the mixture model components is larger than this target a
penalty is applied, making is less likely that that model is selected.
The plots consist of one histogram per (non-parent) population showing the
frequency distribution of the signal ratios of the samples in that population.
The fitted model is shown in green (density and means), and for F1 populations
the samples of parent 1 and 2 are shown as red and blue triangles.
If diplodata are present, a histogram for the diploid samples is plotted
in the top histogram (diploid bars are narrower and gray). The diploid bars
are scaled so the maximum bar is half the maximum polyploid bar. At the
bottom of the plot for the fitted model a rug plot shows the scores of each
sample, while the bottom (red) samples are unscored.
Value
A list with components:
- log
A character vector with the lines of the log text.
- modeldata
A data frame as allmodeldata (see below) with only the
one row with data on the selected model.
- allmodeldata
A data frame with for each tried model one row with
the marker number, marker name, number of samples and (if the marker is
not rejected) data of the fitted model (see below).
- scores
A data frame with the name and data for all samples
(including NA's for the samples that were not selected, see parameter
select), with columns:
marker (the sequential number of the marker (based on alphabetic
order of the marker names in data)
MarkerName
SampleName
ratio (the given ratio from parameter data)
P0 .. P<ploidy> (the probabilities that this sample belongs to each of the
<ploidy+1> mixture components)
maxgeno (0..ploidy, the genotype = mixture component with the highest P
value)
maxP (the P value for this genotype)
geno (the assigned genotype number: same as maxgeno, or NA if
maxP < p.threshold).
- diploscores
A data frame like scores for the samples in the data
frame supplied with argument diplo. If diplo is NA also diploscores will be
NA.
The modeldata and allmodeldata data frames present data on a fitted model.
modeldata presents data on the selected model; allmodeldata lists all
attempted models. Both data frames contain the following columns:
- marker
the sequential number of the marker (based on alphabetic
order of the marker names in data)
- MarkerName
the name of the marker
- m
the number of the fitted model
- model
the type of the fitted model. Possible values are "b1", "b2", "b1,q",
"b2,q", each by itself or followed by "HW" or "pop". The first 4 refer to
the models for the mixture means: b1 and b2 indicate 1 or 2
parameters for signal background, q indicates that a quadratic term in the
signal response was fitted as well. HW and pop refer to the restrictions on
the mixing proportions: HW indicates that the mixing proportions were
constrained according to Hardy-Weinberg equilibrium ratios in case of only
one population, pop indicates that multiple populations were fitted (see
Details section). For more details see Voorrips et al (2011),
doi:10.1186/1471-2105-12-172.
- nsamp
the number of samples for this marker for which select==TRUE,
i.e. the number on which the call rate is based.
- nsel
the number of these samples that have a non-NA ratio value
- npar
the number of free parameters fitted
- iter
the number of iterations to reach convergence
- dip
whether the model had a dip (a smaller peak surrounded by
larger peaks): 0=no, 1=yes
- LL
the log-likelihood of the model
- AIC
Akaike's Information Criterion
- BIC
Bayesian Information Criterion
- selcrit
the selection criterion; the model with the lowest selcrit
is selected. If argument sd.target is NA selcrit is equal to BIC, else
selcrit is larger than BIC if the standard deviation of the mixture
components is larger than sd.target; see Details for details.
- minsepar
a measure of the minimum peak separation. Each difference
of the means of two successive mixture components is divided by the average
of the standard deviations of the two components. The minimum of the
values is reported. All calculations are on the arcsine-square root
transformed scale.
- meanP
For each sample the maximum probability of belonging to any
mixture component is calculated. The average of these P values is reported
in meanP
- P80 .. P99
the fraction of samples that have a probability of at
least 0.80 .. 0.99 to belong to one of the mixture components (by default a
level of 0.99 is required to assign a genotype score to a sample)
- muact0 ..
the actual means of the samples in each of
the mixture components for dosages 0 .. <ploidy> on the transformed scale
- sdact0 ..
the actual standard deviations of the samples in each of
the mixture components for dosages 0 .. <ploidy> on the transformed scale
- mutrans0 ..
the means of the mixture components for dosages 0 ..
<ploidy> on the transformed scale
- sdtrans0 ..
the standard deviations of the mixture components for
dosages 0 .. <ploidy> on the transformed scale
- P0 ..
the mixing proportions of the mixture components for dosages
0 to <ploidy>. If multiple populations are specified there are two
possibilities: (1) the specified population structure is used in the
current model; then for each population the mixing proportions are given
as <npop> sequences of <ploidy+1> fractions, or (2) the population
structure is ignored for the current model, the mixing proportions are given
in the first sequence of <ploidy+1> fractions and all following sequences
are filled with NA. The the item names are adapted to have the
population names between the P and the dosage
- mu0 ..
the model means of the <ploidy+1> mixture components
back-transformed to the original scale
- sd0 ..
the model standard deviations of the <ploidy+1> mixture
components back-transformed to the original scale
- message
if no model was fitted or the model was rejected, the reason
is reported here
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# These examples run for a total of about 9 sec.
data(fitPoly_data)
# triploid, no specified populations
fp <- fitOneMarker(ploidy=3, marker="mrk039",
data=fitPoly_data$ploidy3$dat3x)
# tetraploid, specified populations
# plot of the fitted model saved in tempdir()
fp <- fitOneMarker(ploidy=4, marker=2,
data=fitPoly_data$ploidy4$dat4x,
population=fitPoly_data$ploidy4$pop4x,
pop.parents=fitPoly_data$ploidy4$pop.par4x,
plot="fitted",
plot.dir=paste0(tempdir(),"/fpPlots4x"))
# hexaploid, specified populations, start values for means,
# plot of the fitted model saved in tempdir()
fp <- fitOneMarker(ploidy=6, marker=1,
data=fitPoly_data$ploidy6$dat6x,
population=fitPoly_data$ploidy6$pop6x,
pop.parents=fitPoly_data$ploidy6$pop.par6x,
startmeans=fitPoly_data$ploidy6$startmeans6x,
plot="fitted", plot.dir=paste0(tempdir(),"/fpPlots6x"))
fitPoly documentation built on May 2, 2019, 6:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Examples
Description
This function takes a data frame with allele signal ratios for multiple bi-allelic markers and samples, and fits multiple mixture models to a selected marker. It returns a list, reporting on the performance of these models, selecting the best one based on the BIC criterion, optionally plotting results.
Usage
1 2 3 4 5 6 | fitOneMarker(ploidy, marker, data, diplo=NULL, select=TRUE,
diploselect=TRUE, pop.parents=NULL, population=NULL, parentalPriors=NULL,
samplePriors=NULL, startmeans=NULL, maxiter=40, maxn.bin=200, nbin=200,
sd.threshold=0.1, p.threshold=0.99, call.threshold=0.6, peak.threshold=0.85,
try.HW=TRUE, dip.filter=1, sd.target=NA,
plot="none", plot.type="png", plot.dir, sMMinfo=NULL)
|
Arguments
ploidy |
The ploidy level, 2 or higher: 2 for diploids, 3 for triploids etc. |
marker |
A marker name of number. Used to select the data for one marker, referring to the MarkerName column of parameter data. If a number, the number of the marker based on alphabetic order of the MarkerNames in data. |
data |
A data frame with the polyploid samples, with (at least) columns MarkerName, SampleName and ratio, where ratio is the Y-allele signal divided by the sum of the X- and Y-allele signals: ratio == Y/(X+Y) |
diplo |
NULL or a data frame like data, with the diploid samples and (a subset
of) the same markers as in data. Genotypic scores for diploid samples are
calculated according to the best-fitting model calculated for the polyploid
samples and therefore may range from 0 (nulliplex) to <ploidy>, with the
expected dosages 0 and <ploidy> for the homozygotes and <ploidy/2> for the
heterozygotes. |
select |
A logical vector, recycled if shorter than nrow(data): indicates which rows of data are to be used (default TRUE, i.e. keep all rows) |
diploselect |
A logical vector like select, matching diplo instead of data |
pop.parents |
NULL or a data.frame specifying the population structure. The data frame has 3 columns: the first containing population IDs, the 2nd and 3rd with the population IDs of the parents of these populations (if F1's) or NA (if not). The poopulation IDs should match those in parameter population. If pop.parents is NULL all samples are considered to be in one population, and parameter population should also be NULL (default). |
population |
NULL or a data.frame specifying to which population each sample belongs. The data frame has two columns, the first containing the SampleName (containing all SampleNames occurring in data), the second column containing population IDs that match pop.parents. In both columns NA values are not allowed. Parameters pop.parents and population should both be NULL (default) or both be specified. |
parentalPriors |
NULL or a data frame specifying the prior dosages for
the parental populations. The data frame has one column MarkerName
followed by one column for each F1 parental population. Column names (except
first) are population IDs matching the parental populations in pop.parents.
In case there is just one F1 population in pop.parents, it is possible to
have two columns for both parental populations instead of one (allowing two
specify two different prior dosages); in that case both columns for each
parent have the same caption. Each row specifies the priors for
one marker. The contents of the data frame are dosages, as integers from 0
to <ploidy>; NA values are allowed. |
samplePriors |
NULL or a data.frame specifying prior dosages for individual
samples. The first column called MarkerName is followed by one column per
sample; not all samples in data need to have a column here, only
those samples for which prior dosages for one or more markers are available.
Each row specifies the priors for one marker. The contents of the data frame
are dosages, as integers from 0 to <ploidy>; NA values are allowed. |
startmeans |
NULL or a data.frame specifying the prior means of the mixture distributions. The data frame has one column MarkerName, followed by <ploidy+1> columns with the prior means on the original (untransformed) scale. Each row specifies the means for one marker in strictly ascending order (all means NA is allowed, but markers without start means can also be omitted). |
maxiter |
A single integer, passed to CodomMarker, see there for explanation |
maxn.bin |
A single integer, passed to CodomMarker, see there for explanation |
nbin |
A single integer, passed to CodomMarker, see there for explanation |
sd.threshold |
The maximum value allowed for the (constant) standard deviation of each peak on the arcsine - square root transformed scale, default 0.1. If the optimal model has a larger standard deviation the marker is rejected. Set to a large value (e.g. 1) to disable this filter. |
p.threshold |
The minimum P-value required to assign a genotype (dosage) to a sample; default 0.99. If the P-value for all possible genotypes is less than p.threshold the sample is assigned genotype NA. Set to 1 to disable this filter. |
call.threshold |
The minimum fraction of samples to have genotypes assigned ("called"); default 0.6. If under the optimal model the fraction of "called" samples is less than call.threshold the marker is rejected. Set to 0 to disable this filter. |
peak.threshold |
The maximum allowed fraction of the scored samples that are in one peak; default 0.85. If any of the possible genotypes (peaks in the ratio histogram) contains more than peak.threshold of the samples the marker is rejected (because the remaining samples offers too little information for reliable model fitting). |
try.HW |
Logical: if TRUE (default), try models with and without a constraint on the mixing proportions according to Hardy-Weinberg equilibrium ratios. If FALSE, only try models without this constraint. Even when the HW assumption is not applicable, setting try.HW to TRUE often still leads to a better model. For more details on how try.HW is used see the Details section. |
dip.filter |
if 1 (default), select best model only from models that do not have a dip (a lower peak surrounded by higher peaks: these are not expected under Hardy-Weinberg equilibrium or in cross progenies). If all fitted models have a dip still the best of these is selected. If 2, similar, but if all fitted models have a dip the marker is rejected. If 0, select best model among all fitted models, including those with a dip. |
sd.target |
If the fitted standard deviation of the peaks on the transformed scale is larger than sd.target a penalty is given (see Details); default NA i.e. no penalty is given. |
plot |
String, "none" (default), "fitted" or "all". If "fitted" a plot of the best fitting model and the assigned genotypes is saved with filename <marker number><marker name>.<plot.type>, preceded by "rejected_" if the marker was rejected. If "all", small plots of all models are saved to files (8 per file) with filename <"plots"><marker number><A..F><marker name>.<plot.type> in addition to the plot of the best fitting model. |
plot.type |
String, "png" (default), "emf", "svg" or "pdf". Indicates format for saving the plots. |
plot.dir |
String, the directory where to save the plot files. Must be specified if plot is not "none". Set this to "" to save plot files in the current working directory. |
sMMinfo |
NULL (default), for internal use only. Prevents unneeded checking and recalculation of input parameters when called from saveMarkerModels. |
Details
fitOneMarker fits a series of mixture models for the given marker by repeatedly calling CodomMarker and selects the optimal one. The initial models vary according to the values of try.HW, pop.parents, parentalPriors, samplePriors and startmeans:
no pop.parents, try.HW FALSE: 4 models with different constraints on the means (different or equal X and Y background signal, ratio a linear or quadratic function of dosage), no restrictions on the mixing proportions (the fractions of samples in each dosage peak)
no pop.parents, try.HW TRUE: The previous 4 models are fitted and also 4 models with the same restrictions on the means and the mixing proportions restricted to Hardy-Weinberg ratios (assuming polysomic inheritance)
pop.parents specified, no parentalPriors / samplePriors / startmeans, try.HW FALSE: 4 models are fitted with the same restrictions on the means as above, but with different restrictions on the mixing proportions for each population: no restriction on parental populations, none on accession panels, polysomic F1 segregation ratios on F1 populations. Additionally 4 models are fitted with all samples considered as one population, with the same 4 models for the means and no restrictions on mixing proportions.
pop.parents specified, no parentalPriors / samplePriors / startmeans, try.HW TRUE: 4 models are fitted with the same restrictions on the means as above, but with different restrictions on the mixing proportions for each population: no restriction on parental populations, HW-ratios for accession panels, polysomic F1 segregation ratios on F1 populations. Additionally 4 models are fitted with all samples considered as one population, with the same 4 models for the means and mixing proportions according to HW ratios.
pop.parents and parentalPriors specified, try.HW FALSE: 4 models are fitted with the same restrictions on the means as above, but with different restrictions on the mixing proportions for each population: no restriction on parental populations and the accession panels, polysomic F1 segregation ratios on F1 populations ignoring the parental priors. Additionally 4 models are fitted with the same restrictions on the means and mixing proportions of the accession panels, but where the mixing proportions of the parental populations are set to (almost) 1 for the prior dosage and (almost) 0 for all other dosages, and those for the F1 populations to the polysomic segregation ratios expected for the parental priors.
pop.parents and parentalPriors specified, try.HW TRUE: same as with try.HW FALSE, except that the mixing proportions of accession panels are now restricted to HW ratios.
if parentalPriors and/or samplePriors are specified, these and the signal ratios of the corresponding samples are (also) used to estimate starting values of the mixture component means in the EM algorithm. Alternatively startmeans can be specified directly.
Because convergence to the optimal solution often fails, the models are fitted
with several start values for the <ploidy+1> means of the mixture
distributions: (1) based on initial clustering of the ratios, (2) based on
a uniform distribition from 0.02 to pi/2-0.02 on the asin(sqrt(x)) scale,
and (3) if startmeans are specified or can be calculated from samplePriors
and/or parentalPriors these are used for a third set of model fits.
The main difference between parentalPriors and samplePriors is that
parentalPriors are treated as fixed (and if both parents of an F1 population
have priors, the F1 segregation is also fixed) while samplePriors are only
used to calculate starting ratio means for each dosage. Depending on the
confidence the user has in the prior dosages of the parents they can be
supplied as parentalPriors or samplePriors.
In some cases an additional fit is performed with a modified set of initial
means.
An optimal model is selected based on the Bayesian Information Criterium
(BIC), which takes into account the Log-Likelihood and the number of fitted
parameters of the models. If sd.target is specified and the standard
deviation of the mixture model components is larger than this target a
penalty is applied, making is less likely that that model is selected.
The plots consist of one histogram per (non-parent) population showing the
frequency distribution of the signal ratios of the samples in that population.
The fitted model is shown in green (density and means), and for F1 populations
the samples of parent 1 and 2 are shown as red and blue triangles.
If diplodata are present, a histogram for the diploid samples is plotted
in the top histogram (diploid bars are narrower and gray). The diploid bars
are scaled so the maximum bar is half the maximum polyploid bar. At the
bottom of the plot for the fitted model a rug plot shows the scores of each
sample, while the bottom (red) samples are unscored.
Value
A list with components:
- log
A character vector with the lines of the log text.
- modeldata
A data frame as allmodeldata (see below) with only the one row with data on the selected model.
- allmodeldata
A data frame with for each tried model one row with the marker number, marker name, number of samples and (if the marker is not rejected) data of the fitted model (see below).
- scores
A data frame with the name and data for all samples (including NA's for the samples that were not selected, see parameter select), with columns:
marker (the sequential number of the marker (based on alphabetic order of the marker names in data)
MarkerName
SampleName
ratio (the given ratio from parameter data)
P0 .. P<ploidy> (the probabilities that this sample belongs to each of the <ploidy+1> mixture components)
maxgeno (0..ploidy, the genotype = mixture component with the highest P value)
maxP (the P value for this genotype)
geno (the assigned genotype number: same as maxgeno, or NA if maxP < p.threshold).- diploscores
A data frame like scores for the samples in the data frame supplied with argument diplo. If diplo is NA also diploscores will be NA.
The modeldata and allmodeldata data frames present data on a fitted model. modeldata presents data on the selected model; allmodeldata lists all attempted models. Both data frames contain the following columns:
- marker
the sequential number of the marker (based on alphabetic order of the marker names in data)
- MarkerName
the name of the marker
- m
the number of the fitted model
- model
the type of the fitted model. Possible values are "b1", "b2", "b1,q", "b2,q", each by itself or followed by "HW" or "pop". The first 4 refer to the models for the mixture means: b1 and b2 indicate 1 or 2 parameters for signal background, q indicates that a quadratic term in the signal response was fitted as well. HW and pop refer to the restrictions on the mixing proportions: HW indicates that the mixing proportions were constrained according to Hardy-Weinberg equilibrium ratios in case of only one population, pop indicates that multiple populations were fitted (see Details section). For more details see Voorrips et al (2011), doi:10.1186/1471-2105-12-172.
- nsamp
the number of samples for this marker for which select==TRUE, i.e. the number on which the call rate is based.
- nsel
the number of these samples that have a non-NA ratio value
- npar
the number of free parameters fitted
- iter
the number of iterations to reach convergence
- dip
whether the model had a dip (a smaller peak surrounded by larger peaks): 0=no, 1=yes
- LL
the log-likelihood of the model
- AIC
Akaike's Information Criterion
- BIC
Bayesian Information Criterion
- selcrit
the selection criterion; the model with the lowest selcrit is selected. If argument sd.target is NA selcrit is equal to BIC, else selcrit is larger than BIC if the standard deviation of the mixture components is larger than sd.target; see Details for details.
- minsepar
a measure of the minimum peak separation. Each difference of the means of two successive mixture components is divided by the average of the standard deviations of the two components. The minimum of the values is reported. All calculations are on the arcsine-square root transformed scale.
- meanP
For each sample the maximum probability of belonging to any mixture component is calculated. The average of these P values is reported in meanP
- P80 .. P99
the fraction of samples that have a probability of at least 0.80 .. 0.99 to belong to one of the mixture components (by default a level of 0.99 is required to assign a genotype score to a sample)
- muact0 ..
the actual means of the samples in each of the mixture components for dosages 0 .. <ploidy> on the transformed scale
- sdact0 ..
the actual standard deviations of the samples in each of the mixture components for dosages 0 .. <ploidy> on the transformed scale
- mutrans0 ..
the means of the mixture components for dosages 0 .. <ploidy> on the transformed scale
- sdtrans0 ..
the standard deviations of the mixture components for dosages 0 .. <ploidy> on the transformed scale
- P0 ..
the mixing proportions of the mixture components for dosages 0 to <ploidy>. If multiple populations are specified there are two possibilities: (1) the specified population structure is used in the current model; then for each population the mixing proportions are given as <npop> sequences of <ploidy+1> fractions, or (2) the population structure is ignored for the current model, the mixing proportions are given in the first sequence of <ploidy+1> fractions and all following sequences are filled with NA. The the item names are adapted to have the population names between the P and the dosage
- mu0 ..
the model means of the <ploidy+1> mixture components back-transformed to the original scale
- sd0 ..
the model standard deviations of the <ploidy+1> mixture components back-transformed to the original scale
- message
if no model was fitted or the model was rejected, the reason is reported here
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | # These examples run for a total of about 9 sec.
data(fitPoly_data)
# triploid, no specified populations
fp <- fitOneMarker(ploidy=3, marker="mrk039",
data=fitPoly_data$ploidy3$dat3x)
# tetraploid, specified populations
# plot of the fitted model saved in tempdir()
fp <- fitOneMarker(ploidy=4, marker=2,
data=fitPoly_data$ploidy4$dat4x,
population=fitPoly_data$ploidy4$pop4x,
pop.parents=fitPoly_data$ploidy4$pop.par4x,
plot="fitted",
plot.dir=paste0(tempdir(),"/fpPlots4x"))
# hexaploid, specified populations, start values for means,
# plot of the fitted model saved in tempdir()
fp <- fitOneMarker(ploidy=6, marker=1,
data=fitPoly_data$ploidy6$dat6x,
population=fitPoly_data$ploidy6$pop6x,
pop.parents=fitPoly_data$ploidy6$pop.par6x,
startmeans=fitPoly_data$ploidy6$startmeans6x,
plot="fitted", plot.dir=paste0(tempdir(),"/fpPlots6x"))
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.