qb.mcmc {qtlbim} | R Documentation |
A computationally efficient MCMC algorithm using the Gibbs sampler or Metropolis-Hastings algorithm is used to produce posterior samples for QTL mapping.
qb.mcmc(cross, data, model, mydir = ".", n.iter = 3000, n.thin = 40, n.burnin = 0.01*n.iter*n.thin, algorithm = c("M-H","Gibbs"), genoupdate = TRUE, seed = 0, verbose = TRUE, ...)
cross |
An object of class cross . See read.cross for details. |
data |
the list returned by calling the function qb.data . |
model |
the list returned by calling the function qb.model . |
mydir |
a directory to save output from qb.mcmc in several ‘*.dat’ files.
A directory is created using the trait name and the system time and date.
If no directory is specified, the default directory is the current working
directory. |
n.iter |
number of iterations to be saved in mydir , the default being
3000. |
n.thin |
the thinning number which must be a positive number (default=40) |
n.burnin |
the initial burn-in period, i.e number of iterations to discard at the beginning of the MCMC run default being 0.01*n.iter*n.thin. |
algorithm |
specifies the sampling algorithm for MCMC : Gibbs sampler ("Gibbs") or Metropolis-Hastings algorithm ("M-H") |
genoupdate |
=TRUE will update QTL genotypes and =FALSE will not do so and use the expected value of the QTL genotypes. |
seed |
Specifies the seed for the random number generator. Using the same seed
for two runs of the qb.mcmc function will generate the exact same
output. The seed needs to be an integer. The default value for seed
is the system time. |
verbose |
=TRUE will force periodic output of the number of MCMC iterations saved. The location of the output directory where results are stored and the time taken for the MCMC run will also be displayed to the user. |
... |
Paramters passed to qb.data or
qb.model if data or model , respectively,
is not provided. |
A composite model space approach to develop a Basian model selection framework for identifying interacting QTL for complex traits in experimental crosses from two inbred lines. By placing a liberal constraint on the upper bound of the number of detectable QTL we restrict attention to models of fixed dimension. Either Gibbs sampler or Metroplis-Hastings algorithm can be applied to sample from the posterior distribution.
The qb.mcmc
function returns a list of class qb
, the components of which
contain input parameters from the cross
object, qb.data
and qb.model
.
Since, the parameters have already been described in their respective man pages
we only include the components which have been added on top of these.
However, the qbOject$subset
component is a way to manipulate the
Monte Carlo samples to make it ready for the high-end plotting routines and might not be
of much importance to the average user.
output.dir |
directory used for saving all outputs generated by qb.mcmc . |
subset |
iterdiag is a vector of integers from 1 to n.iter .
mainloci is a vector of length equaling the sum of the number of
QTLs detected in each iteration. This vector is actually
a permutation vector of the mainloci.dat
file stored in mydir sorted with respect to the
iteration number and ties are broken with the chromosome
number and the locus of the putative QTL.
pairloci is a list of
region a data frame storing the first and last position of the
marker map for each chromosome.
|
The following files are saved in output.dir
:
The iteration file saved in mydir
has n.iter
rows and
4 major columns:
column no 1
: iteration number.
column no 2
: number of putative QTLs included.
column no 3
: the overall mean.
column no 4
: the residual variance.
Depending on the type of cross, presence of covariates and epistatic
effects there would be more columns in the following order:
variance of additive effect, variance of dominant effect, variance of
additive-additive interaction, variance of additive-dominant interaction,
variance of dominant-additive interaction, variance of dominant-dominant
interaction, variance of environment-additive interaction, variance of
environmnet-dominant interaction, variance of environemnt effect, total
genetic variance.
The covariate file saved in mydir
has n.iter
rows and
L+M(length(fixcov)+length(rancov)
) columns:
L columns
: Coeffecient of the fixed effect.
M columns
: Variance of the random effect.
If an ordinal trait is analyzed, the cutoff points for the threshold
model are also included in additional columns. There would be C-3
bounded threshold values for an ordinal phenotype with C categories.
The mainloci file saved in mydir
has the N rows (N=sum of number of QTLs
detected in n.iter
iterations) and 4 major columns:
column no 1
: iteration number.
column no 2
: number of putative QTLs included in the model.
column no 3
: the chromosome number on which a putative QTL has been detected.
column no 4
: the QTL position indicator in the grid.
Depending on the type of cross there would be more columns in the
following order: additive effect, dominant effect, variance of additive effect,
and variance of dominant effect.
The pairloci file saved in mydir
has the N rows (N=sum of number
of pairs of QTLs with epistatic effect detected) and 6 major columns:
column no 1
: iteration number.
column no 2
: number of pairs of QTLs detected to have epistatic effect.
column no 3
: the chromosome number for the first one of each pair.
column no 4
: the QTL position for this one.
column no 5
: the chromosome number for the second one of each pair.
column no 6
: the QTL position for this one.
Depending on the type of cross there would be more columns in the following order:
additive-additive interaction effect, additive-dominant interaction effect,
dominant-additive interaction effect, dominant-dominant interaction effect,
variance of additive-additive interaction, variance of additive-dominant interaction,
variance of dominant-additive interaction, variance of dominant-dominant
interaction.
The gbye (Gene by Environment) file saved in mydir
has 5 major columns:
column no 1
: iteration number.
column no 2
: number of GxE interactions.
column no 3
: fixed covariate number.
column no 4
: chromosome number of the putative QTL in the GxE interaction.
column no 5
: position of the corresponding QTL.
Depending on the type of cross there would be more columns in the following order: additive effect, dominant effect, variance of additive effect, and variance of dominant effect.
Nengjun Yi, nyi@ms.ssg.uab.edu
qb.sim.cross
, qb.data
,
qb.model
, qb.mcmc
## Here is a simple example used in the "scan" vignette. ## Simulation map: 201 markers with 1 cM spacing. sim.map <- list(ch1 = 0:200) names(sim.map$ch1) <- paste("M", 0:200, sep = "") ## Simulation model: QTL at 100 cM of size 1.R sim.model <- rbind(qtl.1=c(chr=1,pos=100,effect=1)) ## Set random number seed. set.seed(1234) ## Simulate BC using R/qtl's sim.cross sim <- sim.cross(map = sim.map, model = sim.model, n.ind = 100, type = "bc") ## Set up genotype probabilities. sim <- qb.genoprob(sim,step=2) ## Create MCMC samples. Vignette had n.iter = 3000. qbSim <- qb.mcmc(sim, epistasis=FALSE, n.iter = 300, seed = 1616) ## Summary of qb object "qbSim". summary(qbSim) ## Remove MCMC samples. qb.remove(qbSim) ## Remove simulated data. remove(sim, sim.model, sim.map) ## Not run: ## Here is the example as save in the R/qtlbim package. ## Simulate F2 cross. cross <- qb.sim.cross(len = rep(100,20), n.mar = 11, eq.spacing =FALSE, n.ind = 500, type = "f2", ordinal = c(0.3,0.3,0.2,0.2), missing.geno = 0.03, missing.pheno = 0.03, qtl.pos = rbind(qtl.1=c(chr=1,pos=15),qtl.2=c(1,45),qtl.3=c(3,12), qtl.4=c(5,15),qtl.5=c(7,15),qtl.6=c(10,15),qtl.7=c(12,35),qtl.8=c(19,15)), qtl.main = rbind(main.1=c(qtl=1,add=0.5,dom=0),main.2=c(2,0,0.7), main3=c(3,-0.5,0),main4=c(4,0.5,-0.5)), qtl.epis = rbind(epis1=c(qtl.a=4,qtl.b=5,aa=-0.7,ad=0,da=0,dd=0), epis2=c(6,8,0,1.2,0,0)), covariate = c(fix.cov=0.5,ran.cov=0.07), gbye = rbind(GxE.1=c(qtl=7,add=0.8,dom=0)) ) ## Calculate grids and genotypic probabilites. cross <- qb.genoprob(cross, step=2) ## Create MCMC samples ## First line as qb.data options; second line has qb.model options. qbExample <- qb.mcmc(cross, pheno.col = 3, rancov = 2, fixcov = 1, chr.nqtl = rep(3,nchr(cross)), prop = c(0.5, 1, 0.5), intcov = 1, n.iter = 3000, n.thin = 10, n.burnin = 1000, genoupdate = FALSE) ## End(Not run)