Description Usage Arguments Details Author(s) Examples
It is used to perform genome-wide association studies on individuals that are both related and have repeated measurements. The function computes score statistic based p-values for a linear mixed model including random polygenic effects and a random effect for repeated measurements. A p-value is computed for each marker and the null hypothesis tested is a zero additive marker effect.
1 2 |
formula.FixedEffects |
Formula including the response variable and cofactors as fixed effects. |
genabel.data |
An GenABEL object including marker information. This object has one observtion per individuals. |
phenotype.data |
A data frame including the repeated observations and IDs. |
id.name |
The column name of the IDs in phen.data |
GRM |
An optional genetic relationship matrix (GRM) can be included as input. Otherwise the GRM is computed within the function. |
V |
An optional (co)variance matrix can be included as input. Otherwise it is computed using the hglm function. |
memory |
Used to optimize computations. The maximum number of elements in a matrix that can be stored efficiently. |
A generalized squares (GLS) is fitted for each marker given a (co)variance matrix V.
The computations are made fast by transforming the GLS to
an ordinary least-squares (OLS) problem using an eigen-decomposition of V.
The OLS are computed using QR-factorization. If V is not specified then a model
including random polygenic effects and permanent environmental effects is
fitted (using the hglm package) to compute V. A GenABEL object (scan.gwaa class)
is returned (including also the hglm
results).
Let e.g. GWAS1 be an object returned by the rGLS
function.
Then a Manhattan plot can be produced by calling plot(GWAS1)
and
the top SNPs using summary(GWAS1)
. Both of these functions are
generic GenABEL functions.
The results from the fitted linear mixed model without any SNP effect included
are produced by calling summary(GWAS1@call$hglm)
.
Lars Ronnegard
1 2 3 4 5 6 7 | data(Phen.Data) #Phenotype data with repeated observations
data(gen.data) #GenABEL object including IDs and marker genotypes
GWAS1 <- rGLS(y ~ age + sex, genabel.data = gen.data, phenotype.data = Phen.Data)
plot(GWAS1, main="")
summary(GWAS1)
#Summary for variance component estimation without SNP effects
summary(GWAS1@call$hglm)
|
Loading required package: hglm
Loading required package: Matrix
Loading required package: MASS
Loading required package: hglm.data
hglm: Hierarchical Generalized Linear Models
Version 2.1-1 (2015-08-28) installed
Authors: Moudud Alam, Lars Ronnegard, Xia Shen
Maintainer: Xia Shen <xia.shen@ki.se>
Use citation("hglm") to know how to cite our work.
Discussion: https://r-forge.r-project.org/forum/?group_id=558
BugReports: https://r-forge.r-project.org/tracker/?group_id=558
VideoTutorials: http://www.youtube.com/playlist?list=PLn1OmZECD-n15vnYzvJDy5GxjNpVV5Jr8
Loading required package: GenABEL
Loading required package: GenABEL.data
[1] "GRM ready"
[1] "Variance component estimation ready"
[1] "Rotation matrix ready"
[1] "Rotate LMM started"
[1] "Rotate LMM ready"
Summary for top 10 results, sorted by P1df
Chromosome Position Strand A1 A2 effB se_effB chi2.1df
rs120315 1 3725352 + T C 2.1696316 0.4258979 NA
rs9670687 1 4779800 - C A -2.0412127 0.5351582 NA
rs891586 2 8024318 + A G -1.6669267 0.4478499 NA
rs1352451 2 8022651 - A C 1.6169429 0.4564078 NA
rs1252282 2 8167984 - A T -0.9982089 0.2825403 NA
rs7633966 1 3721952 + C T 1.8586037 0.5325691 NA
rs9922492 1 3729983 + C T -1.8586037 0.5325691 NA
rs4378234 1 4800348 - G T -1.9012073 0.5509388 NA
rs7499832 3 10242953 - G C -1.0897248 0.3182487 NA
rs6561272 1 3715538 - T G -1.9073305 0.5726170 NA
P1df Pc1df
rs120315 3.501204e-07 NA
rs9670687 1.366120e-04 NA
rs891586 1.975997e-04 NA
rs1352451 3.959643e-04 NA
rs1252282 4.109058e-04 NA
rs7633966 4.832322e-04 NA
rs9922492 4.832322e-04 NA
rs4378234 5.588239e-04 NA
rs7499832 6.167705e-04 NA
rs6561272 8.656541e-04 NA
Call:
hglm.default(X = X, y = y, Z = Z, maxit = 200, RandC = c(ncol(Z.GRM),
ncol(Z.indx)))
----------
MEAN MODEL
----------
Summary of the fixed effects estimates:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 102.05662 0.78674 129.721 < 2e-16 ***
age 0.06504 0.01451 4.483 1.06e-05 ***
sex 0.85160 0.38016 2.240 0.0258 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: P-values are based on 290 degrees of freedom
Summary of the random effects estimates:
Estimate Std. Error
Z.1 -2.1532 0.7458
Z.2 0.0267 0.7619
Z.3 -0.2035 0.7559
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
Summary of the random effects estimates:
Estimate Std. Error
Z.101 0.3941 0.9913
Z.102 0.0778 0.8626
Z.103 1.5205 0.8765
...
NOTE: to show all the random effects, use print(summary(hglm.object), print.ranef = TRUE).
----------------
DISPERSION MODEL
----------------
NOTE: h-likelihood estimates through EQL can be biased.
Dispersion parameter for the mean model:
[1] 3.090352
Model estimates for the dispersion term:
Link = log
Effects:
Estimate Std. Error
1.1283 0.0830
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
Dispersion parameter for the random effects:
[1] 1.376 1.354
Dispersion model for the random effects:
Link = log
Effects:
.|Random1
Estimate Std. Error
0.3194 0.2495
.|Random2
Estimate Std. Error
0.3032 0.2407
Dispersion = 1 is used in Gamma model on deviances to calculate the standard error(s).
EQL estimation converged in 2 iterations.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.