[Genabel-commits] r2012 - in pkg/MultiABEL: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Aug 15 18:32:37 CEST 2015
Author: shenxia
Date: 2015-08-15 18:32:37 +0200 (Sat, 15 Aug 2015)
New Revision: 2012
Modified:
pkg/MultiABEL/ChangeLog
pkg/MultiABEL/DESCRIPTION
pkg/MultiABEL/NAMESPACE
pkg/MultiABEL/R/MultiLoad.R
pkg/MultiABEL/R/MultiSummary.R
pkg/MultiABEL/R/Multivariate.R
pkg/MultiABEL/R/load.summary.R
pkg/MultiABEL/src/MultiSummaryLoop.f90
Log:
* Updated to version 1.1-3
* load.summary() case not sensitive anymore for file headers
* load.summary() stringsAsFactors bug fixed
* MultiSummary() coefficients standard errors provided for outbred HWE populations
Modified: pkg/MultiABEL/ChangeLog
===================================================================
--- pkg/MultiABEL/ChangeLog 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/ChangeLog 2015-08-15 16:32:37 UTC (rev 2012)
@@ -1,3 +1,12 @@
+2015-08-15 18:28 xia
+
+ * Updated to version 1.1-3
+ * load.summary() case not sensitive anymore
+ for file headers
+ * load.summary() stringsAsFactors bug fixed
+ * MultiSummary() coefficients standard errors
+ provided for outbred HWE populations
+
2015-07-31 15:54 xia
* Updated to version 1.1-2
Modified: pkg/MultiABEL/DESCRIPTION
===================================================================
--- pkg/MultiABEL/DESCRIPTION 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/DESCRIPTION 2015-08-15 16:32:37 UTC (rev 2012)
@@ -1,8 +1,8 @@
Package: MultiABEL
Type: Package
Title: Multi-Trait Genome-Wide Association Analyses
-Version: 1.1-2
-Date: 2015-07-31
+Version: 1.1-3
+Date: 2015-08-15
Author: Xia Shen
Maintainer: Xia Shen <xia.shen at ki.se>
Description: Multivariate genome-wide association analyses. The analysis can be performed on individual-level data or multiple single-trait genome-wide summary statistics.
@@ -14,4 +14,4 @@
DatABEL
License: GPL (>= 2)
LazyLoad: yes
-Packaged: 2015-07-31 16:03:18 CET; xia
+Packaged: 2015-08-15 18:30:18 CET; xia
Modified: pkg/MultiABEL/NAMESPACE
===================================================================
--- pkg/MultiABEL/NAMESPACE 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/NAMESPACE 2015-08-15 16:32:37 UTC (rev 2012)
@@ -1,4 +1,4 @@
useDynLib(MultiABEL)
import("svMisc")
-export("Multivariate", "MultiRep", "MultiMeta", "MultiLoad","load.summary", "MultiSummary")
+export("Multivariate", "MultiRep", "MultiMeta", "MultiLoad","load.summary", "MultiSummary", "StratSummary")
exportClasses("multi.summary", "multi.loaded")
\ No newline at end of file
Modified: pkg/MultiABEL/R/MultiLoad.R
===================================================================
--- pkg/MultiABEL/R/MultiLoad.R 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/R/MultiLoad.R 2015-08-15 16:32:37 UTC (rev 2012)
@@ -56,7 +56,8 @@
#' @keywords multivariate, multiload, load
#'
`MultiLoad` <- function(gwaa.data = NULL, phenofile = NULL, genofile = NULL,
- trait.cols, covariate.cols = NULL, cuts = 20, impute = TRUE, ...) {
+ trait.cols, covariate.cols = NULL, cuts = 20,
+ impute = TRUE, gaussianize = TRUE, ...) {
set.seed(911)
if (!require(GenABEL) | !require(DatABEL)) {
stop('GenABEL and DatABEL packages required!')
@@ -89,7 +90,11 @@
okidx <- which(!is.na(rowSums(Y)))
Y <- Y[okidx,]
n <- length(okidx)
- for (i in 1:ncol(Y)) Y[,i] <- zscore(Y[,i])
+ if (gaussianize) {
+ for (i in 1:ncol(Y)) {
+ Y[,i] <- zscore(Y[,i])
+ }
+ }
X <- matrix(1, length(okidx), 1)
if (!is.null(covariate.cols)) {
X <- cbind(X, as.matrix(pheno[okidx,covariate.cols]))
Modified: pkg/MultiABEL/R/MultiSummary.R
===================================================================
--- pkg/MultiABEL/R/MultiSummary.R 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/R/MultiSummary.R 2015-08-15 16:32:37 UTC (rev 2012)
@@ -65,12 +65,13 @@
n <- x$gwa[,2*m + 2]
k <- nrow(x$gwa)
betamat <- x$gwa[,seq(1, 2*m, 2)]
+ dimnames(betamat) <- list(NULL, NULL)
res <- data.frame(Marker = rownames(x$gwa), Freq = f, N = n, Beta.S = NA, SE = NA, P = NA)
rownames(res) <- rownames(x$gwa)
if (type == 'outbred') {
scan <- .Fortran('MultiSummaryLoop', k = as.integer(k), m = as.integer(m), nn = n,
f = f, betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno),
- E0 = diag(m), H = diag(m), coef = betamat,
+ E0 = diag(m), H = diag(m), invE0 = diag(m), s = betamat, coef = betamat,
pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
} else if (type == 'inbred') {
scan <- .Fortran('MultiSummaryLoopInbred', k = as.integer(k), m = as.integer(m), nn = n,
@@ -89,7 +90,9 @@
res$Beta.S <- scan$pil
res$P <- pf(scan$fstat, m, n - m - 1, lower.tail = FALSE)
res$SE <- res$Beta.S/sqrt(qchisq(res$P, 1, lower.tail = FALSE))
- coef <- scan$coef
- colnames(coef) <- paste('Coef.Trait', 1:m, sep = '.')
+ coef <- matrix(NA, k, m*2)
+ coef[,seq(1, 2*m, 2)] <- scan$coef
+ coef[,seq(2, 2*m, 2)] <- scan$s
+ colnames(coef) <- colnames(x$gwa)[1:(2*m)]
return(data.frame(res, coef))
}
Modified: pkg/MultiABEL/R/Multivariate.R
===================================================================
--- pkg/MultiABEL/R/Multivariate.R 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/R/Multivariate.R 2015-08-15 16:32:37 UTC (rev 2012)
@@ -5,6 +5,7 @@
#' analysis of variance (MANOVA).
#'
#' @param x An object created by \code{\link{MultiLoad}}.
+#' @param trait.idx A vector giving the indices of traits to be analyzed.
#' @param ... not used.
#'
#' @note Either \code{gwaa.data} (for GenABEL data format) or the combination of
@@ -34,7 +35,7 @@
#' data(ge03d2ex.clean)
#'
#' ## running multivariate GWAS for 3 traits: height, weight, bmi
-#' loaded <- Multivariate(gwaa.data = ge03d2ex.clean, trait.cols = c(5, 6, 8),
+#' loaded <- MultiLoad(gwaa.data = ge03d2ex.clean, trait.cols = c(5, 6, 8),
#' covariate.cols = c(2, 3))
#'
#' ## running the multivariate GWAS again
@@ -43,10 +44,22 @@
#' @aliases Multivariate, multivariate
#' @keywords multivariate
#'
-`Multivariate` <- function(x, ...) {
+`Multivariate` <- function(x, trait.idx = NULL, ...) {
if (class(x) != 'multi.loaded') {
stop('wrong data type!')
}
+ if (!is.null(trait.idx)) {
+ idx <- as.numeric(t(cbind(trait.idx*2 - 1, trait.idx*2)))
+ x$gwa <- x$gwa[,idx]
+ x$cor.pheno <- x$cor.pheno[trait.idx,trait.idx]
+ }
+ if (any(is.na(rowSums(x$gwa)))) {
+ naidx <- which(is.na(rowSums(x$gwa)))
+ x$gwa <- x$gwa[-naidx,]
+ x$cvg <- x$cvg[-naidx]
+ x$snp <- x$snp[-naidx]
+ x$nsnp <- x$nsnp - length(naidx)
+ }
npheno <- ncol(x$cor.pheno)
cat('initializing analysis ...')
res <- matrix(NA, x$nsnp, 5)
Modified: pkg/MultiABEL/R/load.summary.R
===================================================================
--- pkg/MultiABEL/R/load.summary.R 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/R/load.summary.R 2015-08-15 16:32:37 UTC (rev 2012)
@@ -74,7 +74,7 @@
data <- c()
fn <- rev(files)
for (i in m:1) {
- dd <- read.table(fn[i], header = TRUE)
+ dd <- read.table(fn[i], header = TRUE, stringsAsFactors = FALSE)
colnames(dd) <- tolower(colnames(dd))
idx <- which(duplicated(dd$snp))
if (length(idx) > 0) {
@@ -153,6 +153,9 @@
cor.pheno <- cor(z, use = 'pairwise.complete.obs')
cat('done.\n')
}
+ dimnames(cor.pheno) <- list(files, files)
+ gwanames <- c(paste(rep(files, each = 2), rep(c('.beta', '.se'), m), sep = ''), 'f', 'n')
+ colnames(gwa0) <- gwanames
dd <- list(gwa = gwa0, cor.pheno = cor.pheno)
class(dd) <- 'multi.summary'
return(dd)
Modified: pkg/MultiABEL/src/MultiSummaryLoop.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoop.f90 2015-08-10 13:39:25 UTC (rev 2011)
+++ pkg/MultiABEL/src/MultiSummaryLoop.f90 2015-08-15 16:32:37 UTC (rev 2012)
@@ -4,27 +4,37 @@
! Created by Xia Shen on 5/29/15.
!
-subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, E0, H, coef, pil, fstat)
+subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, E0, H, invE0, s, coef, pil, fstat)
implicit none
integer :: i, j, k, m
- double precision, dimension(m, m) :: R, invR, E0, H, HinvE0
- double precision, dimension(1, m) :: beta
+ double precision, dimension(m, m) :: R, invR, E0, H, HinvE0, D, invE0
+ double precision, dimension(1, m) :: beta, b
double precision, dimension(m, 1) :: tmp
double precision, dimension(k) :: nn, f, pil, fstat
- double precision, dimension(k, m) :: betamat, coef
+ double precision, dimension(k, m) :: betamat, coef, s
double precision, dimension(m) :: lambda
- double precision :: vg
+ double precision, dimension(1, 1) :: bDbeta0
+ double precision :: vg, gg, bDbeta
do j = 1, k
vg = 2*f(j)*(1 - f(j))
E0 = nn(j)*R
+ invE0 = invR/nn(j)
do i = 1, m
beta(1,i) = betamat(j,i)
+ D(i,i) = E0(i,i)
end do
tmp = matmul(invR, transpose(beta))*vg
do i = 1, m
coef(j,i) = tmp(i,1)
+ b(1,i) = tmp(i,1)
end do
+ gg = vg*nn(j)
+ bDbeta0 = matmul(matmul(b, D), transpose(beta))
+ bDbeta = bDbeta0(1,1)
+ do i = 1, m
+ s(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)
+ end do
H = matmul(transpose(beta), beta)*nn(j)*vg
HinvE0 = matmul(H, invR)/nn(j)
do i = 1, m
More information about the Genabel-commits
mailing list