[Genabel-commits] r782 - in pkg/GenABEL: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 11 21:37:12 CEST 2011
Author: yurii
Date: 2011-09-11 21:37:12 +0200 (Sun, 11 Sep 2011)
New Revision: 782
Modified:
pkg/GenABEL/CHANGES.LOG
pkg/GenABEL/DESCRIPTION
pkg/GenABEL/R/estlambda.R
pkg/GenABEL/R/merge.snp.data.R
pkg/GenABEL/R/mlreg.R
pkg/GenABEL/R/mlreg.p.R
pkg/GenABEL/R/polygenic.R
pkg/GenABEL/R/ss.R
pkg/GenABEL/R/zzz.R
pkg/GenABEL/man/estlambda.Rd
Log:
Deal with exceptional situation in merge.snp.data / monomorphic part;
added case of no overlap in SNPs (skip monomorphic staff then)
Modifications in 'estlambda': plot=FALSE by default, added option
to estimate Lambda with median method; added option 'filter' for
filtering SNPs with 0-statistics
Modified: pkg/GenABEL/CHANGES.LOG
===================================================================
--- pkg/GenABEL/CHANGES.LOG 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/CHANGES.LOG 2011-09-11 19:37:12 UTC (rev 782)
@@ -1,3 +1,12 @@
+*** v. 1.7-0 (2011.09.10)
+
+Deal with exceptional situation in merge.snp.data / monomorphic part;
+added case of no overlap in SNPs (skip monomorphic staff then)
+
+Modifications in 'estlambda': plot=FALSE by default, added option
+to estimate Lambda with median method; added option 'filter' for
+filtering SNPs with 0-statistics
+
*** v. 1.6-9 (2011.08.30)
Change the name of 'maintainer' to genabel.project at gmail.com to avoid spamming
Modified: pkg/GenABEL/DESCRIPTION
===================================================================
--- pkg/GenABEL/DESCRIPTION 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/DESCRIPTION 2011-09-11 19:37:12 UTC (rev 782)
@@ -1,8 +1,8 @@
Package: GenABEL
Type: Package
Title: genome-wide SNP association analysis
-Version: 1.6-9
-Date: 2011-08-30
+Version: 1.7-0
+Date: 2011-09-10
Author: GenABEL developers
Maintainer: GenABEL developers <genabel.project at gmail.com>
Depends: R (>= 2.10.0), methods, MASS
Modified: pkg/GenABEL/R/estlambda.R
===================================================================
--- pkg/GenABEL/R/estlambda.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/estlambda.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -1,5 +1,5 @@
"estlambda" <-
-function(data,plot=TRUE,proportion=1.0, ...) {
+ function(data,plot = FALSE, proportion=1.0, method="regression", filter=TRUE,...) {
data <- data[which(!is.na(data))]
if (proportion>1.0 || proportion<=0) stop("proportion argument should be greater then zero and less than or equal to one")
ntp <- round(proportion*length(data))
@@ -18,6 +18,10 @@
# }
data <- qchisq(data,1,low=FALSE)
}
+ if (filter)
+ {
+ data[which(abs(data)<1e-8)] <- NA
+ }
data <- sort(data)
ppoi <- ppoints(data)
ppoi <- sort(qchisq(1-ppoi,1))
@@ -25,7 +29,6 @@
ppoi <- ppoi[1:ntp]
# s <- summary(lm(data~offset(ppoi)))$coeff
# bug fix thanks to Franz Quehenberger
- s <- summary(lm(data~0+ppoi))$coeff
if (plot) {
lim <- c(0,max(data,ppoi,na.rm=T))
# plot(ppoi,data,xlim=lim,ylim=lim,xlab="Expected",ylab="Observed", ...)
@@ -34,7 +37,15 @@
abline(a=0,b=(s[1,1]),col="red")
}
out <- list()
- out$estimate <- s[1,1]
- out$se <- s[1,2]
+ if (method=="regression") {
+ s <- summary(lm(data~0+ppoi))$coeff
+ out$estimate <- s[1,1]
+ out$se <- s[1,2]
+ } else if (method=="median") {
+ out$estimate <- median(data,na.rm=TRUE)/qchisq(0.5,1)
+ out$se <- NA
+ } else {
+ stop("'method' should be either 'regression' or 'median'!")
+ }
out
}
Modified: pkg/GenABEL/R/merge.snp.data.R
===================================================================
--- pkg/GenABEL/R/merge.snp.data.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/merge.snp.data.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -141,6 +141,14 @@
}
#x_filtered_mono_x_not_mono_y <- which(monos(coding(x)[which_snp_intersect_in_x]) &
# !monos(coding(y)[which_snp_intersect_in_y]))
+#print(which_snp_intersect_in_y)
+#print(coding(y))
+#print(coding(y)[which_snp_intersect_in_y])
+#print(!monos(coding(y)[which_snp_intersect_in_y]))
+#
+# do this only if there is any intersection between SNPs
+#
+if (length(which_snp_intersect_in_y)>0 || length(which_snp_intersect_in_x)>0) {
if_mono_x_not_mono_y <- monos(coding(x)[which_snp_intersect_in_x]) &
!monos(coding(y)[which_snp_intersect_in_y])
which_mono_x_not_mono_y <- which_snp_intersect_in_x[if_mono_x_not_mono_y]
@@ -199,6 +207,7 @@
coding(y)[iii] <- newC
# print(newC)
}
+}
##### END taking care of monomorphics ###################
if(intersected_snps_only)
@@ -231,7 +240,8 @@
# stop("For intersected SNPs several values in chromosome vector for x is not concur with values in y. Check your data sets.\n")
}
-if(length(levels(chromosome_logic_vec_factor)) == 0)
+anyIntersectingSNPs <- !(length(levels(chromosome_logic_vec_factor)) == 0)
+if(!anyIntersectingSNPs)
{
cat("There are no intersecting SNPs\n")
if(intersected_snps_only) return(list(data=NULL, id=NULL, snp=NULL))
Modified: pkg/GenABEL/R/mlreg.R
===================================================================
--- pkg/GenABEL/R/mlreg.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/mlreg.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -71,7 +71,7 @@
if (gtdata at nids != n)
stop("incompatible dimensions")
if (any(gtdata at chromosome == "X") & !any(colnames(x)=="sex"))
- warning("You analysed X chromosome without adjhstment for sex")
+ warning("You analysed X chromosome without adjustment for sex")
if (ttype==1)
chi2 <- .C("linreg_gwaa",as.double(y),as.double(x),as.raw(gtdata at gtps),as.integer(n),
as.integer(ncol(x)),as.integer(gtdata at nsnps),as.integer(gtm),chi2=double(3*gtdata at nsnps))$chi2
Modified: pkg/GenABEL/R/mlreg.p.R
===================================================================
--- pkg/GenABEL/R/mlreg.p.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/mlreg.p.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -73,7 +73,7 @@
if (gtdata at nids != n)
stop("incompatible dimensions")
if (any(gtdata at chromosome == "X") & !any(colnames(x)=="sex"))
- warning("You analysed X chromosome without adjhstment for sex")
+ warning("You analysed X chromosome without adjustment for sex")
if (ttype==1)
chi2 <- .C("linreg_gwaa",as.double(y),as.double(x),as.raw(gtdata at gtps),as.integer(n),
as.integer(ncol(x)),as.integer(gtdata at nsnps),as.integer(gtm),chi2=double(3*gtdata at nsnps))$chi2
Modified: pkg/GenABEL/R/polygenic.R
===================================================================
--- pkg/GenABEL/R/polygenic.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/polygenic.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -71,7 +71,7 @@
#' from previous regression, these are expected to change little from the
#' initial estimate. The default value of 1000 proved to work rather well under a
#' range of conditions.
-#' @param quiet If FALSE (default), details of optimisation process are reported.
+#' @param quiet If FALSE (default), details of optimisation process are reported
#' @param steptol steptal parameter of "nlm"
#' @param gradtol gradtol parameter of "nlm"
#' @param optimbou fixed effects boundary scale parameter for 'optim'
@@ -168,7 +168,7 @@
#'
"polygenic" <-
function(formula,kinship.matrix,data,fixh2,starth2=0.3,trait.type="gaussian",
- opt.method="nlm",scaleh2=1,quiet=FALSE,
+ opt.method="nlm",scaleh2=1,quiet=FALSE,
steptol=1e-8, gradtol = 1e-8, optimbou = 8,
fglschecks=TRUE,maxnfgls=8,maxdiffgls=1e-4, patchBasedOnFGLS = TRUE,
llfun = "polylik_eigen", ...) {
Modified: pkg/GenABEL/R/ss.R
===================================================================
--- pkg/GenABEL/R/ss.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/ss.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -314,7 +314,7 @@
definition = function(x,i,j,drop)
{
res <- results(x)
- if (missing(j)) drop=FALSE
+ if (missing(j)) drop <- FALSE
return(res[i,j,drop=drop])
}
);
Modified: pkg/GenABEL/R/zzz.R
===================================================================
--- pkg/GenABEL/R/zzz.R 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/R/zzz.R 2011-09-11 19:37:12 UTC (rev 782)
@@ -1,6 +1,6 @@
.onLoad <- function(lib, pkg) {
- GenABEL.version <- "1.6-9"
- cat("GenABEL v.",GenABEL.version,"(August 30, 2011) loaded\n")
+ GenABEL.version <- "1.7-0"
+ cat("GenABEL v.",GenABEL.version,"(September 10, 2011) loaded\n")
# check for updates and news
address <- c(
Modified: pkg/GenABEL/man/estlambda.Rd
===================================================================
--- pkg/GenABEL/man/estlambda.Rd 2011-09-02 09:31:59 UTC (rev 781)
+++ pkg/GenABEL/man/estlambda.Rd 2011-09-11 19:37:12 UTC (rev 782)
@@ -7,7 +7,7 @@
visualise the distribution of P-values coming from other tests.
}
\usage{
-estlambda(data, plot = TRUE, proportion = 1.0, ... )
+estlambda(data, plot = FALSE, proportion = 1.0, method="regression", filter=TRUE, ... )
}
\arguments{
\item{data}{A vector of reals. If all are <=1, it is assumed that this
@@ -16,6 +16,10 @@
\item{plot}{Wether the prot should be presented}
\item{proportion}{The proportion of lowest P (Chi2) to be used when
estimating the inflation factor Lambda}
+ \item{method}{Either "regression" or "median", the method
+ to be used for Lambda estimation}
+ \item{filter}{if the test statistics with 0-value of chi-square
+ should be excluded prior to estimation of Lambda}
\item{...}{arguments passed to \code{\link{plot}} function}
}
%\details{
More information about the Genabel-commits
mailing list