[Genabel-commits] r2070 - in pkg/MultiABEL: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Jun 15 00:56:27 CEST 2017


Author: shenxia
Date: 2017-06-15 00:56:27 +0200 (Thu, 15 Jun 2017)
New Revision: 2070

Modified:
   pkg/MultiABEL/DESCRIPTION
   pkg/MultiABEL/NAMESPACE
   pkg/MultiABEL/R/load.summary.R
   pkg/MultiABEL/man/load.summary.Rd
Log:


Modified: pkg/MultiABEL/DESCRIPTION
===================================================================
--- pkg/MultiABEL/DESCRIPTION	2017-06-14 22:35:32 UTC (rev 2069)
+++ pkg/MultiABEL/DESCRIPTION	2017-06-14 22:56:27 UTC (rev 2070)
@@ -10,7 +10,8 @@
     statistics.
 Depends:
     R (>= 2.10),
-    svMisc
+    svMisc,
+    data.table
 Suggests:
     GenABEL,
     DatABEL

Modified: pkg/MultiABEL/NAMESPACE
===================================================================
--- pkg/MultiABEL/NAMESPACE	2017-06-14 22:35:32 UTC (rev 2069)
+++ pkg/MultiABEL/NAMESPACE	2017-06-14 22:56:27 UTC (rev 2070)
@@ -2,6 +2,6 @@
 importFrom("stats", "cor", "density", "lm", "na.omit", "pchisq", "pf",
              "qchisq", "qnorm", "quantile")
 importFrom("utils", "packageDescription", "read.table", "write.table")
-import("svMisc")
+import("svMisc", "data.table")
 export("Multivariate", "MultiRep", "MultiMeta", "MultiLoad","load.summary", "MultiSummary")
 exportClasses("multi.summary", "multi.loaded")
\ No newline at end of file

Modified: pkg/MultiABEL/R/load.summary.R
===================================================================
--- pkg/MultiABEL/R/load.summary.R	2017-06-14 22:35:32 UTC (rev 2069)
+++ pkg/MultiABEL/R/load.summary.R	2017-06-14 22:56:27 UTC (rev 2070)
@@ -5,7 +5,8 @@
 #' 
 #' @param files A vector of file names as strings. Each file name should contain summary statistics of
 #' one trait to be included in the multi-trait analysis. The columns of the summary statistics have to
-#' contain (uppercase or lowercase does not matter) \code{'snp'} (marker ID), \code{'a1'} (the first allele), \code{'freq'} 
+#' contain (uppercase or lowercase does not matter) \code{'snp'} (marker ID), \code{'a1'} (the first allele),
+#' \code{'a2'} (the second allele), \code{'freq'},
 #' (frequency of the first allele), \code{'beta'} (effect size), \code{'se'} (standard error), and 
 #' \code{'n'} (sample size).
 #' @param cor.pheno A #traits x #traits matrix of correlation matrix of the phenotypes, to be used to 
@@ -21,7 +22,7 @@
 #' If \code{TRUE}, each phenotypic variance will be estimated to adjust the summary statistics, so that
 #' the corresponding phenoypic variance is 1.
 #' @param columnNames A vector with names of columns containing necessary information in the input file;
-#' default values are c('snp','a1','freq','beta','se','n'). The values are case-insensitive. Note: check
+#' default values are c('snp','a1','a2','freq','beta','se','n'). The values are case-insensitive. Note: check
 #' your allele definitions for different traits are based on the same strand!
 #' @param fixedN sample size to assume across all analyses, when provided, this number will be used
 #' (instead of the ones specified in the input files)
@@ -30,14 +31,12 @@
 #' (the cleaned data to be processed in multi-trait GWAS), \code{cor.pheno} (user input or estimated), and 
 #' \code{var.pheno} (default or estimated).
 #' 
-#' @author Xia Shen, Yurii S. Aulchenko
+#' @author Xia Shen, Yakov, Tsepilov, Yurii S. Aulchenko
 #' 
 #' @references 
-#' Xia Shen, Zheng Ning, Yakov Tsepilov, Masoud Shirali, 
-#' Generation Scotland, Blair H. Smith, Lynne J. Hocking, Sandosh Padmanabhan, Caroline Hayward, 
-#' David J. Porteous, Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2015).
-#' Simple multi-trait analysis identifies novel loci 
-#' associated with growth and obesity measures. \emph{Submitted}.
+#' Xia Shen, Yakov Tsepilov, ..., Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2017).
+#' Discovery, replication, and in silico functional investigation of 
+#' 22 new pleiotropic anthropometric loci. \emph{Submitted}.
 #' 
 #' @seealso 
 #' \code{MultiSummary}
@@ -66,66 +65,91 @@
 #' @aliases load.summary
 #' @keywords multivariate, meta-analysis
 #' 
-`load.summary` <- function(files, cor.pheno = NULL, indep.snps = NULL, est.var = FALSE, 
-		                   columnNames = c('snp', 'a1', 'freq', 'beta', 'se', 'n'), fixedN = NULL) {
+`load.summary` <- function(files, cor.pheno = NULL, indep.snps = NULL, 
+		est.var = FALSE, 
+		columnNames = c('snp', 'a1', 'a2', 'freq', 'beta', 'se', 'n'),
+		fixedN = NULL) {
+	
+	### I. Sanity checks ###
+	
+	# require(data.table)
+	
 	if (!all(is.character(files))) {
 		stop('files should be given as strings!')
 	}
+	
 	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) {
 			stop('wrong dimensions of cor.pheno!')
 		}
 	}
+	
 	columnNames <- tolower(columnNames)
+	
 	if (!is.null(fixedN)) if (fixedN <= 0) {
 			stop('fixedN should be a positive number')
 		}
-	if (is.null(fixedN)) { colNamLen = 6 } else { colNamLen = 5 }
+	
+	if (is.null(fixedN)) { colNamLen = 7 } else { colNamLen = 6 }
+	
 	if (!is.character(columnNames)) {
 		stop('columnNames should be character')
 	}
+	
 	if (length(columnNames) != colNamLen) {
-		cat('columnNames should be a vector with', colNamLen, 'elements')
+		cat('columnNames should be a vector with',colNamLen,'elements')
 		stop('... exiting')
 	}
+	
 	if (length(unique(columnNames)) != colNamLen) {
 		stop('elements of columnNames must be unique')
 	}
-	# column Names Translation
+	
 	cNT = list(
-		'snp' = columnNames[1],
-		'a1'  = columnNames[2],
-		'freq'= columnNames[3],
-		'beta'= columnNames[4],
-		'se'  = columnNames[5],
-		'n'   = columnNames[6]
+			'snp' = columnNames[1],
+			'a1'  = columnNames[2],
+			'a2' = columnNames[3],
+			'freq'= columnNames[4],
+			'beta'= columnNames[5],
+			'se'  = columnNames[6],
+			'n'   = columnNames[7]
 	)
-	cat('loading data ...\n')
+	
+	### II. Data load ###
+	
+	cat("loading data ...\n")
+	
 	data <- c()
 	fn <- files # rev(files)
 	vys <- rep(1, m)
 	for (i in m:1) {
-		dd <- read.table(fn[i], header = TRUE, stringsAsFactors = FALSE)
+		
+		# Reserved for possible problems with data.table package
+		# dd <- read.table(fn[i], header = TRUE, stringsAsFactors = FALSE)
+		
+		dd <- data.table::fread(fn[i], header = TRUE, stringsAsFactors = FALSE,verbose = FALSE, showProgress = FALSE, data.table = FALSE)
+		
 		colnames(dd) <- tolower(colnames(dd))
 		currentColNames <- colnames(dd)
 		if (any(!(columnNames %in% currentColNames))) {
 			stop('file column names do not match columnNames in ', fn[i], '... exiting!')
 		}
-		dd[,cNT$a1] <- toupper(dd[,cNT$a1])
-		idx <- which(duplicated(dd[,cNT$snp]))
+		idx <- which(duplicated(dd[, cNT$snp]))
 		if (length(idx) > 0) {
 			data[[i]] <- dd[-idx,]
-			rownames(data[[i]]) <- dd[-idx,cNT$snp]
+			rownames(data[[i]]) <- dd[ -idx , cNT$snp]
 		} else {
 			data[[i]] <- dd
-			rownames(data[[i]]) <- dd[,cNT$snp]
+			rownames(data[[i]]) <- dd[, cNT$snp]
 		}
 		if (est.var) {
 			if (is.null(fixedN)) {
@@ -142,43 +166,79 @@
 	}
 	cat('\n')
 	if (est.var) cat('phenotypic variances are:', vys, '\n')
-	cat('checking markers ...\n')
-	snps <- data[[1]][, cNT$snp]
-	for (i in 2:m) {
-		snps <- data[[i]][ data[[i]][, cNT$snp] %in% snps, cNT$snp]
-		progress(i/m*100)
+	
+	### III. SNP cheking and cleaning ###
+	
+	cat("checking markers ...\n")
+	
+	snps <- data[[1]][,cNT$snp]
+	for (i in 1:m) {
+		data[[i]]=na.omit(data[[i]][,unlist(cNT)])
+		snps=intersect(snps,data[[i]][, cNT$snp])
+		progress(i/m * 100)
 	}
+	
 	snps <- unique(snps)
-	cat('\n')
-	cat('cleaning data ...\n')
+	cat("\n")
+	
+	cat("cleaning data ...\n")
 	for (i in 1:m) {
-		data[[i]] <- data[[i]][snps,]
-		progress(i/m*100)
+		data[[i]] <- data[[i]][snps, ]
+		progress(i/m * 100)
 	}
-	cat('\n')
-	cat('correcting parameters ...\n')
+	cat("\n")
+	cat("correcting parameters ...\n")
+	
+	for (i in 1:m) {
+		data[[i]][, cNT[["a1"]]] = toupper(data[[i]][, cNT[["a1"]]])
+		data[[i]][, cNT[["a2"]]] = toupper(data[[i]][, cNT[["a2"]]])
+	}
+	
+	d1a12 <- paste(data[[1]][,c(cNT$a1)],data[[1]][,c(cNT$a2)], sep = "")
+	d1a21 <- paste(data[[1]][,c(cNT$a2)],data[[1]][,c(cNT$a1)], sep = "")
+	
+	
+	snps=rownames(data[[1]])
+	
+	i <- 2
 	for (i in 2:m) {
-		if (any(data[[i]][, cNT[['a1']]] != data[[1]][, cNT[['a1']]])) {
-			adj <- 2*as.numeric(data[[i]][, cNT[['a1']]] == data[[1]][, cNT[['a1']]]) - 1
-			data[[i]][, cNT$beta] <- data[[i]][, cNT$beta]*adj
-			data[[i]][, cNT$freq] <- (adj == 1)*data[[i]][, cNT$freq] + (adj == -1)*(1 - data[[i]][,cNT$freq])
+		
+		dia12 <- paste(data[[i]][,c(cNT$a1)],data[[i]][,c(cNT$a2)], sep = "")
+		
+		ind=which(((dia12==d1a12) + (dia12==d1a21))==1)
+		snps=intersect(rownames(data[[i]])[ind],snps)
+		
+		if (any(data[[i]][, cNT[["a1"]]] != data[[1]][, cNT[["a1"]]])) {
+			adj <- 2 * as.numeric(data[[i]][, cNT[["a1"]]] == data[[1]][, cNT[["a1"]]]) - 1
+			data[[i]][, cNT$beta] <- data[[i]][, cNT$beta] *adj
+			data[[i]][, cNT$freq] <- (adj == 1) * data[[i]][,cNT$freq] + (adj == -1) * (1 - data[[i]][, cNT$freq])
 		}
-		progress(i/m*100)
+		progress(i/m * 100)
 	}
-	cat('\n')
-	cat('adjusting sample size ... ')
+	
+	for (i in 1:m) {
+		data[[i]] <- data[[i]][snps, ]
+		#progress(i/m * 100)
+	}
+	
+	cat("\n")
+	
+	### IV. Sample size and finalisation ###
+	
+	cat("adjusting sample size ... ")
 	n0 <- matrix(NA, nrow(data[[1]]), m)
 	if (is.null(fixedN)) {
 		for (i in 1:m) {
-			n0[,i] <- data[[i]][,cNT$n]
+			n0[, i] <- data[[i]][, cNT$n]
 		}
-	} else {
+	}else {
 		for (i in 1:m) {
-			n0[,i] <- fixedN
-		}            
+			n0[, i] <- fixedN
+		}
 	}
 	n <- apply(n0, 1, "min")
-	cat('done.\n')
+	cat("done.\n")
+	
 	cat('finalizing summary statistics ...\n')
 	gwa0 <- matrix(NA, nrow(data[[1]]), 2*m + 2)
 	for (i in 1:m) {
@@ -191,6 +251,7 @@
 	rownames(gwa0) <- data[[1]][,cNT$snp]
 	gwa0 <- na.omit(gwa0)
 	cat('\n')
+	
 	if (is.null(cor.pheno)) {
 		n.ratio <- diag(m)
 		for (i in 1:(m - 1)) {
@@ -211,10 +272,16 @@
 		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, var.pheno = vys)
-	class(dd) <- 'multi.summary'
+	
+	A12 <- data[[1]][rownames(gwa0), c(cNT[["a1"]],cNT[["a2"]])]
+	colnames(A12) <- c("A1","A2")
+	
+	dd <- list(gwa = gwa0, cor.pheno = cor.pheno, var.pheno = vys, alleles = A12)
+	class(dd) <- "multi.summary"
 	return(dd)
+	
 }

Modified: pkg/MultiABEL/man/load.summary.Rd
===================================================================
--- pkg/MultiABEL/man/load.summary.Rd	2017-06-14 22:35:32 UTC (rev 2069)
+++ pkg/MultiABEL/man/load.summary.Rd	2017-06-14 22:56:27 UTC (rev 2070)
@@ -5,12 +5,14 @@
 \title{Loading multiple summary statistics from genome-wide association studies}
 \usage{
 load.summary(files, cor.pheno = NULL, indep.snps = NULL, est.var = FALSE,
-  columnNames = c("snp", "a1", "freq", "beta", "se", "n"), fixedN = NULL)
+  columnNames = c("snp", "a1", "a2", "freq", "beta", "se", "n"),
+  fixedN = NULL)
 }
 \arguments{
 \item{files}{A vector of file names as strings. Each file name should contain summary statistics of
 one trait to be included in the multi-trait analysis. The columns of the summary statistics have to
-contain (uppercase or lowercase does not matter) \code{'snp'} (marker ID), \code{'a1'} (the first allele), \code{'freq'} 
+contain (uppercase or lowercase does not matter) \code{'snp'} (marker ID), \code{'a1'} (the first allele),
+\code{'a2'} (the second allele), \code{'freq'},
 (frequency of the first allele), \code{'beta'} (effect size), \code{'se'} (standard error), and 
 \code{'n'} (sample size).}
 
@@ -30,7 +32,7 @@
 the corresponding phenoypic variance is 1.}
 
 \item{columnNames}{A vector with names of columns containing necessary information in the input file;
-default values are c('snp','a1','freq','beta','se','n'). The values are case-insensitive. Note: check
+default values are c('snp','a1','a2','freq','beta','se','n'). The values are case-insensitive. Note: check
 your allele definitions for different traits are based on the same strand!}
 
 \item{fixedN}{sample size to assume across all analyses, when provided, this number will be used
@@ -68,17 +70,15 @@
 }
 }
 \references{
-Xia Shen, Zheng Ning, Yakov Tsepilov, Masoud Shirali, 
-Generation Scotland, Blair H. Smith, Lynne J. Hocking, Sandosh Padmanabhan, Caroline Hayward, 
-David J. Porteous, Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2015).
-Simple multi-trait analysis identifies novel loci 
-associated with growth and obesity measures. \emph{Submitted}.
+Xia Shen, Yakov Tsepilov, ..., Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2017).
+Discovery, replication, and in silico functional investigation of 
+22 new pleiotropic anthropometric loci. \emph{Submitted}.
 }
 \seealso{
 \code{MultiSummary}
 }
 \author{
-Xia Shen, Yurii S. Aulchenko
+Xia Shen, Yakov, Tsepilov, Yurii S. Aulchenko
 }
 \keyword{meta-analysis}
 \keyword{multivariate,}



More information about the Genabel-commits mailing list