[Genabel-commits] r1189 - in pkg/GenABEL: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 9 17:37:06 CEST 2013
Author: yurii
Date: 2013-04-09 17:37:06 +0200 (Tue, 09 Apr 2013)
New Revision: 1189
Modified:
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/R/cocohet.R
pkg/GenABEL/man/cocohet.Rd
Log:
Fixing bug [#2672] (GenABEL::cocohet default graph output)
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2013-04-09 15:17:56 UTC (rev 1188)
+++ pkg/GenABEL/CHANGES.LOG 2013-04-09 15:37:06 UTC (rev 1189)
@@ -1,5 +1,12 @@
*** v. 1.7-5
+(2013.04.09)
+Fixing bug [#2672] (GenABEL::cocohet default graph output). Now by
+default the function does not produce graphics (makePlot = FALSE); if
+graphics is requested (makePlot = TRUE) it is produced using a default
+plot device (so a user can redirect the output to a file if (s)he
+desires and in the way (s)he desires).
+
(2013.04.04)
Fixing bug [#2525] (http://r-forge.r-project.org/tracker/index.php?func=detail&aid=2525&group_id=505&atid=2058)
Thanks to Vladimir Naumov for submitting the patch!
Modified: pkg/GenABEL/R/cocohet.R
===================================================================
--- pkg/GenABEL/R/cocohet.R 2013-04-09 15:17:56 UTC (rev 1188)
+++ pkg/GenABEL/R/cocohet.R 2013-04-09 15:37:06 UTC (rev 1189)
@@ -1,3 +1,5 @@
+#' Test for compound heterozygote effects
+#'
#' Detecting rare recessive and compound heterozygote alleles in genome wide
#' association.
#'
@@ -29,8 +31,8 @@
#' element is a chisq obtained as a maximum chisq between a given SNP and SNPs
#' on the right within a window. If TRUE then return also a matrix where
#' chisq's for all tests are stored. Object of class \code{\link{logical}}
-#' @param manhettan_filename File name where manhettan plot will be saved after
-#' analysis. Object of class \code{\link{character}}
+#' @param makePlot whether the Manhattan-type plot should be produced
+#' (TRUE or FALSE)
#' @param test Name of the test to be performed. Available tests are "CHI2",
#' "YATES" (chi2 with Yates correction), and "FISHER". Object of class
#' \code{\link{character}}
@@ -57,11 +59,11 @@
#' @keywords manip
#' @examples
#' data(srdta)
-#' chis2_nocorrection <- cocohet(data=srdta at gtdata, trait=srdta at phdata$bt, window=3, test="CHI2")
+#' chis2_nocorrection <- cocohet(data=gtdata(srdta),
+#' trait=phdata(srdta)$bt, window=3, test="CHI2")
#'
-"cocohet" <- function(data, trait, window, return_all_result=TRUE,
- manhettan_filename="manhettan_plot.jpeg", test="CHI2",
- min_expected_cut_off=-1)
+"cocohet" <- function(data, trait, window, return_all_result = TRUE,
+ makePlot = FALSE, test = "CHI2", min_expected_cut_off = -1)
{
#Correction is performed in case expected value one of the margin in contigency table less 5. Available values for chi2_correction:
# "CHI2"
@@ -78,21 +80,21 @@
if(dim(table(trait)) != 2) {
stop("\"trait\" must be binary")
}
-
+
if(!is(window, "integer") && !is(window, "numeric")) {
stop("Wrong data class: the argument \"window\" should has type \"integer\" or \"numeric\".")
}
snp_number <- data at nsnps
-
+
if(window >= snp_number) {
print(paste("warning: \"window\" value ", window, " is bigger than total snp number. Automaticly reduced to maximum possible."))
window <- snp_number - 1
print(paste("New window value is ", window, ".\n"))
}
-
+
if(window==0) {stop("Input parameter window can not be zero.\n")}
-
+
output <- .Call("interaction_rare_recesive_allele_C_", as.raw(data at gtps), as.integer(data at nids), as.integer(snp_number),
as.integer(trait),
as.integer(window),
@@ -104,21 +106,17 @@
if(return_all_result) {
results$chi2_all <- matrix(output[data at nsnps:length(output)], ncol=window, nrow=data at nsnps-1, byrow=T)
}
-
- print("plotting manhettan plot picture...")
- bitmap(manhettan_filename, type="jpeg")
- results_pval <- pchisq(results$chi2_max, df=1, lower.tail=F)
- results_list <- list()
- results_list$P1df <- results_pval
- results_list$map <- data at map[1:(data at nsnps-1)]
- results_list$chromosome <- data at chromosome[1:(data at nsnps-1)]
- class(results_list) <- "scan.gwaa"
- print("1")
- try(plot.scan.gwaa(results_list))
- print("2")
- dev.off()
- print(paste("The plot is saved into file ", manhettan_filename, sep=""))
+ if (makePlot) {
+ results_pval <- pchisq(results$chi2_max, df=1, lower.tail=F)
+ results_list <- list()
+ results_list$P1df <- results_pval
+ results_list$map <- data at map[1:(data at nsnps-1)]
+ results_list$chromosome <- data at chromosome[1:(data at nsnps-1)]
+ class(results_list) <- "scan.gwaa"
+ plot.scan.gwaa(results_list)
+ }
+
snp_position_top <- which.max(results$chi2_max)
snp_name_top <- data at snpnames[snp_position_top]
snp_chromosome_top <- data at chromosome[snp_position_top]
@@ -135,77 +133,78 @@
results
}
-#TODO it seems the function below is never used - should it be deleted from the code?
-"look_at_snp_chi2s_interactions_rare_recesive_alleles" <- function(
- chi2_all, data, snp_position,
- manhettan_filename="one_snp_manhettan.jpeg")
-{
+#TODO it seems the function below is never used - should it be deleted from the code?
+# Commenting it out for the time-being
+#"look_at_snp_chi2s_interactions_rare_recesive_alleles" <- function(
+# chi2_all, data, snp_position,
+# manhettan_filename="one_snp_manhettan.jpeg")
+#{
+#
+# if(!is(data,"snp.data")) {
+# stop("Wrong data class: the argument \"data\" should has type \"snp.data\".")
+# }
+#
+# if(!is(chi2_all,"matrix")) {
+# stop("Wrong data class: the argument \"chi2_all\" should has type \"matrix\".")
+# }
+#
+# window <- dim(chi2_all)[2]
+# snp_number <- dim(chi2_all)[1] + 1
+#
+# if(window == 0 | snp_number==0) stop("Wrong chi2 matrix.")
+#
+# results <- list()
+# results$right_window_chi2s <- chi2_all[snp_position,]
+# results$left_window_chi2s <- c()
+#
+# left_border <- snp_position - window
+# if(left_border <=0 ) left_border <- 1
+#
+# right_border <- snp_position + window
+# if(right_border > snp_number-1) right_border <- snp_number-1
+#
+# for(i in 1:left_border)
+# {
+# results$left_window_chi2s <- c(results$left_window_chi2s, chi2_all[snp_position-i, i])
+# }
+#
+# bitmap(manhettan_filename, type="jpeg")
+#
+# print(snp_position)
+# results_chi2 <- c(results$left_window_chi2s, results$right_window_chi2s)
+# print(length(results$left_window_chi2s))
+# print(length(results$right_window_chi2s))
+# print(length(results_chi2))
+# print("________________")
+# results_pval <- pchisq(results_chi2, df=1, lower.tail=F)
+#
+# max_position_local <- which.max(results_chi2)
+# max_position <- max_position_local + snp_position
+#
+# results_list <- list()
+# results_list$P1df <- results_pval
+# results_list$map <- c(data at map[left_border:(snp_position-1)], data at map[(snp_position+1):right_border])
+# results_list$chromosome <- as.factor(c(as.character(data at chromosome[left_border:(snp_position-1)]),
+# as.character(data at chromosome[(snp_position+1):right_border])))
+# class(results_list) <- "scan.gwaa"
+#
+#
+# print(length(results_list$P1df))
+# print(length(results_list$map))
+# print(length(results_list$chromosome))
+#
+# try(plot.scan.gwaa(results_list))
+#
+# dev.off()
+#
+# print(paste("The plot is saved into file ", manhettan_filename, sep=""))
+#
+# print(paste("Maximum chi2 is ", results_chi2[max_position_local], " for a snp in position ", max_position, sep=""))
+#
+#
+# results_chi2
+#}
- if(!is(data,"snp.data")) {
- stop("Wrong data class: the argument \"data\" should has type \"snp.data\".")
- }
-
- if(!is(chi2_all,"matrix")) {
- stop("Wrong data class: the argument \"chi2_all\" should has type \"matrix\".")
- }
-
- window <- dim(chi2_all)[2]
- snp_number <- dim(chi2_all)[1] + 1
-
- if(window == 0 | snp_number==0) stop("Wrong chi2 matrix.")
-
- results <- list()
- results$right_window_chi2s <- chi2_all[snp_position,]
- results$left_window_chi2s <- c()
-
- left_border <- snp_position - window
- if(left_border <=0 ) left_border <- 1
-
- right_border <- snp_position + window
- if(right_border > snp_number-1) right_border <- snp_number-1
-
- for(i in 1:left_border)
- {
- results$left_window_chi2s <- c(results$left_window_chi2s, chi2_all[snp_position-i, i])
- }
-
- bitmap(manhettan_filename, type="jpeg")
-
- print(snp_position)
- results_chi2 <- c(results$left_window_chi2s, results$right_window_chi2s)
- print(length(results$left_window_chi2s))
- print(length(results$right_window_chi2s))
- print(length(results_chi2))
- print("________________")
- results_pval <- pchisq(results_chi2, df=1, lower.tail=F)
-
- max_position_local <- which.max(results_chi2)
- max_position <- max_position_local + snp_position
-
- results_list <- list()
- results_list$P1df <- results_pval
- results_list$map <- c(data at map[left_border:(snp_position-1)], data at map[(snp_position+1):right_border])
- results_list$chromosome <- as.factor(c(as.character(data at chromosome[left_border:(snp_position-1)]),
- as.character(data at chromosome[(snp_position+1):right_border])))
- class(results_list) <- "scan.gwaa"
-
-
- print(length(results_list$P1df))
- print(length(results_list$map))
- print(length(results_list$chromosome))
-
- try(plot.scan.gwaa(results_list))
-
- dev.off()
-
- print(paste("The plot is saved into file ", manhettan_filename, sep=""))
-
- print(paste("Maximum chi2 is ", results_chi2[max_position_local], " for a snp in position ", max_position, sep=""))
-
-
- results_chi2
-}
-
# original header moved here because the roxygen doc is not at the beginning
#=====================================================================================
#
Modified: pkg/GenABEL/man/cocohet.Rd
===================================================================
--- pkg/GenABEL/man/cocohet.Rd 2013-04-09 15:17:56 UTC (rev 1188)
+++ pkg/GenABEL/man/cocohet.Rd 2013-04-09 15:37:06 UTC (rev 1189)
@@ -1,12 +1,11 @@
\name{cocohet}
\alias{chi2_CG}
\alias{cocohet}
-\title{Detecting rare recessive and compound heterozygote alleles in genome wide
-association.}
+\title{Test for compound heterozygote effects}
\usage{
cocohet(data, trait, window, return_all_result = TRUE,
- manhettan_filename = "manhettan_plot.jpeg",
- test = "CHI2", min_expected_cut_off = -1)
+ makePlot = FALSE, test = "CHI2",
+ min_expected_cut_off = -1)
}
\arguments{
\item{data}{Genotype data for analysis. Object of class
@@ -26,9 +25,8 @@
chisq's for all tests are stored. Object of class
\code{\link{logical}}}
- \item{manhettan_filename}{File name where manhettan plot
- will be saved after analysis. Object of class
- \code{\link{character}}}
+ \item{makePlot}{whether the Manhattan-type plot should be
+ produced (TRUE or FALSE)}
\item{test}{Name of the test to be performed. Available
tests are "CHI2", "YATES" (chi2 with Yates correction),
@@ -46,6 +44,10 @@
A list is returned.
}
\description{
+ Detecting rare recessive and compound heterozygote
+ alleles in genome wide association.
+}
+\details{
The function is an inplementation of the method aimed to
detect a gene-phenotype association caused by recessive
and compound heterozygote genotype states of multiple
@@ -56,8 +58,7 @@
van Duijn, Albert Hofman, Andre G. Uitterlinden, Yurii S.
Aulchenko, and Manfred Kayser. PLoS ONE 6(11): e28145.
doi:10.1371/journal.pone.0028145
-}
-\details{
+
The three tests are implemented: Pearson's chi-square
test, Pearson's chi-square test with Yates correction,
Fisher exact test. In case when the input parameter
@@ -71,7 +72,8 @@
}
\examples{
data(srdta)
-chis2_nocorrection <- cocohet(data=srdta
+chis2_nocorrection <- cocohet(data=gtdata(srdta),
+trait=phdata(srdta)$bt, window=3, test="CHI2")
}
\author{
Maksim Struchalin
More information about the Genabel-commits
mailing list