[Genabel-commits] r2014 - in pkg/MultiABEL: . R src

[noreply at r-forge.r-project.org]

Author: shenxia
Date: 2015-08-20 09:11:59 +0200 (Thu, 20 Aug 2015)
New Revision: 2014

Modified:
   pkg/MultiABEL/ChangeLog
   pkg/MultiABEL/R/MultiSummary.R
   pkg/MultiABEL/R/load.summary.R
   pkg/MultiABEL/src/MultiSummaryLoop.f90
Log:
Reverse traits order bug fixed for load.summary. MultiSummary can selected a part of loaded traits for analysis.

Modified: pkg/MultiABEL/ChangeLog
===================================================================
--- pkg/MultiABEL/ChangeLog	2015-08-18 10:09:34 UTC (rev 2013)
+++ pkg/MultiABEL/ChangeLog	2015-08-20 07:11:59 UTC (rev 2014)
@@ -1,11 +1,14 @@
-2015-08-15 18:28 xia
+2015-08-20 09:09 xia
 
     * Updated to version 1.1-3
     * load.summary() case not sensitive anymore
       for file headers
-    * load.summary() stringsAsFactors bug fixed
+    * load.summary() stringsAsFactors bug fixed;
+      Reversed order of traits bug fixed
     * MultiSummary() coefficients standard errors
-      provided for outbred HWE populations
+      provided for outbred HWE populations;
+      Traits can be partially selected for multi-
+      trait scans.
 
 2015-07-31 15:54 xia
 

Modified: pkg/MultiABEL/R/MultiSummary.R
===================================================================
--- pkg/MultiABEL/R/MultiSummary.R	2015-08-18 10:09:34 UTC (rev 2013)
+++ pkg/MultiABEL/R/MultiSummary.R	2015-08-20 07:11:59 UTC (rev 2014)
@@ -55,11 +55,17 @@
 #' @aliases MultiSummary, multi.summary
 #' @keywords multivariate, meta-analysis
 #' 
-`MultiSummary` <- function(x, type = 'outbred', vars = NULL) {
+`MultiSummary` <- function(x, index = NULL,  type = 'outbred', vars = NULL) {
 	if (class(x) != 'multi.summary') {
 		stop('wrong data type!')
 	}
 	cat('multi-trait genome scan ... ')
+	if (!is.null(index)) {
+		m <- nrow(x$cor.pheno)
+		x$cor.pheno <- x$cor.pheno[index,index]
+		idx <- which(!(1:m %in% index))
+		x$gwa <- x$gwa[,-c(idx*2 - 1, idx*2)]
+	}
 	m <- nrow(x$cor.pheno)
 	f <- x$gwa[,2*m + 1]
 	n <- x$gwa[,2*m + 2]
@@ -69,9 +75,9 @@
 	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), invE0 = diag(m), s = betamat, coef = betamat, 
+		scan <- .Fortran('MultiSummaryLoop', k = as.integer(k), m = as.integer(m), nn = as.numeric(n), 
+			f = as.numeric(f), betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno), 
+			E0 = diag(m), H = diag(m), invE0 = diag(m), D = 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, 
@@ -92,7 +98,7 @@
 	res$SE <- res$Beta.S/sqrt(qchisq(res$P, 1, lower.tail = FALSE))
 	coef <- matrix(NA, k, m*2)
 	coef[,seq(1, 2*m, 2)] <- scan$coef
-	coef[,seq(2, 2*m, 2)] <- scan$s
+	coef[,seq(2, 2*m, 2)] <- sqrt(scan$s)
 	colnames(coef) <- colnames(x$gwa)[1:(2*m)]
 	return(data.frame(res, coef))
 }

Modified: pkg/MultiABEL/R/load.summary.R
===================================================================
--- pkg/MultiABEL/R/load.summary.R	2015-08-18 10:09:34 UTC (rev 2013)
+++ pkg/MultiABEL/R/load.summary.R	2015-08-20 07:11:59 UTC (rev 2014)
@@ -64,6 +64,9 @@
 	if (sum(file.exists(files)) < 2) {
 		stop('number of traits has to be more than 2!')
 	}
+	if (is.null(cor.pheno) & is.null(indep.snps)) {
+		stop('indep.snps required for estimating cor.pheno!')
+	}
 	m <- length(files)
 	if (!is.null(cor.pheno)) {
 		if (nrow(cor.pheno) != m | ncol(cor.pheno) != m) {
@@ -72,7 +75,7 @@
 	}
 	cat('loading data ...\n')
 	data <- c()
-	fn <- rev(files)
+	fn <- files # rev(files)
 	for (i in m:1) {
 		dd <- read.table(fn[i], header = TRUE, stringsAsFactors = FALSE)
 		colnames(dd) <- tolower(colnames(dd))
@@ -131,9 +134,6 @@
 	gwa0 <- na.omit(gwa0)
 	cat('\n')
 	if (is.null(cor.pheno)) {
-		if (is.null(indep.snps)) {
-			stop('indep.snps required for estimating cor.pheno!')
-		}
 		n.ratio <- diag(m)
 		for (i in 1:(m - 1)) {
 			for (j in (i + 1):m) {

Modified: pkg/MultiABEL/src/MultiSummaryLoop.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoop.f90	2015-08-18 10:09:34 UTC (rev 2013)
+++ pkg/MultiABEL/src/MultiSummaryLoop.f90	2015-08-20 07:11:59 UTC (rev 2014)
@@ -4,7 +4,7 @@
 !	Created by Xia Shen on 5/29/15.
 !	
 
-subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, E0, H, invE0, s, coef, pil, fstat)
+subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, E0, H, invE0, D, s, coef, pil, fstat)
 
   implicit none
   integer :: i, j, k, m
@@ -31,7 +31,7 @@
      end do
      gg = vg*nn(j)
      bDbeta0 = matmul(matmul(b, D), transpose(beta))
-     bDbeta = bDbeta0(1,1)
+     bDbeta = bDbeta0(1,1)*vg
      do i = 1, m
         s(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)
      end do