[Genabel-commits] r690 - in pkg/VariABEL: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 8 12:07:55 CET 2011
Author: maksim
Date: 2011-03-08 12:07:55 +0100 (Tue, 08 Mar 2011)
New Revision: 690
Added:
pkg/VariABEL/R/var_test_gwaa.R
Removed:
pkg/VariABEL/R/var.test.gwaa.R
pkg/VariABEL/data/
pkg/VariABEL/demo/
Modified:
pkg/VariABEL/NAMESPACE
pkg/VariABEL/cleanup
pkg/VariABEL/man/var_test_gwaa.Rd
Log:
fix errors from R CMD check
Modified: pkg/VariABEL/NAMESPACE
===================================================================
--- pkg/VariABEL/NAMESPACE 2011-03-07 23:41:59 UTC (rev 689)
+++ pkg/VariABEL/NAMESPACE 2011-03-08 11:07:55 UTC (rev 690)
@@ -3,7 +3,7 @@
export(
# var.meta.gwaa,
- var.test.gwaa#,
+ var_test_gwaa#,
#test_databel
)
Deleted: pkg/VariABEL/R/var.test.gwaa.R
===================================================================
--- pkg/VariABEL/R/var.test.gwaa.R 2011-03-07 23:41:59 UTC (rev 689)
+++ pkg/VariABEL/R/var.test.gwaa.R 2011-03-08 11:07:55 UTC (rev 690)
@@ -1,230 +0,0 @@
-#=====================================================================================
-#
-# Filename: var.test.gwaa.R
-#
-# Description: VarABEL functions: test each SNP in gwaa data for homogeneity of trait's variance among genotypes.
-#
-# Version: 0.1
-# Created: ---
-# Revision: none
-# last modification: 26-Apr-2010
-#
-# Author: Maksim V. Struchalin
-# Modified by: Maksim V. Struchalin, 11-Jan-2009
-# Company: ErasmusMC, Epidemiology Department, The Netherlands.
-# Email: m.struchalin at erasmusmc.nl
-# license: GPL (>=2)
-#
-#=====================================================================================
-
-
-
-
-
-#_____________________________________________________________________________________________________________________________
-#Available for user function for genome-wide testing variance homogeneity of trait's distribution.
-#
-"var_test_gwaa" <- function(formula, genodata, phenodata, genodata_info=NULL, testname="sqlm", analysis_type="AAvsABvsBB")
-{
-
-
-if(testname == "bartlett") {testname <- 0}
-else if(testname == "levene") {testname <- 1}
-else if(testname == "likelihood") {testname <- 2}
-else if(testname == "kolmogorov_smirnov") {testname <- 3}
-else if(testname == "sqlm") {testname <- 4}
-else {stop(paste(testname, "is unsupportive type of test. Only levene, and sqlm are supported."))}
-
-
-if(analysis_type == "AAvsABvsBB") {analysis_type <- 0}
-else if(analysis_type == "AAvsABandBB") {analysis_type <- 1}
-else if(analysis_type == "ABvsAAandBB") {analysis_type <- 2}
-else if(analysis_type == "BBvsAAandAB") {analysis_type <- 3}
-else {stop(paste(analysis_type, "is unsupportive type of analysis. Only AAvsABvsBB, AAvsABandBB, ABvsAAandBB, and BBvsAAandAB are supported."))}
-
-MAR <- 2
-OUT <- "R"
-FUN <- "variance_homogeneity_test_C_wrapper"
-
-
-idnum <- 0
-
-if(is(formula, "formula"))
- {
- design_matrix <- model.matrix(formula, data = phenodata)
-# design_matrix <- design_matrix_
- trat_name <- as.character(formula[[2]])
- trait <- phenodata[,trat_name]
-# idnum <- length(trait)
- design_matrix_df <- data.frame(design_matrix)
-
-#Put NAs back to the matrix
- if(dim(design_matrix_df)[1] != dim(phenodata)[1])
- {
- design_matrix_df_with_nas <- matrix(rep(NA, dim(phenodata)[1]*dim(design_matrix_df)[2]), ncol=dim(design_matrix_df)[2], nrow=dim(phenodata)[1])
- design_matrix_df_with_nas <- data.frame(design_matrix_df_with_nas)
- formula_vars <- all.vars(formula)
- na_bool <- rep(T, dim(phenodata)[1])
- for(term in formula_vars)
- {
- na_bool <- na_bool & !is.na(phenodata[,term])
- }
- colnames(design_matrix_df_with_nas) <- colnames(design_matrix_df)
- design_matrix_df_with_nas[na_bool,] <- design_matrix_df
- design_matrix_df <- design_matrix_df_with_nas
- }
-
- design_matrix_df$snp <- 0 #will be filled into iterator
- idnum <- dim(design_matrix_df)[1]
-
-
- }
-# else if(is(formula, "numeric") || is(formula, "integer") || is(formula, "double"))
- else if(is.character(formula))
- {
- trait <- phenodata[,formula]
- idnum <- length(trait)
- design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
- }
-
-#design_matrix_geno_means_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
-
-
-p <- dim(design_matrix_df)[2] #covariates number + 1
-
-gtNrow <- NA
-gtNcol <- NA
-
-genodata_info_df <- data.frame()
-
-
-if(class(genodata) == "snp.data")
- {
- print("reading genotype info...")
- genodata_info_df <- data.frame(name=genodata at snpnames,
- coding=as.character(genodata at coding),
- strand=as.character(genodata at strand),
- chromosome=as.character(genodata at chromosome),
- map=genodata at map
- )
-# gtNrow <- genodata at nsnps
-# gtNcol <- genodata at nids
- gtNrow <- genodata at nids
- gtNcol <- genodata at nsnps
- genodata <- as.raw(genodata at gtps)
- }
- else if (class(genodata)=="databel")
- {
- gtNrow <- dim(genodata)[1]
- gtNcol <- dim(genodata)[2]
- genodata <- genodata at data
-
- if(!is.null(genodata_info))
- {
- if(file.exists(genodata_info))
- {
- print("reading genotype's info...")
- genodata_info_df <- read.table(genodata_info, strings=F, header=T)
-# print(c(dim(genodata_info_df)[1], gtNcol))
- if(dim(genodata_info_df)[1] != gtNcol)
- {
- print(paste("warning: File ", genodata_info, " contains information about ", dim(genodata_info_df)[1],
- " SNPs but genotype data contains ", gtNcol, ". Ignoring of the info file.", sep=""))
- genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
- }
- }
- else
- {
- print(paste("warning: File ", genodata_info, "does not exist. Trying to run analysis withouth SNP's information..."))
- genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
- }
- }
- else
- {
- genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
- }
-
- }
- else
- {
- stop(paste("genodata class not recognized ('",class(genodata),"')",sep=""))
- }
-
-
-if(testname == 1) #Levene's
- {
- if(is(formula, "formula"))
- {
- trait_regression <- lm(formula, data=phenodata)
- trait <- trait_regression$residuals
- }
- }
-
-
-#print(design_matrix_df)
-#print(dim(design_matrix_df))
-
-
-print("Start variance analysis...")
-results_C <- .Call("iterator", genodata,
- as.integer(gtNrow), as.integer(gtNcol),
- as.character(FUN),
- as.character(OUT),
- as.integer(MAR),
- as.integer(1),
- as.integer(18), #iterator additional inputa parameters number
- as.double(trait),
- as.double(data.matrix(design_matrix_df)),
- as.double(data.matrix(design_matrix_df)),
- as.integer(p),
- as.integer(analysis_type),
- as.integer(testname),
- double(p+1),#betas
- double(p+1),#se
- double(1),#chi2
- integer(1),#df
- double(idnum),#residuals
- double(idnum),#qty
- integer(p),#jpvt
- double(p),#qraux
- double(2*p),#work
- double(p*p),#v
- double(p*p),#x_for_ch2inv
- integer(idnum)
- )
-
-print("Variance analysis done.")
-
-#print(genodata_info_df)
-#print(results_C)
-
-results_df <- data.frame(genodata_info_df, results_C)
-#print(results_df)
-
-cov_names <- colnames(design_matrix_df)
-
-#print(cov_names)
-
-output_column_names <- c(colnames(genodata_info_df), "chisq", "df", "Intercept_effect", "Intercept_sd")
-
-
-#print(output_column_names)
-
-for(i in 2:p)
- {
- output_column_names <- c(output_column_names, paste(cov_names[i], "_eff", sep=""))
- output_column_names <- c(output_column_names, paste(cov_names[i], "_se", sep=""))
- }
-
-output_column_names <- c(output_column_names, "snp_eff_dispertion")
-output_column_names <- c(output_column_names, "snp_se_dispertion")
-
-colnames(results_df) <- output_column_names
-
-
-return(results_df)
-}
-
-
-
-
Added: pkg/VariABEL/R/var_test_gwaa.R
===================================================================
--- pkg/VariABEL/R/var_test_gwaa.R (rev 0)
+++ pkg/VariABEL/R/var_test_gwaa.R 2011-03-08 11:07:55 UTC (rev 690)
@@ -0,0 +1,230 @@
+#=====================================================================================
+#
+# Filename: var.test.gwaa.R
+#
+# Description: VarABEL functions: test each SNP in gwaa data for homogeneity of trait's variance among genotypes.
+#
+# Version: 0.1
+# Created: ---
+# Revision: none
+# last modification: 26-Apr-2010
+#
+# Author: Maksim V. Struchalin
+# Modified by: Maksim V. Struchalin, 11-Jan-2009
+# Company: ErasmusMC, Epidemiology Department, The Netherlands.
+# Email: m.struchalin at erasmusmc.nl
+# license: GPL (>=2)
+#
+#=====================================================================================
+
+
+
+
+
+#_____________________________________________________________________________________________________________________________
+#Available for user function for genome-wide testing variance homogeneity of trait's distribution.
+#
+"var_test_gwaa" <- function(formula, genodata, phenodata, genodata_info=NULL, testname="sqlm", analysis_type="AAvsABvsBB")
+{
+
+
+if(testname == "bartlett") {testname <- 0}
+else if(testname == "levene") {testname <- 1}
+else if(testname == "likelihood") {testname <- 2}
+else if(testname == "kolmogorov_smirnov") {testname <- 3}
+else if(testname == "sqlm") {testname <- 4}
+else {stop(paste(testname, "is unsupportive type of test. Only levene, and sqlm are supported."))}
+
+
+if(analysis_type == "AAvsABvsBB") {analysis_type <- 0}
+else if(analysis_type == "AAvsABandBB") {analysis_type <- 1}
+else if(analysis_type == "ABvsAAandBB") {analysis_type <- 2}
+else if(analysis_type == "BBvsAAandAB") {analysis_type <- 3}
+else {stop(paste(analysis_type, "is unsupportive type of analysis. Only AAvsABvsBB, AAvsABandBB, ABvsAAandBB, and BBvsAAandAB are supported."))}
+
+MAR <- 2
+OUT <- "R"
+FUN <- "variance_homogeneity_test_C_wrapper"
+
+
+idnum <- 0
+
+if(is(formula, "formula"))
+ {
+ design_matrix <- model.matrix(formula, data = phenodata)
+# design_matrix <- design_matrix_
+ trat_name <- as.character(formula[[2]])
+ trait <- phenodata[,trat_name]
+# idnum <- length(trait)
+ design_matrix_df <- data.frame(design_matrix)
+
+#Put NAs back to the matrix
+ if(dim(design_matrix_df)[1] != dim(phenodata)[1])
+ {
+ design_matrix_df_with_nas <- matrix(rep(NA, dim(phenodata)[1]*dim(design_matrix_df)[2]), ncol=dim(design_matrix_df)[2], nrow=dim(phenodata)[1])
+ design_matrix_df_with_nas <- data.frame(design_matrix_df_with_nas)
+ formula_vars <- all.vars(formula)
+ na_bool <- rep(T, dim(phenodata)[1])
+ for(term in formula_vars)
+ {
+ na_bool <- na_bool & !is.na(phenodata[,term])
+ }
+ colnames(design_matrix_df_with_nas) <- colnames(design_matrix_df)
+ design_matrix_df_with_nas[na_bool,] <- design_matrix_df
+ design_matrix_df <- design_matrix_df_with_nas
+ }
+
+ design_matrix_df$snp <- 0 #will be filled into iterator
+ idnum <- dim(design_matrix_df)[1]
+
+
+ }
+# else if(is(formula, "numeric") || is(formula, "integer") || is(formula, "double"))
+ else if(is.character(formula))
+ {
+ trait <- phenodata[,formula]
+ idnum <- length(trait)
+ design_matrix_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
+ }
+
+#design_matrix_geno_means_df <- data.frame(X.Intercept=rep(1,idnum), snp=rep(0,idnum))
+
+
+p <- dim(design_matrix_df)[2] #covariates number + 1
+
+gtNrow <- NA
+gtNcol <- NA
+
+genodata_info_df <- data.frame()
+
+
+if(class(genodata) == "snp.data")
+ {
+ print("reading genotype info...")
+ genodata_info_df <- data.frame(name=genodata at snpnames,
+ coding=as.character(genodata at coding),
+ strand=as.character(genodata at strand),
+ chromosome=as.character(genodata at chromosome),
+ map=genodata at map
+ )
+# gtNrow <- genodata at nsnps
+# gtNcol <- genodata at nids
+ gtNrow <- genodata at nids
+ gtNcol <- genodata at nsnps
+ genodata <- as.raw(genodata at gtps)
+ }
+ else if (class(genodata)=="databel")
+ {
+ gtNrow <- dim(genodata)[1]
+ gtNcol <- dim(genodata)[2]
+ genodata <- genodata at data
+
+ if(!is.null(genodata_info))
+ {
+ if(file.exists(genodata_info))
+ {
+ print("reading genotype's info...")
+ genodata_info_df <- read.table(genodata_info, strings=F, header=T)
+# print(c(dim(genodata_info_df)[1], gtNcol))
+ if(dim(genodata_info_df)[1] != gtNcol)
+ {
+ print(paste("warning: File ", genodata_info, " contains information about ", dim(genodata_info_df)[1],
+ " SNPs but genotype data contains ", gtNcol, ". Ignoring of the info file.", sep=""))
+ genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
+ }
+ }
+ else
+ {
+ print(paste("warning: File ", genodata_info, "does not exist. Trying to run analysis withouth SNP's information..."))
+ genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
+ }
+ }
+ else
+ {
+ genodata_info_df <- data.frame(name=sub("^", "snp", as.character(1:gtNcol)))
+ }
+
+ }
+ else
+ {
+ stop(paste("genodata class not recognized ('",class(genodata),"')",sep=""))
+ }
+
+
+if(testname == 1) #Levene's
+ {
+ if(is(formula, "formula"))
+ {
+ trait_regression <- lm(formula, data=phenodata)
+ trait <- trait_regression$residuals
+ }
+ }
+
+
+#print(design_matrix_df)
+#print(dim(design_matrix_df))
+
+
+print("Start variance analysis...")
+results_C <- .Call("iterator", genodata,
+ as.integer(gtNrow), as.integer(gtNcol),
+ as.character(FUN),
+ as.character(OUT),
+ as.integer(MAR),
+ as.integer(1),
+ as.integer(18), #iterator additional inputa parameters number
+ as.double(trait),
+ as.double(data.matrix(design_matrix_df)),
+ as.double(data.matrix(design_matrix_df)),
+ as.integer(p),
+ as.integer(analysis_type),
+ as.integer(testname),
+ double(p+1),#betas
+ double(p+1),#se
+ double(1),#chi2
+ integer(1),#df
+ double(idnum),#residuals
+ double(idnum),#qty
+ integer(p),#jpvt
+ double(p),#qraux
+ double(2*p),#work
+ double(p*p),#v
+ double(p*p),#x_for_ch2inv
+ integer(idnum)
+ )
+
+print("Variance analysis done.")
+
+#print(genodata_info_df)
+#print(results_C)
+
+results_df <- data.frame(genodata_info_df, results_C)
+#print(results_df)
+
+cov_names <- colnames(design_matrix_df)
+
+#print(cov_names)
+
+output_column_names <- c(colnames(genodata_info_df), "chisq", "df", "Intercept_effect", "Intercept_sd")
+
+
+#print(output_column_names)
+
+for(i in 2:p)
+ {
+ output_column_names <- c(output_column_names, paste(cov_names[i], "_eff", sep=""))
+ output_column_names <- c(output_column_names, paste(cov_names[i], "_se", sep=""))
+ }
+
+output_column_names <- c(output_column_names, "snp_eff_dispertion")
+output_column_names <- c(output_column_names, "snp_se_dispertion")
+
+colnames(results_df) <- output_column_names
+
+
+return(results_df)
+}
+
+
+
+
Modified: pkg/VariABEL/cleanup
===================================================================
--- pkg/VariABEL/cleanup 2011-03-07 23:41:59 UTC (rev 689)
+++ pkg/VariABEL/cleanup 2011-03-08 11:07:55 UTC (rev 690)
@@ -1,3 +1,5 @@
rm -f configure.Rout
rm -rf .RData R/.RData
rm -f src/*.o
+rm -f src/*.cpp
+rm -f src/*.h
Modified: pkg/VariABEL/man/var_test_gwaa.Rd
===================================================================
--- pkg/VariABEL/man/var_test_gwaa.Rd 2011-03-07 23:41:59 UTC (rev 689)
+++ pkg/VariABEL/man/var_test_gwaa.Rd 2011-03-08 11:07:55 UTC (rev 690)
@@ -39,7 +39,7 @@
% \code{\link{var.meta}},
%}
\examples{
- library(GenABEL)
+ require(GenABEL)
data(srdta)
result1 <- var_test_gwaa(bt~qt1+qt2, genodata=srdta at gtdata, phenodata=srdta at phdata)
More information about the Genabel-commits
mailing list