[Genabel-commits] r950 - pkg/GenABEL/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 31 09:14:35 CEST 2012
Author: nd_0001
Date: 2012-08-31 09:14:34 +0200 (Fri, 31 Aug 2012)
New Revision: 950
Modified:
pkg/GenABEL/R/convert.snp.text.R
pkg/GenABEL/R/descriptives.marker.R
pkg/GenABEL/R/egscore.R
pkg/GenABEL/R/egscore.old.R
pkg/GenABEL/R/estlambda.R
pkg/GenABEL/R/export.merlin.R
pkg/GenABEL/R/findRelatives.R
pkg/GenABEL/R/generateOffspring.R
pkg/GenABEL/R/getLogLikelihoodGivenRelation.R
pkg/GenABEL/R/grammar.R
pkg/GenABEL/R/grammar.old.R
pkg/GenABEL/R/load.gwaa.data.R
Log:
Some little changes in GenABEL
Modified: pkg/GenABEL/R/convert.snp.text.R
===================================================================
--- pkg/GenABEL/R/convert.snp.text.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/convert.snp.text.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -3,7 +3,7 @@
ifile <- file(infile,"r")
ofile <- file(outfile,"w")
header <- scan(file=ifile,what=character(),nlines=1,quiet=TRUE)
- vver <- grep(x=header,pat="version")
+ vver <- grep(x=header,pattern="version")
if (length(vver)>0) {ver <- as.numeric(header[vver+1]);} else {ver <- 0;}
if (is.na(ver)) warning("Incorrect data format version number")
if (ver > 0) {ids <- scan(file=ifile,what=character(),nlines=1,quiet=TRUE);}
Modified: pkg/GenABEL/R/descriptives.marker.R
===================================================================
--- pkg/GenABEL/R/descriptives.marker.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/descriptives.marker.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -29,10 +29,10 @@
mhet <- mean(s[,"Het"],na.rm=T)
sdhet <- sd(s[,"Het"],na.rm=T)
rm(s);gc(verbose=FALSE)
- out[["Minor allele frequency distribution"]] <- catable(maf,cat=mafc,digits=digits)
- out[["Cumulative distr. of number of SNPs out of HWE, at different alpha"]] <- catable(hwe,cat=hwec,cumulative=TRUE,digits=digits)
- out[["Distribution of proportion of successful genotypes (per person)"]] <- catable(ill,cat=idcc,digits=digits)
- out[["Distribution of proportion of successful genotypes (per SNP)"]] <- catable(cll,cat=snpc,digits=digits)
+ out[["Minor allele frequency distribution"]] <- catable(maf,categories=mafc,digits=digits)
+ out[["Cumulative distr. of number of SNPs out of HWE, at different alpha"]] <- catable(hwe,categories=hwec,cumulative=TRUE,digits=digits)
+ out[["Distribution of proportion of successful genotypes (per person)"]] <- catable(ill,categories=idcc,digits=digits)
+ out[["Distribution of proportion of successful genotypes (per SNP)"]] <- catable(cll,categories=snpc,digits=digits)
out[["Mean heterozygosity for a SNP"]] <- smh
out[["Standard deviation of the mean heterozygosity for a SNP"]] <- ssdh
out[["Mean heterozygosity for a person"]] <- mhet
Modified: pkg/GenABEL/R/egscore.R
===================================================================
--- pkg/GenABEL/R/egscore.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/egscore.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -97,7 +97,7 @@
lambda$se <- NA
} else {
if (lenn<100) warning("Number of observations < 100, Lambda estimate is unreliable")
- lambda <- estlambda(chi2.1df,plot=FALSE,prop=propPs)
+ lambda <- estlambda(chi2.1df,plot=FALSE,proportion=propPs)
if (lambda$estimate<1.0 && clambda==TRUE) {
warning("Lambda estimated < 1, set to 1")
lambda$estimate <- 1.0
@@ -147,9 +147,9 @@
Pc1df <- pr.c1df/times
Pc1df <- replace(Pc1df,(Pc1df==0),1/(1+times))
} else {
- P1df <- pchisq(chi2.1df,1,lower=F)
- P2df <- pchisq(chi2.2df,actdf,lower=F)
- Pc1df <- pchisq(chi2.c1df,1,lower=F)
+ P1df <- pchisq(chi2.1df,1,lower.tail=FALSE)
+ P2df <- pchisq(chi2.2df,actdf,lower.tail=FALSE)
+ Pc1df <- pchisq(chi2.c1df,1,lower.tail=FALSE)
}
#out$lambda <- lambda
#out$effB <- effB
Modified: pkg/GenABEL/R/egscore.old.R
===================================================================
--- pkg/GenABEL/R/egscore.old.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/egscore.old.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -95,7 +95,7 @@
lambda$se <- NA
} else {
if (lenn<100) warning("Number of observations < 100, Lambda estimate is unreliable")
- lambda <- estlambda(chi2.1df,plot=FALSE,prop=propPs)
+ lambda <- estlambda(chi2.1df,plot=FALSE,proportion=propPs)
if (lambda$estimate<1.0 && clambda==TRUE) {
warning("Lambda estimated < 1, set to 1")
lambda$estimate <- 1.0
@@ -145,9 +145,9 @@
out$Pc1df <- pr.c1df/times
out$Pc1df <- replace(out$Pc1df,(out$Pc1df==0),1/(1+times))
} else {
- out$P1df <- pchisq(chi2.1df,1,lower=F)
- out$P2df <- pchisq(chi2.2df,actdf,lower=F)
- out$Pc1df <- pchisq(chi2.c1df,1,lower=F)
+ out$P1df <- pchisq(chi2.1df,1,lower.tail=FALSE)
+ out$P2df <- pchisq(chi2.2df,actdf,lower.tail=FALSE)
+ out$Pc1df <- pchisq(chi2.c1df,1,lower.tail=FALSE)
}
out$lambda <- lambda
out$effB <- effB
Modified: pkg/GenABEL/R/estlambda.R
===================================================================
--- pkg/GenABEL/R/estlambda.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/estlambda.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -54,7 +54,7 @@
# warning(paste("Some probabilities < 1.e-16; set to 1.e-16"))
# data[lt16] <- 1.e-16
# }
- data <- qchisq(data,1,low=FALSE)
+ data <- qchisq(data,1,lower.tail=FALSE)
}
if (filter)
{
Modified: pkg/GenABEL/R/export.merlin.R
===================================================================
--- pkg/GenABEL/R/export.merlin.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/export.merlin.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -72,21 +72,21 @@
inf <- rbind(inf0,inf)
}
if (!is.null(datafile)) {
- write.table(inf,file=datafile,col.n=FALSE,row.n=FALSE,quote=FALSE)
+ write.table(inf,file=datafile,col.names=FALSE,row.names=FALSE,quote=FALSE)
}
if (format=="merlin") {
map <- data.frame(chromosome=as.character(data at gtdata@chromosome),markername=data at gtdata@snpnames,position=data at gtdata@map)
- write.table(map,file=mapfile,col.n=TRUE,row.n=FALSE,quote=FALSE)
+ write.table(map,file=mapfile,col.names=TRUE,row.names=FALSE,quote=FALSE)
} else if (format=="plink") {
map <- data.frame(chromosome=as.character(chromosome(data)),
markername=snpnames(data),gpos=0,position=map(data))
- write.table(map,file=mapfile,col.n=FALSE,row.n=FALSE,quote=FALSE)
+ write.table(map,file=mapfile,col.names=FALSE,row.names=FALSE,quote=FALSE)
} else {
stop("non-formalized format")
}
if (extendedmap) {
map <- data.frame(chromosome=as.character(data at gtdata@chromosome),markername=data at gtdata@snpnames,position=data at gtdata@map,strand=as.character(data at gtdata@strand),coding=as.character(data at gtdata@coding))
- write.table(map,file=paste(mapfile,".ext",sep=""),col.n=TRUE,row.n=FALSE,quote=FALSE)
+ write.table(map,file=paste(mapfile,".ext",sep=""),col.names=TRUE,row.names=FALSE,quote=FALSE)
}
}
@@ -108,7 +108,7 @@
} else {
x <- data.frame(seq(fromid,toid),ids,rep(0,nids),rep(0,nids),sx,x)
}
- write.table(x,file=pedfile,col.n=FALSE,row.n=FALSE,quote=FALSE,append=append,sep=" ")
+ write.table(x,file=pedfile,col.names=FALSE,row.names=FALSE,quote=FALSE,append=append,sep=" ")
}
dump.piece.New <- function(data,fromid,toid,traits,pedfile,append,format="merlin") {
Modified: pkg/GenABEL/R/findRelatives.R
===================================================================
--- pkg/GenABEL/R/findRelatives.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/findRelatives.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -259,7 +259,7 @@
}
if (!quiet) cat("\n")
#print(out)
- meiMtmp <- apply(out,MAR=1,FUN=function(x){return(which.max(x))})
+ meiMtmp <- apply(out,MARGIN=1,FUN=function(x){return(which.max(x))})
#print(meiMtmp)
if (is.null(vsIDs)) {
meiM <- matrix(ncol=dim(gtdata)[1],nrow=dim(gtdata)[1])
@@ -274,9 +274,9 @@
colnames(meiM) <- idnam[notVsIDs]
}
#print(meiM)
- firstBestLik <- apply(out,MAR=1,FUN=function(x){return(max(x))})
- which_firstBestLik <- apply(out,MAR=1,FUN=function(x){return(which.max(x))})
- secondBestLik <- apply(out,MAR=1,FUN=function(x){x[order(x,decreasing=TRUE)[2]]})
+ firstBestLik <- apply(out,MARGIN=1,FUN=function(x){return(max(x))})
+ which_firstBestLik <- apply(out,MARGIN=1,FUN=function(x){return(which.max(x))})
+ secondBestLik <- apply(out,MARGIN=1,FUN=function(x){x[order(x,decreasing=TRUE)[2]]})
cndX <- !(exp(firstBestLik-secondBestLik)>OddsVsNextBest
& exp(firstBestLik)>OddsVsNull & (which_firstBestLik != (length(meiTab)-1)))
if (is.null(vsIDs)) {
@@ -292,7 +292,7 @@
out <- data.frame(outN1,outN2,out,stringsAsFactors=FALSE)
names(out) <- c("id1","id2",names(meiTab))
tmp <- guess; diag(tmp) <- NA;
- todrop <- apply(tmp,MAR=1,FUN=function(x){return(all(is.na(x)))})
+ todrop <- apply(tmp,MARGIN=1,FUN=function(x){return(all(is.na(x)))})
compressedGuess <- guess[!todrop,!todrop]
finalOut <- list(call=match.call(),profile=out,estimatedNmeioses=meiM,
firstBestLik=firstBestLik,
Modified: pkg/GenABEL/R/generateOffspring.R
===================================================================
--- pkg/GenABEL/R/generateOffspring.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/generateOffspring.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -59,6 +59,6 @@
}
# generate offspring genotypes
g1g2 <- t(matrix(c(g1+1,g2+1),ncol=2))
- res <- apply(g1g2,MAR=2,FUN=genSOG)
+ res <- apply(g1g2,MARGIN=2,FUN=genSOG)
res
}
Modified: pkg/GenABEL/R/getLogLikelihoodGivenRelation.R
===================================================================
--- pkg/GenABEL/R/getLogLikelihoodGivenRelation.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/getLogLikelihoodGivenRelation.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -68,7 +68,7 @@
bGenotype1[2,] * TransitionMatrix[2,3,] +
bGenotype1[3,] * TransitionMatrix[3,3,]
x <- mm * bGenotype2
- x <- apply(x,FUN=sum,MAR=2)
+ x <- apply(x,FUN=sum,MARGIN=2)
out <- list(locusLik = x, logLik = sum(log(x)))
out
}
Modified: pkg/GenABEL/R/grammar.R
===================================================================
--- pkg/GenABEL/R/grammar.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/grammar.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -88,18 +88,18 @@
out at results[,"chi2.1df"] <- out[,"chi2.1df"]/polyObject$grammarGamma$Test
out at results[,"effB"] <- out[,"chi2.1df"]/polyObject$grammarGamma$Beta
# recompute p-values
- out at results[,"P1df"] <- pchisq(out[,"chi2.1df"],df=1,low=FALSE)
+ out at results[,"P1df"] <- pchisq(out[,"chi2.1df"],df=1,lower.tail=FALSE)
# recompute Lambda
- out at lambda <- estlambda(out[,"chi2.1df"],plot=FALSE,prop=propPs)
+ out at lambda <- estlambda(out[,"chi2.1df"],plot=FALSE,proportion=propPs)
if (out at lambda$estimate <= 1) {
warning("estimate of Lambda < 1, constraining to 1")
out at lambda$estimate <- 1.0
out at lambda$se <- NA
}
} else if (method == "gc") {
- out <- qtscore(polyObject$pgres,data=data,clambda=FALSE, prop=propPs, ... )
+ out <- qtscore(polyObject$pgres,data=data,clambda=FALSE, propPs=propPs, ... )
} else if (method == "raw") {
- out <- qtscore(polyObject$pgres,data=data,clambda=TRUE, prop=propPs, ... )
+ out <- qtscore(polyObject$pgres,data=data,clambda=TRUE, propPs=propPs, ... )
} else {
stop("method should be one of 'gamma','gc', or 'raw'")
}
Modified: pkg/GenABEL/R/grammar.old.R
===================================================================
--- pkg/GenABEL/R/grammar.old.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/grammar.old.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -73,7 +73,7 @@
lambda$se <- NA
} else {
if (lenn<100) warning("Number of observations < 100, Lambda estimate is unreliable")
- lambda <- estlambda(chi2.1df,plot=FALSE,prop=propPs)
+ lambda <- estlambda(chi2.1df,plot=FALSE,proportion=propPs)
if (lambda$estimate<1.0 && clambda==TRUE) {
warning("Lambda estimated < 1, set to 1")
lambda$estimate <- 1.0
@@ -123,9 +123,9 @@
out$Pc1df <- pr.c1df/times
# out$Pc1df <- replace(out$Pc1df,(out$Pc1df==0),1/(1+times))
} else {
- out$P1df <- pchisq(chi2.1df,1,lower=F)
-# out$P2df <- pchisq(chi2.2df,actdf,lower=F)
- out$Pc1df <- pchisq(chi2.c1df,1,lower=F)
+ out$P1df <- pchisq(chi2.1df,1,lower.tail=F)
+# out$P2df <- pchisq(chi2.2df,actdf,lower.tail=F)
+ out$Pc1df <- pchisq(chi2.c1df,1,lower.tail=F)
}
out$lambda <- lambda
out$effB <- effB
Modified: pkg/GenABEL/R/load.gwaa.data.R
===================================================================
--- pkg/GenABEL/R/load.gwaa.data.R 2012-08-30 09:43:00 UTC (rev 949)
+++ pkg/GenABEL/R/load.gwaa.data.R 2012-08-31 07:14:34 UTC (rev 950)
@@ -35,7 +35,7 @@
# read in genotypic data
ifile <- file(genofile,"r")
header <- scan(file=ifile,what=character(),nlines=1,quiet=TRUE)
- vver <- grep(x=header,pat="version")
+ vver <- grep(x=header,pattern="version")
if (length(vver)>0) {ver <- as.numeric(header[vver+1]);} else {ver <- 0;}
if (is.na(ver)) warning("Incorrect data format version number")
if (ver > 0) {ids <- scan(file=ifile,what=character(),nlines=1,quiet=TRUE);}
More information about the Genabel-commits
mailing list