diff -Nru r-other-mott-happy-2.1/DESCRIPTION r-other-mott-happy-2.4/DESCRIPTION --- r-other-mott-happy-2.1/DESCRIPTION 2008-03-07 10:20:49.000000000 +0000 +++ r-other-mott-happy-2.4/DESCRIPTION 2013-12-19 16:21:35.000000000 +0000 @@ -1,10 +1,11 @@ -Package: happy -Version: 2.1 -Date: 2008-01-31 -Title: Quantitative Trait Locus genetic analysis in Heterogeneous Stocks +Package: happy.hbrem +Version: 2.4 +Date: 2012-07-20 +Title: Quantitative Trait Locus genetic analysis in Heterogeneous + Stocks Author: Richard Mott Maintainer: Richard Mott -Depends: R (>= 2.6.0), g.data +Depends: R (>= 2.6.0), g.data, multicore Description: happy is an R interface into the HAPPY C package for fine-mapping Quantitative Trait Loci (QTL) in Heterogenous Stocks (HS). An HS is an advanced intercross between (usually @@ -14,6 +15,6 @@ probability of descent from each of the founders, at each locus position, but the happy packager allows a much richer range of models to be fit to the data. -License: GPL version 2 or newer +License: GPL (>=2) URL: http://www.r-project.org, http://www.well.ox.ac.uk/happy -Packaged: Fri Mar 7 10:20:49 2008; rmott +Packaged: 2013-12-19 16:21:35 UTC; rmott diff -Nru r-other-mott-happy-2.1/INDEX r-other-mott-happy-2.4/INDEX --- r-other-mott-happy-2.1/INDEX 2006-10-21 20:33:55.000000000 +0000 +++ r-other-mott-happy-2.4/INDEX 1970-01-01 00:00:00.000000000 +0000 @@ -1,18 +0,0 @@ -epistasis Analysis of Epistasis between Markers -gfit Fit a Gaussian Mixture Model to an object - returned by happy() -happy-internal Internal Happy Functions -happyplot Plotting functions for happy model fits -hdesign Extract design matrix or genotypes for a - specific marker interval from a happy object -hfit Fit a model to an object returned by happy() -introduction Quantitative Trait Locus analysis in - Heterogeneous Stocks -mergelist Create an object descrbing how to merge strains - together -mergematrices Construct matrices used to merge together - founder strains -mergeprepare Perform tests to determine whether individual - polymorphisms could have given rise to a QTL -save.genome Save HAPPY design matrices and genotypes to - disk for rapid reloading diff -Nru r-other-mott-happy-2.1/NAMESPACE r-other-mott-happy-2.4/NAMESPACE --- r-other-mott-happy-2.1/NAMESPACE 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/NAMESPACE 2013-12-19 14:26:44.000000000 +0000 @@ -0,0 +1,4 @@ +exportPattern("^[^\\.]") +exportPattern("^[[:alpha:]]+") +useDynLib(happy.hbrem) + diff -Nru r-other-mott-happy-2.1/R/happy.R r-other-mott-happy-2.4/R/happy.R --- r-other-mott-happy-2.1/R/happy.R 2008-02-01 16:46:05.000000000 +0000 +++ r-other-mott-happy-2.4/R/happy.R 2013-12-19 16:12:52.000000000 +0000 @@ -1,8 +1,8 @@ -.packageName <- "happy_2.0.1" +.packageName <- "happy.hbrem_2.4" library(MASS) library(g.data) - +library(multicore) # C interface to read in .data and .alleles files, perform DP and create a happy object # # a happy object is a list with the following attributes: @@ -14,56 +14,75 @@ # happy$handle integer handle which maps the R happy object to the corresponding C QTL object -happy<- function( datafile, allelesfile, generations=200, - standardise=FALSE, phase="unknown", +happy <- function( datafile, allelesfile, generations=200, phase="unknown", file.format="happy", missing.code="NA", do.dp=TRUE, min.dist=1.0e-5, mapfile=NULL, ancestryfile=NULL, haploid=FALSE ) { gen <- as.numeric(generations)+0 if ( phase=="estimate" ) file.format <- "ped" - do.dp = ifelse( do.dp, 1, 0 ) - hap <- ifelse( haploid, 1, 0 ); - h <- .Call( "happy", datafile, allelesfile, gen, phase, file.format, missing.code, do.dp, min.dist=min.dist, haploid=hap, ancestryfile=ancestryfile, PACKAGE="happy" ) - if ( standardise ) { - m <- mean(h$phenotype) - s <- sd(h$phenotype) - h$phenotype <- (h$phenotype-m)/s - cat('\nPhenotypes have been standardised to mean 0, variance 1\n') - } + + h <- .Call( "happy", datafile, allelesfile, gen, phase, file.format, missing.code, do.dp=as.integer(do.dp), min.dist=min.dist, haploid=as.integer(haploid), ancestryfile=ancestryfile, PACKAGE="happy.hbrem" ) h$phase <- phase h$haploid <- haploid - h$nam <- make.names(h$strains) - nam2 <- c() - ns <- length(h$nam) - for(i in 1:ns) { - for(j in 1:i) { - nam2 <- c( nam2, paste(h$nam[i], ".", h$nam[j], sep="")) - } - } - h$nam2 <- nam2 - - nam3 <- c() - for(i in 1:ns) { - for(j in 1:ns) { - nam3 <- c( nam3, paste(h$nam[i], ".", h$nam[j], sep="")) - } - } - h$nam3 <- nam3 + strain.names <- make.names(h$strains) + num.strains <- length(h$strains) + h$names.additive <- strain.names + + diplotype.names <- matrix(kronecker(strain.names, strain.names, paste, sep="."), nrow=num.strains) + h$names.full.symmetric <- c( diag(diplotype.names), diplotype.names[upper.tri(diplotype.names, diag=FALSE)]) + h$names.full.asymmetric <- c( t(diplotype.names) ) # assumes row major order, ie, (row1, row2, etc), in C object + + # for backwards compatibility + h$nam <- strain.names + h$nam2 <- h$names.full.symmetric + h$nam3 <- h$names.full.asymmetric + h$bp = h$map + if ( ! is.null(mapfile)) { - map <- read.table(mapfile, "\t", header=TRUE) + map <- read.delim(mapfile) if ( !is.null(map$marker) && ! is.null(map$bp)) { hbp <- data.frame(markers=h$markers,bp=rep(NA,length(h$markers))) h$bp <- map$bp[match(hbp$markers,map$marker)] } else - warning( "incorrect column names found in mapfile ", mapfile , "\n") + stop( "incorrect column names found in mapfile ", mapfile , "\n") } - + + h$additive = list() + nm = length(h$markers)-1 + h$additive$genome <- data.frame( + marker = I(as.character(h$markers))[1:nm], + map = as.numeric(h$map)[1:nm], + bp = as.numeric(h$bp)[1:nm], + chromosome = I(as.character(h$chromosome))[1:nm]) + + h$full = list() + h$full$genome <- data.frame( + marker = I(as.character(h$markers))[1:nm], + map = as.numeric(h$map)[1:nm], + bp = as.numeric(h$bp)[1:nm], + chromosome = I(as.character(h$chromosome))[1:nm]) + + mother = vector( mode="character", length=length(h$mother)) + w = which(h$mother>0) + mother[w] = h$subjects[h$mother[w]] + w = which(h$mother<=0) + mother[w] = NA + h$mother = mother + + father = vector( mode="character", length=length(h$father)) + w = which(h$father>0) + father[w] = h$subjects[h$father[w]] + w = which(h$father<=0) + father[w] = NA + h$father = father + return(h) } + happy.matrices <- function( h ) { if ( class(h) == "happy" ) { @@ -127,12 +146,13 @@ else { # data are in C memory handle <- as.numeric(h$handle)+0 if ( h$haploid ) { - d <- .Call( "haploid_happydesign", handle, marker, PACKAGE="happy") + d <- .Call( "haploid_happydesign", handle, marker, PACKAGE="happy.hbrem") if ( ! is.null(d) ) - colnames(d) <- h$nam + if ( model == 'additive' ) colnames(d) <- h$nam + else colnames(d) <- h$nam2 } else { - d <- .Call( "happydesign", handle, marker, model, PACKAGE="happy") + d <- .Call( "happydesign", handle, marker, model, PACKAGE="happy.hbrem") if ( ! is.null(d) ) { if ( model == 'additive' ) colnames(d) <- h$nam else colnames(d) <- h$nam2 @@ -170,12 +190,12 @@ nm <- length(h$markers)-1 r <- vector("numeric",length=nm) for(m in 1:nm) { - x <- .Call( "happynonrecomb", handle, m, PACKAGE="happy") + x <- .Call( "happynonrecomb", handle, m, PACKAGE="happy.hbrem") r[m] <- 0.5*mean(x) } } else { - r <- .Call( "happynonrecomb", handle, marker, PACKAGE="happy") + r <- .Call( "happynonrecomb", handle, marker, PACKAGE="happy.hbrem") r <- 0.5*r if ( do.mean ) r <- mean(r) @@ -187,22 +207,30 @@ hprob <- function( h, marker=NULL ) { handle <- as.numeric(h$handle)+0 - p <- .Call( "happyprobs", handle, marker, PACKAGE="happy") + p <- .Call( "happyprobs", handle, marker, PACKAGE="happy.hbrem") colnames(p) <- h$nam2 return(p) } -hprob2 <- function( h, marker=NULL ) { +hprob2 <- function( h, marker=NULL, symmetrize=FALSE ) { handle <- as.numeric(h$handle)+0 - p <- .Call( "happyprobs2", handle, marker, PACKAGE="happy") - colnames(p) <- h$nam3 + if ( symmetrize==TRUE ) + symmetrize=1 + else + symmetrize=0 + + p <- .Call( "happyprobs2", handle, marker, symmetrize, PACKAGE="happy.hbrem") + if ( symmetrize==1 ) + colnames(p) <- h$nam2 + else + colnames(p) <- h$nam3 return(p) } h.sum.prob2 <- function( h, marker=NULL ) { # for Caroline - the sum of squares of the probabilities handle <- as.numeric(h$handle)+0 if ( ! is.null( marker ) ) { - p <- .Call( "happyprobs2", handle, marker, PACKAGE="happy") + p <- .Call( "happyprobs2", handle, marker, PACKAGE="happy.hbrem") p2 <- apply( p*p, 1, sum ) return(p2) } @@ -210,7 +238,7 @@ nm <- length(h$markers)-1 mat <- matrix( nrow=length(h$subjects), ncol=nm ) for(i in 1:nm ) { - p <- .Call( "happyprobs2", handle, i, PACKAGE="happy") + p <- .Call( "happyprobs2", handle, i, PACKAGE="happy.hbrem") mat[,i] <- apply( p*p, 1, sum ) } rownames(mat) <- h$subjects @@ -221,16 +249,23 @@ hgenotype <- function( h, marker=NULL, collapse=FALSE, sep="" ) { - handle <- as.numeric(h$handle)+0 - g <- .Call( "happygenotype", handle, marker, PACKAGE="happy") - if ( collapse ) { - y <- paste( g[,1], g[,2], sep=sep ) - g <- ifelse ( y == "NANA", NA, y ) - } - else { - colnames(g) <- c("allele1", "allele2") - } - return(g) + + if ( class(h) == "happy" ) { + handle <- as.numeric(h$handle)+0 + g <- .Call( "happygenotype", handle, marker, PACKAGE="happy.hbrem") + } + else if ( class(h) == "happy.genome" ) { # delayed data package + loaded.markers <- load.markers( h, c(marker), model="genotype" ) + g <- loaded.markers[[1]] + } + if ( collapse ) { + y <- paste( g[,1], g[,2], sep=sep ) + g <- ifelse ( y == "NANA", NA, y ) + } + else { + colnames(g) <- c("allele1", "allele2") + } + return(g) } @@ -243,13 +278,28 @@ # return value is a table giving the logP values for each marker tested hfit <- function( h, markers=NULL, model='additive', mergematrix=NULL, covariatematrix=NULL, verbose=FALSE, phenotype=NULL, family='gaussian', permute=0 ) { - map <- h$map - if ( is.null(markers) ) - markers <- h$markers[1:length(h$markers)-1] + if ( class(h) == "happy.genome" ) { + if ( !is.null( h[[model]] ) ) + map <- h[[model]]$map + if ( is.null(markers) ) { + nm <- length(h[[model]]$markers)-1 + markers <- h[[model]]$markers[1:nm] + } + if ( is.null(phenotype)) { + stop( "phenotype must be set\n") + } + } + else { + + map <- h$map + + if ( is.null(markers) ) + markers <- h$markers[1:length(h$markers)-1] - if ( is.null(phenotype) ) - phenotype = h$phenotype + if ( is.null(phenotype) ) + phenotype = h$phenotype + } if ( model == 'partial' || model == 'full' ) { lp <- matrix( ncol=8, nrow=length(map)-1) @@ -832,7 +882,7 @@ mergematrix <- mergematrices( h$strains, mergelist=mlist) # the matrices representing the merge if ( design ) { # return the design matrix marker <- prepmerge$interval[ind] - retval <- hdesign( h, marker, merge=mergematrix, model=model ) + retval <- hdesign( h, marker, mergematrix=mergematrix, model=model ) } else # return the mergematrix retval <- mergematrix @@ -1291,7 +1341,7 @@ ma1 <- markers1[m1] print(ma1) for( m2 in 1:(m1-1) ) { - results[r,] <- epistasispair( h, ma1, markers1[m2], merge1=merge1, merge2=merge2, mode=mode, verbose=verbose, d1=d[[m1]], d2=d[[m2]], main1=main[[m1]], main2=main[[m2]] ) + results[r,] <- epistasispair( h, ma1, markers1[m2], merge1=merge1, merge2=merge2, model=model, verbose=verbose, d1=d[[m1]], d2=d[[m2]], main1=main[[m1]], main2=main[[m2]] ) r <- r+1 } } @@ -1317,7 +1367,7 @@ } epistasispair<- function( h, marker1, marker2, merge1=NULL, merge2=NULL, model='additive', verbose=FALSE, d1=NULL, d2=NULL, main1=0, main2=0, family='gaussian' ) { - + logten <- log(10) @@ -1636,7 +1686,7 @@ # Genome Cache functions save.genome <- function ( gdir, sdir, prefix, chrs=NULL, - file.format="ped", mapfile=NULL, ancestryfile=NULL, generations=50, phase="unknown", haploid=FALSE ) { + file.format="ped", mapfile=NULL, ancestryfile=NULL, generations=50, phase="unknown", haploid=FALSE, mc.cores=1 ) { if ( is.null(chrs) ) @@ -1652,17 +1702,33 @@ dir.create(additive) genotype <- paste(sdir, "/genotype/", sep="") dir.create(genotype) + + + if ( ! require(multicore)) mc.cores = 1 + + if ( mc.cores <=1 ) { + lapply( chrs, save.happy.internal, gdir, prefix, file.format, ancestryfile, generations, mapfile, phase, haploid, additive, full, genotype) + } + else { + mclapply( chrs, save.happy.internal, gdir, prefix, file.format, ancestryfile, generations, mapfile, phase, haploid, additive, full, genotype, mc.cores=mc.cores) + } +} - for ( chr in chrs ) { - h <- happy( paste( gdir, chr, prefix, ".data", sep="" ), - paste( gdir, chr, prefix, ".alleles", sep="" ), - file.format=file.format, - ancestryfile=ancestryfile, - generations=generations, do.dp=TRUE, mapfile=mapfile, phase=phase, haploid=haploid ) +save.happy.internal <- function( chr, gdir, prefix, file.format, ancestryfile, generations, mapfile, phase, haploid, additive, full, genotype) { + h <- happy( paste( gdir, chr, prefix, ".data", sep="" ), + paste( gdir, chr, prefix, ".alleles", sep="" ), + file.format=file.format, + ancestryfile=ancestryfile, + generations=generations, do.dp=TRUE, mapfile=mapfile, phase=phase, haploid=haploid ) save.happy( h, chr, dir=additive, model="additive" ) if ( h$haploid == FALSE ) save.happy( h, chr, dir=full, model="full" ) save.happy( h, chr, dir=genotype, model="genotype" ) - } + delete.happy.cobject(h) +} + +delete.happy.cobject <- function(h) +{ + cat("delete.happy() called\n") } @@ -1684,91 +1750,106 @@ nm <- length(h$markers) else nm <- length(h$markers) -1 - +# markers.safe = as.character(h$markers[1:nm]) + markers.safe = make.names(h$markers[1:nm]) + + assign("markers", h$markers[1:nm], 2) + assign("markers.safe",markers.safe,2) assign("map", h$map[1:nm], 2 ) assign("chromosome", h$chromosome[1:nm], 2 ) assign("subjects", h$subjects, 2 ); assign("strains", h$strains, 2 ) + assign("haploid", h$haploid, 2) print ("saving strains") if ( !is.null(h$bp)) assign("bp", h$bp[1:nm], 2); if ( model == "genotype" ) { for( m in 1:nm ) { - assign(h$markers[m], hgenotype( h, m, collapse=TRUE ), 2) + assign(markers.safe[m], hgenotype( h, m, collapse=FALSE ), 2) } - g.data.save(ddp, compress=TRUE) + g.data.save(ddp) } else { for( m in 1:nm ) { - assign(h$markers[m], hdesign( h, m, model=model ), 2) + assign(markers.safe[m], hdesign( h, m, model=model ), 2) } - g.data.save(ddp, compress=TRUE) + g.data.save(ddp) } return(ddp) } -load.genome <- function (dir, + +load.genome <- function (sdir, use.X = TRUE, - chr = the.chromosomes(use.X=use.X), + chrs = the.chromosomes(use.X=use.X), + n.chr=NA, models=c("additive", "full", "genotype")) # CHANGED { g <- list() - old.subjects <- NULL - old.strains <- NULL - for (model in models) { - pkgs <- paste(dir, model, chr, sep = "/") # CHANGED - markers <- c() - chromosome <- c() - map <- c() - pkgname <- c() - bp <- c() - for (p in pkgs) { - chromosome <- c(chromosome, g.data.get("chromosome", - p)) - m <- g.data.get("markers", p) - markers <- c(markers, m) - map <- c(map, g.data.get("map", p)) - bp <- c(bp, g.data.get("bp", p)) - pkgname <- c(pkgname, rep(p, length(m))) - subjects <- g.data.get("subjects", p) - strains <- g.data.get("strains", p) - - if ( is.null(old.subjects)) { - old.subjects <- subjects - } - if ( subjects != old.subjects ) { - cat( "ERROR - subject names are inconsistent for chromosome ", chromosome , "\n") - stop( "FATAL HAPPY ERROR") - } - - if ( is.null(old.strains)) { - old.strains <- strains - } - if ( strains != old.strains ) { - cat( "ERROR - strain names are inconsistent for chromosome ", chromosome , "\n") - stop( "FATAL HAPPY ERROR") - } - } - genome <- data.frame( - marker = I(as.character(markers)), - map = as.numeric(map), - bp = as.numeric(bp), - ddp = I(as.character(pkgname)), - chromosome = I(as.character(chromosome))) - g[[model]] <- list( - genome = genome, - subjects = subjects, - strains = strains, - markers = as.character(genome$marker), - chromosome = as.character(genome$chromosome), - map = genome$map) - } - g$subjects <- g$genotype$subjects - g$strains <- g$additive$strains - g$markers <- g$genotype$markers - class(g) <- "happy.genome" - return(g) + old.subjects <- NULL + old.strains <- NULL + if ( is.integer(n.chr) ) + chrs = paste("chr", 1:n.chr, sep="") + + for (model in models) { + if ( file.exists( paste(sdir, model, sep = "/") )) { + pkgs <- paste(sdir, model, chrs, sep = "/") # CHANGED + pkgs = pkgs[file.exists(pkgs)] + markers <- c() + chromosome <- c() + map <- c() + pkgname <- c() + bp <- c() + for (p in pkgs) { +# chromosome <- c(chromosome, g.data.get("chromosome", p)) + chromosome <- c(chromosome, happy.load.data("chromosome", p)) + m <- happy.load.data("markers", p) + markers <- c(markers, m) + map <- c(map, happy.load.data("map", p)) + bp <- c(bp, happy.load.data("bp", p)) + pkgname <- c(pkgname, rep(p, length(m))) + subjects <- happy.load.data("subjects", p) + strains <- happy.load.data("strains", p) + + if ( is.null(old.subjects)) { + old.subjects <- subjects + } + if ( !identical(subjects,old.subjects )) { + cat( "ERROR - subject names are inconsistent for chromosome ", chromosome[1] , "\n", subjects, "\n", old.subjects, "\n") + stop( "FATAL HAPPY ERROR") + } + + if ( is.null(old.strains)) { + old.strains <- strains + } + if ( ! identical( strains, old.strains) ) { + cat( "ERROR - strain names are inconsistent for chromosome ", chromosome[1] , "\n", strains, "\n", old.strains, "\n") + stop( "FATAL HAPPY ERROR") + } + } + genome <- data.frame( + marker = I(as.character(markers)), + map = as.numeric(map), + bp = as.numeric(bp), + ddp = I(as.character(pkgname)), + chromosome = I(as.character(chromosome))) + g[[model]] <- list( + genome = genome, + subjects = subjects, + strains = strains, + markers = as.character(genome$marker), + chromosome = as.character(genome$chromosome), + map = genome$map) + } + } + g$subjects <- g$genotype$subjects + g$strains <- g$additive$strains + g$markers <- g$genotype$markers + g$chromosome <- g$genotype$chromosome + g$map <- g$genotype$map + class(g) <- "happy.genome" + return(g) } @@ -1781,6 +1862,7 @@ marker.list <- list() model.list <- list() + for(i in 1:length(markers)) { genome.model <- genome[[model[i]]] @@ -1792,9 +1874,12 @@ if ( length(rows) > 0 ) { r <- rows[1] + m <- as.character(genome.model$genome[r,"marker"]) +# m.names = make.names(m) + pkg <- as.character(genome.model$genome[r,"ddp"]) - marker.list[[m]] <- g.data.get( m, pkg) + marker.list[[m]] <- happy.load.data( m, pkg) ### model.list[[m]] <- model[i] } @@ -1805,3 +1890,38 @@ return( marker.list ) } +happy.load.data <- function (item, dir) # replaces calls to g.data.get, to make things backwards compatible. +{ + env <- new.env() + + # determine which version of g.data was used to save the data + filename.pre2009 <- file.path(dir, "data", paste(item, "RData", sep = ".")) + if (file.exists(filename.pre2009)) + { + load(filename.pre2009, env) + return ( get(item, envir = env) ) + } + + # assume 2009 version of g.data was used + filename.post2009 <- file.path(dir, + paste(gsub("([[:upper:]])", "@\\1", item), "RData", sep = ".") + ) + if (file.exists(filename.post2009)) + { + load(filename.post2009, env) + return ( get(item, envir = env ) ) + } + mm = make.names(item) + filename.make.names = file.path(dir,paste(gsub("([[:upper:]])", "@\\1", mm), "RData", sep = ".")) + + if (file.exists(filename.make.names)) + { + load(filename.make.names, env) + return ( get(mm, envir = env ) ) + } + stop("Could not find data for ", item, " in package ", dir) +} + + + + diff -Nru r-other-mott-happy-2.1/R/hbrem.R r-other-mott-happy-2.4/R/hbrem.R --- r-other-mott-happy-2.1/R/hbrem.R 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/R/hbrem.R 2013-12-19 14:26:44.000000000 +0000 @@ -0,0 +1,317 @@ +.packageName <- "happy.hbrem_2.2" + + +hbrem <- function( RX, HaploidInd, Ndip, Nstrain, Nind, Npost=2000, Nbin, Ry ) { + + brem <- .Call( "hbrem", RX, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin, Ry, PACKAGE="happy.hbrem" ) + + return(brem) +} + + +hbrem.true <- function( RX, HaploidInd, Ndip, Nstrain, Nind, Npost=2000, Nbin, Ry ) { + + brem.true <- .Call( "hbrem_true", RX, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin, Ry, PACKAGE="happy.hbrem" ) + + return(brem.true) +} + + +hbrem.locus <- function(m, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) { + + if ( class(g) == "condensed.happy" ) { + + cum.mark <- 0 + chr <- 0 + while ( m > cum.mark) { + chr = ( chr + 1 ) + cum.mark = ( cum.mark + length(g[[model]]$chr[[chr]]) ) + } + + cum.gen <- ( cum.mark - length(g[[model]]$chr[[chr]]) ) + chr.mark <- ( m - cum.gen ) + + d <- g[[model]]$chr[[chr]][chr.mark][[1]][[1]] + d.cc <- d[cc,] + + } else { + d <- hdesign(g, m, model=model) + d.cc <- d[cc,] + } + + hb <- hbrem(RX=d.cc, HaploidInd=HaploidInd, Ndip=Ndip, Nstrain=Nstrain, Nind=Nind, Npost=Npost, Nbin=Nbin, Ry=Ry) + + reg.full.lm <- lm(Ry ~ d.cc) + reg.null.lm <- lm(Ry ~ 1) + Ftest <- anova(reg.null.lm, reg.full.lm) + pval <- Ftest$"Pr(>F)"[2] + + # cat( m, -log10(pval), "\n" ) + return( c( -log10(pval), hb[[1]], hb[[2]], hb[[3]], hb[[4]], hb[[5]] )) +} + + +hbrem.perm.locus <- function(m, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) { + + if ( class(g) == "condensed.happy" ) { + + cum.mark <- 0 + chr <- 0 + while ( m > cum.mark) { + chr = ( chr + 1 ) + cum.mark = ( cum.mark + length(g[[model]]$chr[[chr]]) ) + } + + cum.gen <- ( cum.mark - length(g[[model]]$chr[[chr]]) ) + chr.mark <- ( m - cum.gen ) + + d <- g[[model]]$chr[[chr]][chr.mark][[1]][[1]] + d.cc <- d[cc,] + + } else { + d <- hdesign(g, m, model=model) + d.cc <- d[cc,] + } + + nrow.d <- length(d.cc[,1]) + ncol.d <- length(d.cc[1,]) + d.best <- matrix(rep(0, nrow.d*ncol.d), nrow.d, ncol.d) + + best.cc <- c() + for ( i in 1:nrow.d ) { + + best <- which( d.cc[i,] == max(d.cc[i,]) ) + if ( length(best) == 1 ) { + best.cc[i] = best + } else { + best.cc[i] = sample(best, 1) + } + + d.best[ i, best.cc[i] ] = 1 + } + + hb <- hbrem.true(RX=d.best, HaploidInd=HaploidInd, Ndip=Ndip, Nstrain=Nstrain, Nind=Nind, Npost=Npost, Nbin=Nbin, Ry=Ry) + + reg.full.lm <- lm(Ry ~ d.cc) + reg.null.lm <- lm(Ry ~ 1) + Ftest <- anova(reg.null.lm, reg.full.lm) + pval <- Ftest$"Pr(>F)"[2] + + # cat( m, -log10(pval), "\n" ) + return( c( -log10(pval), hb[[1]], hb[[2]], hb[[3]], hb[[4]], hb[[5]] )) +} + + +hbrem.merge.locus <- function(m, sdp, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) { + + d <- hdesign(g, m, model=model) + d.cc <- d[cc,] + + if ( model == "additive" ) { + + all.0 <- numeric(length(d[,1])) + all.1 <- numeric(length(d[,1])) + for ( i in 1:Nstrain ) { + if ( sdp[i] == 0 ) { + all.0 = ( all.0 + d[,i] ) + } else if ( sdp[i] == 1 ) { + all.1 = ( all.1 + d[,i] ) + } else { + cat("sdp not 0 or 1\n") + } + } + + d.merge <- cbind(all.0, all.1) + + } else if ( model == "full" ) { + + sdp.matrix <- matrix(kronecker(sdp, sdp, paste, sep=""), nrow=Nstrain) + sdp.vector <- c( diag(sdp.matrix), sdp.matrix[upper.tri(sdp.matrix, diag=FALSE)]) + sdp.full <- rep(1,36) + sdp.full[sdp.vector == "00"] = 0 + sdp.full[sdp.vector == "11"] = 2 + + dip.0 <- numeric(d[,1]) + dip.1 <- numeric(d[,1]) + dip.2 <- numeric(d[,1]) + for ( i in 1:Ndip ) { + if ( sdp.full[i] == 0 ) { + dip.0 = ( dip.0 + d[,i] ) + } else if ( sdp.full[i] == 1 ) { + dip.1 = ( dip.1 + d[,i] ) + } else if ( sdp.full[i] == 2 ) { + dip.2 = ( dip.2 + d[,i] ) + } else { + cat("sdp not 0,1 or 2\n") + } + } + + d.merge <- cbind(dip.0,dip.1,dip.2) + + } else { + cat("model not specified: must be one of additive, full\n") + } + + + hb <- hbrem(RX=d.merge, HaploidInd=HaploidInd, Ndip=3, Nstrain=2, Nind=Nind, Npost=Npost, Nbin=Nbin, Ry=Ry) + + reg.full.lm <- lm(Ry ~ d.merge) + reg.null.lm <- lm(Ry ~ 1) + Ftest <- anova(reg.null.lm, reg.full.lm) + pval <- Ftest$"Pr(>F)"[2] + + # cat( m, -log10(pval), "\n" ) + return( c( -log10(pval), hb[[1]], hb[[2]], hb[[3]], hb[[4]], hb[[5]] )) + +} + + +hbrem.region <- function(g, markers, Ry, cc, HaploidInd, Npost, Nbin, Nperm=1000, thres.quick=c(0.5, 0.1, 0.05), thres.precise=c( 0.05, (1/length(markers)) ), thres.method="none", mc.cores=1) { + + if ( HaploidInd == 1 ) { + model = "additive" + Ndip = length(g$strains) + Nstrain = Ndip + Nind = length(Ry) + marker.names <- g$additive$genome$marker[markers] + } + else if ( HaploidInd == 0 ) { + model = "full" + ns = length(g$strains) + Ndip = ns*(ns+1)/2 + Nstrain = ns + Nind = length(Ry) + marker.names <- g$full$genome$marker[markers] + } + + nmark <- length(markers) + + + if ( mc.cores == 1 ) + res = t(sapply ( markers, hbrem.locus, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) ) + else { + res.list=mclapply ( markers, hbrem.locus, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin, mc.cores=mc.cores) + res = t(do.call( "cbind", res.list ) ) + } + + mark.pars.df <- data.frame(res[,1:46]) + names(mark.pars.df) = c( "F.logPval", "Hbar", "sd.Ni", "BIC.qtl", "BIC.null", "BF", "logBF", "DIC.qtl", "DIC.null", "DIC.diff", "pd.qtl", "pd.null", "mode.k", "ga", "gb", "mode.var", "med.k", "med.mu", "med.var", "mean.k", "mean.mu", "mean.var", "hpd.k.50.lower", "hpd.k.50.upper", "hpd.mu.50.lower", "hpd.mu.50.upper", "hpd.var.50.lower", "hpd.var.50.upper", "hpd.k.75.lower", "hpd.k.75.upper", "hpd.mu.75.lower", "hpd.mu.75.upper", "hpd.var.75.lower", "hpd.var.75.upper", "hpd.k.95.lower", "hpd.k.95.upper", "hpd.mu.95.lower", "hpd.mu.95.upper", "hpd.var.95.lower", "hpd.var.95.upper", "hpd.k.99.lower", "hpd.k.99.upper", "hpd.mu.99.lower", "hpd.mu.99.upper", "hpd.var.99.lower", "hpd.var.99.upper") + offset = ncol(mark.pars.df) + mark.pars.df$Name=as.character(marker.names) + if ( class(g) == "happy.genome" ) { + idx = match( marker.names, g[[model]]$genome$marker ) + mark.pars.df$Chr = g[[model]]$genome$chromosome[idx] + mark.pars.df$Bp = g[[model]]$genome$bp[idx] + + bp2 = mark.pars.df$Bp[2:length(mark.pars.df$Bp)] + bidx = which(bp2 < mark.pars.df$Bp[1:(length(mark.pars.df$Bp)-1)])-1 + + mark.pars.df$CumBp = rep(0,nrow(mark.pars.df)) + mark.pars.df$CumBp[bidx+2] = mark.pars.df$Bp[bidx+1] + mark.pars.df$CumBp = cumsum(mark.pars.df$CumBp) + mark.pars.df$Bp + + } + + cat("max logP = ", max(mark.pars.df$F.logPval), "\n") + cat("max mode(k) = ", max(mark.pars.df$mode.k), "\n") + + + hap.means.df <- data.frame( res[,(offset+1):(offset+Ndip)]) + hap.sdevs.df <- data.frame( res[,(offset+Ndip+1):(offset+2*Ndip)]) + hap.avNis.df <- data.frame( res[,(offset+2*Ndip+1):(offset+3*Ndip)]) + strain.means.df <- data.frame( res[,(offset+3*Ndip+1):(offset+3*Ndip+Nstrain)]) + + if ( HaploidInd == 1 ) { + names(hap.means.df) <- g$strains + names(hap.sdevs.df) <- g$strains + names(hap.avNis.df) <- g$strains + names(strain.means.df) <- g$strains + } + else if ( HaploidInd == 0 ) { + strain.names = g$strains + num.strains = length(g$strains) + diplotype.names <- matrix(kronecker(strain.names, strain.names, paste, sep="."), nrow=num.strains) + names.full.symmetric <- c( diag(diplotype.names), diplotype.names[upper.tri(diplotype.names, diag=FALSE)]) + names(hap.means.df) <- names.full.symmetric + names(hap.sdevs.df) <- names.full.symmetric + names(hap.avNis.df) <- names.full.symmetric + names(strain.means.df) <- g$strains + } + + + permuted.y=NULL + if ( is.numeric(Nperm) & Nperm>0 ) { + permuted.y = replicate( Nperm, sample(Ry) ) + } + FlogP.thres=NULL + modek.thres=NULL + + if ( thres.method == "precise" ) { + + FlogP.region = matrix( rep(0, nmark*Nperm), Nperm, nmark ) + modek.region = matrix( rep(0, nmark*Nperm), Nperm, nmark ) + + for ( i in 1:Nperm ) { + + cat("perm ", i, "\n") + + if ( mc.cores == 1 ) + res = t(sapply ( markers, hbrem.locus, g, model, permuted.y[,i], cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) ) + else { + res.list=mclapply ( markers, hbrem.locus, g, model, permuted.y[,i], cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin, mc.cores=mc.cores) + res = t(do.call( "cbind", res.list ) ) + } + + FlogP.vector <- res[,1] + modek.vector <- res[,13] + for ( j in 1:nmark ) { + FlogP.region[i,j] = FlogP.vector[j] + modek.region[i,j] = modek.vector[j] + } + + } + + for ( j in 1:nmark ) { + FlogP.region[,j] <- sort.list(-FlogP.region[,j]) + modek.region[,j] <- sort.list(-modek.region[,j]) + } + + FlogP.thres <- FlogP.region[Nperm*thres.precise,] + modek.thres <- modek.region[Nperm*thres.precise,] + + } else if ( thres.method == "quick" ) { + + FlogP.regionwide.distn <- numeric(Nperm) + modek.regionwide.distn <- numeric(Nperm) + + for ( i in 1:Nperm ) { + + cat("perm ", i, "\n") + + if ( mc.cores == 1 ) + res = t(sapply ( markers, hbrem.perm.locus, g, model, permuted.y[,i], cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) ) + else { + res.list=mclapply ( markers, hbrem.perm.locus, g, model, permuted.y[,i], cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin, mc.cores=mc.cores) + res = t(do.call( "cbind", res.list ) ) + } + + FlogP.vector <- res[,1] + modek.vector <- res[,12] + max.FlogP <- max(FlogP.vector) + max.modek <- max(modek.vector) + + FlogP.regionwide.distn[i] = max.FlogP + modek.regionwide.distn[i] = max.modek + } + + FlogP.sort <- FlogP.regionwide.distn[ sort.list(-FlogP.regionwide.distn) ] + modek.sort <- modek.regionwide.distn[ sort.list(-modek.regionwide.distn) ] + + FlogP.thres <- FlogP.sort[Nperm*thres.quick] + modek.thres <- modek.sort[Nperm*thres.quick] + } + + hbrem.region.list <- list(Summary.Parameters=mark.pars.df, Diplo.Means=hap.means.df, Diplo.StDevs=hap.sdevs.df, Diplo.ExpCounts=hap.avNis.df, Strain.Means=strain.means.df, F.logP.thres=FlogP.thres, Mode.k.thres=modek.thres) + + return(hbrem.region.list) +} diff -Nru r-other-mott-happy-2.1/R/zzz.R r-other-mott-happy-2.4/R/zzz.R --- r-other-mott-happy-2.1/R/zzz.R 2005-02-21 12:47:08.000000000 +0000 +++ r-other-mott-happy-2.4/R/zzz.R 2013-12-19 14:52:58.000000000 +0000 @@ -1,5 +1,2 @@ -.packageName <- "happy_2.0.1" +.packageName <- "happy.hbrem_2.4" -.First.lib <- function(lib, pkg) { - library.dynam('happy', pkg, lib) -} diff -Nru r-other-mott-happy-2.1/debian/changelog r-other-mott-happy-2.4/debian/changelog --- r-other-mott-happy-2.1/debian/changelog 2013-05-20 10:16:24.000000000 +0000 +++ r-other-mott-happy-2.4/debian/changelog 2014-01-24 16:45:31.000000000 +0000 @@ -1,8 +1,31 @@ -r-other-mott-happy (2.1-7build1) saucy; urgency=low +r-other-mott-happy (2.4-1) unstable; urgency=low - * Rebuild against R 3.0.0. + * New upstream version (Closes: #709190). + * Revived the debug package. - -- Colin Watson Mon, 20 May 2013 11:16:24 +0100 + -- Steffen Moeller Fri, 24 Jan 2014 15:48:04 +0100 + +r-other-mott-happy (2.3-1) unstable; urgency=low + + * Team upload. + + [ Charles Plessy ] + * New upstream version, renaming happy to happy.hbrem. + * Simplified debian/rules following changes in r-base-dev 2.14.2~20120222-1 + * Echo a dummy NAMESPACE file (Closes: #709190). + * Use Debhelper 9. + * Refreshed quilt patches. + * Added Debhelper's misc:Depends substitution variable. + * Normalised VCS URLs. + * Normalised debian/control with config-model-edit. + * Converted Debian copyright file to machine-readable format. + * Conforms to Debian Policy version 3.9.4. + * renamed debian/upstream-metadata.yaml to debian/upstream. + + [ Andreas Tille ] + * debian/upstream: Fixed spelling of journal field + + -- Charles Plessy Sun, 30 Jun 2013 11:25:27 +0900 r-other-mott-happy (2.1-7) unstable; urgency=low diff -Nru r-other-mott-happy-2.1/debian/compat r-other-mott-happy-2.4/debian/compat --- r-other-mott-happy-2.1/debian/compat 2011-12-17 12:15:37.000000000 +0000 +++ r-other-mott-happy-2.4/debian/compat 2013-08-31 15:52:00.000000000 +0000 @@ -1 +1 @@ -8 +9 diff -Nru r-other-mott-happy-2.1/debian/control r-other-mott-happy-2.4/debian/control --- r-other-mott-happy-2.1/debian/control 2011-12-17 19:03:26.000000000 +0000 +++ r-other-mott-happy-2.4/debian/control 2014-01-24 16:38:56.000000000 +0000 @@ -1,19 +1,32 @@ Source: r-other-mott-happy -Section: gnu-r -Priority: optional Maintainer: Debian Med Packaging Team Uploaders: Steffen Moeller , - Andreas Tille -Build-Depends: debhelper (>= 8), cdbs, r-base-dev (>= 2.12), r-cran-g.data, r-cran-multicore, r-cran-mass -Standards-Version: 3.9.2 + Andreas Tille +Section: gnu-r +Priority: optional +Build-Depends: debhelper (>= 9), + cdbs, + r-base-dev (>= 2.12), + r-cran-g.data, + r-cran-multicore, + r-cran-mass +Standards-Version: 3.9.5 +Vcs-Browser: http://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-other-mott-happy/trunk/ +Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-other-mott-happy/trunk/ Homepage: http://www.well.ox.ac.uk/happy/happyR.shtml -Vcs-Browser: http://svn.debian.org/wsvn/debian-med/trunk/packages/R/r-other-mott-happy/trunk/ -Vcs-Svn: svn://svn.debian.org/debian-med/trunk/packages/R/r-other-mott-happy/trunk/ -Package: r-other-mott-happy +Package: r-other-mott-happy.hbrem Architecture: any -Depends: ${shlibs:Depends}, ${R:Depends}, r-cran-g.data, r-cran-mass -Recommends: r-cran-vr, r-cran-multicore +Depends: ${shlibs:Depends}, + ${R:Depends}, + ${misc:Depends}, + r-cran-g.data, + r-cran-mass +Recommends: r-cran-vr, + r-cran-multicore +Provides: r-other-mott-happy +Conflicts: r-other-mott-happy, r-other-mott-happy.hbrem-dbg (<< ${binary:Version}) +Replaces: r-other-mott-happy Description: GNU R package for fine-mapping complex diseases Happy is an R interface into the HAPPY C package for fine-mapping Quantitative Trait Loci (QTL) in Heterogenous Stocks (HS). An HS is @@ -30,3 +43,16 @@ . Read /usr/share/doc/r-other-mott-happy/README.Debian for a more detailed explanation. + +Package: r-other-mott-happy.hbrem-dbg +Architecture: any +Section: debug +Depends: ${misc:Depends}, r-other-mott-happy.hbrem (= ${binary:Version}) +Description: Debug information for Happy R package native library + The documentation of Happy is fairly complete, but the runtime itself + does not inform as verbosely as users would like it on errors in + configuration files. This debug package may complement nicely the + on or other missing check on data integrity. + . + With this package installed, a debugger like 'gdb' indicates the + position in the C implementation of R routines that spawned an error. diff -Nru r-other-mott-happy-2.1/debian/copyright r-other-mott-happy-2.4/debian/copyright --- r-other-mott-happy-2.1/debian/copyright 2009-09-03 10:49:37.000000000 +0000 +++ r-other-mott-happy-2.4/debian/copyright 2013-08-31 15:52:00.000000000 +0000 @@ -1,16 +1,21 @@ -The R library combinat was written by Richard Mott - and is also maintained by him. +Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/ +Source: http://www.well.ox.ac.uk/happy/happy.hbrem_2.3.tar.gz -This package was created by Steffen Moeller . - -The sources were downloaded from http://www.well.ox.ac.uk/happy/happyR.shtml. - -The package was renamed from its upstream name 'happy' to -'r-other-mott-happy' in harmony with the R packaging policy to indicate -that the package is external to the CRAN or BioC repositories. - -happy is copyright 1995-2008 of Richard Mott, and released under the -GNU General Public License (GPL) version 2 or later versions. - -On a Debian GNU/Linux system, the GPL license version 2 is included in the file -/usr/share/common-licenses/GPL-2. +Files: * +Copyright: © Richard Mott +License: GPL-2+ + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + . + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + . + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +Comment: On a Debian GNU/Linux system, the GPL license version 2 is included + in the file /usr/share/common-licenses/GPL-2. diff -Nru r-other-mott-happy-2.1/debian/dirs r-other-mott-happy-2.4/debian/dirs --- r-other-mott-happy-2.1/debian/dirs 2011-12-17 13:21:46.000000000 +0000 +++ r-other-mott-happy-2.4/debian/dirs 1970-01-01 00:00:00.000000000 +0000 @@ -1 +0,0 @@ -usr/lib/R/site-library diff -Nru r-other-mott-happy-2.1/debian/patches/ELT_String.patch r-other-mott-happy-2.4/debian/patches/ELT_String.patch --- r-other-mott-happy-2.1/debian/patches/ELT_String.patch 2011-12-17 12:37:10.000000000 +0000 +++ r-other-mott-happy-2.4/debian/patches/ELT_String.patch 1970-01-01 00:00:00.000000000 +0000 @@ -1,32 +0,0 @@ ---- happy.orig/src/Rhappy.c -+++ happy/src/Rhappy.c -@@ -195,21 +195,21 @@ - - - SET_VECTOR_ELT(ans,0,strains); -- SET_VECTOR_ELT(names,0,mkChar("strains")); -+ SET_STRING_ELT(names,0,mkChar("strains")); - SET_VECTOR_ELT(ans,1,markers); -- SET_VECTOR_ELT(names,1,mkChar("markers")); -+ SET_STRING_ELT(names,1,mkChar("markers")); - SET_VECTOR_ELT(ans,2,map); -- SET_VECTOR_ELT(names,2,mkChar("map")); -+ SET_STRING_ELT(names,2,mkChar("map")); - SET_VECTOR_ELT(ans,3,subjects); -- SET_VECTOR_ELT(names,3,mkChar("subjects")); -+ SET_STRING_ELT(names,3,mkChar("subjects")); - SET_VECTOR_ELT(ans,4,phenotypes); -- SET_VECTOR_ELT(names,4,mkChar("phenotypes")); -+ SET_STRING_ELT(names,4,mkChar("phenotypes")); - SET_VECTOR_ELT(ans,5,handle); -- SET_VECTOR_ELT(names,5,mkChar("handle")); -+ SET_STRING_ELT(names,5,mkChar("handle")); - SET_VECTOR_ELT(ans,6,chromosome); -- SET_VECTOR_ELT(names,6,mkChar("chromosome")); -+ SET_STRING_ELT(names,6,mkChar("chromosome")); - SET_VECTOR_ELT(ans,7,family); -- SET_VECTOR_ELT(names,7,mkChar("family")); -+ SET_STRING_ELT(names,7,mkChar("family")); - - setAttrib( ans, R_NamesSymbol, names ); - diff -Nru r-other-mott-happy-2.1/debian/patches/happyMatrices.patch r-other-mott-happy-2.4/debian/patches/happyMatrices.patch --- r-other-mott-happy-2.1/debian/patches/happyMatrices.patch 2011-12-17 12:37:05.000000000 +0000 +++ r-other-mott-happy-2.4/debian/patches/happyMatrices.patch 2013-08-31 15:52:00.000000000 +0000 @@ -1,6 +1,8 @@ ---- happy.orig/R/happy.R -+++ happy/R/happy.R -@@ -72,10 +72,12 @@ +Index: r-other-mott-happy-2.3/R/happy.R +=================================================================== +--- r-other-mott-happy-2.3.orig/R/happy.R ++++ r-other-mott-happy-2.3/R/happy.R +@@ -91,10 +91,12 @@ for( m in h$markers) { add <- hdesign( h, m, model='additive' ) full <- hdesign( h, m, model='full' ) diff -Nru r-other-mott-happy-2.1/debian/patches/series r-other-mott-happy-2.4/debian/patches/series --- r-other-mott-happy-2.1/debian/patches/series 2011-12-17 11:58:58.000000000 +0000 +++ r-other-mott-happy-2.4/debian/patches/series 2014-01-24 14:56:11.000000000 +0000 @@ -1,2 +1 @@ happyMatrices.patch -ELT_String.patch diff -Nru r-other-mott-happy-2.1/debian/r-cran_without_simple-patchsys.mk r-other-mott-happy-2.4/debian/r-cran_without_simple-patchsys.mk --- r-other-mott-happy-2.1/debian/r-cran_without_simple-patchsys.mk 2011-12-17 17:00:48.000000000 +0000 +++ r-other-mott-happy-2.4/debian/r-cran_without_simple-patchsys.mk 1970-01-01 00:00:00.000000000 +0000 @@ -1,108 +0,0 @@ -#!/usr/bin/make -f -# -*- makefile -*- -# -# Generic debian/rules file for the Debian/GNU Linux r-cran-* packages -# -# Should be sufficient for Debianization of CRAN (http://cran.r-project.org) -# packages. Note that you still need to provide the other files in debian/*, -# in particular control, changelog and copyright. -# -# Copyright 2003-2008 by Dirk Eddelbuettel - -include /usr/share/cdbs/1/rules/debhelper.mk -include /usr/share/cdbs/1/class/langcore.mk -## include /usr/share/cdbs/1/rules/dpatch.mk -# Check whether source format 3.0 (quilt) is used. If yes, do not include the conflicting simple-patchsys.mk -formatfile := $(CURDIR)/debian/source/format -format_3_quilt = $(shell if [ -f $(formatfile) ] ; then if grep -q '3.0[[:space:]]*(quilt)' $(formatfile) ; then echo 1 ; else echo 0 ; fi else echo 0 ; fi ) -ifeq ($(format_3_quilt),0) - include /usr/share/cdbs/1/rules/simple-patchsys.mk -endif - -# awk command to extract word after Package or Bundle, not lowercased -awkString := "'/^(Package|Bundle):/ {print $$2 }'" - -# apply it to the upstream meta-info file DESCRIPTION, also generate a lc version -cranNameOrig := $(shell awk "$(awkString)" DESCRIPTION) -cranName := $(shell echo "$(cranNameOrig)" | tr A-Z a-z) - -## if no debRreposname is known, set default to cran -- thanks, Steffen! -ifeq ($(debRreposname),) - debRreposname := cran -endif - -## we can define additional flags for R's make, eg "CXXFLAGS=-g0" for -## RQuantLib but the default is empty -## makeFlags := -## if makeFlags are defined, then we'll use them in this variable -## which would otherwise be empty -ifneq ($(makeFlags),) - makeFlagsCall := MAKEFLAGS=$(makeFlags) -endif - -## and use the results to build the Debian'ized package name -package := r-$(debRreposname)-$(cranName) - -## awk command to extract word after Priority -prioritystr := "'/^Priority:/ {print tolower($$2) }'" -priority := $(shell awk "$(prioritystr)" DESCRIPTION) - -ifeq ($(priority),recommended) - debRdir := usr/lib/R/library -else - debRdir := usr/lib/R/site-library -endif - -## we use these results for the to-be-installed-in directory -debRlib := $(CURDIR)/debian/$(package)/$(debRdir) - -## optional installation of a lintian silencer -lintiandir := $(CURDIR)/debian/$(package)/usr/share/lintian/overrides - -common-install-indep:: R_any_arch -common-install-arch:: R_any_arch - -R_any_arch: - ## create the target directory - dh_installdirs $(debRdir) - ## - ## call R to install the sources we're looking at - ## use this inside xvfb-run if this wrapper is installed - if test -f /usr/bin/xvfb-run; then \ - $(makeFlagsCall) xvfb-run -a \ - R CMD INSTALL -l $(debRlib) --clean \ - $(extraInstallFlags) . ; \ - else \ - $(makeFlagsCall) R CMD INSTALL -l $(debRlib) \ - --clean $(extraInstallFlags) . ;\ - fi - ## remove extra files which are present in some packages - rm -vf $(debRlib)/R.css \ - $(debRlib)/$(cranNameOrig)/COPYING \ - $(debRlib)/$(cranNameOrig)/LICENSE.txt - ## if we have an overrides file for lintian, install it - if test -f debian/overrides; then \ - install -d $(lintiandir) ; \ - install -m 0644 debian/overrides \ - $(lintiandir)/$(package); \ - fi - -## clean target from patch by Steffen Moeller on 16 May 2009 -clean:: - ## the re-invocation of a build process should not - ## leave a footprint in Debian's diff.gz. - if test -d src; then \ - find src -regex ".*\..*o" | \ - xargs --no-run-if-empty -r rm; \ - fi - rm -f config.log config.status - ## the configure file is provided by upstream but - ## could be recreated by a call to 'autoconf'. - #if [ -r configure.in ]; then \ - # rm -f configure \ - #fi - ## - # if [ -r src/Makevars.in ]; then \ - # rm -f src/Makevars; \ - # fi - diff -Nru r-other-mott-happy-2.1/debian/rules r-other-mott-happy-2.4/debian/rules --- r-other-mott-happy-2.1/debian/rules 2011-12-17 17:17:08.000000000 +0000 +++ r-other-mott-happy-2.4/debian/rules 2014-01-10 13:39:07.000000000 +0000 @@ -1,14 +1,18 @@ #!/usr/bin/make -f -# -*- makefile -*- -# debian/rules file for the Debian/GNU Linux r-other-mott-happy package -# Copyright 2008 by Andreas Tille debRreposname := other-mott -## FIXME: Dirty workaround for bug #652456 -# include /usr/share/R/debian/r-cran.mk -include debian/r-cran_without_simple-patchsys.mk - -# Require a number equal or superior than the R version the package was built with. -install/$(package):: - echo "R:Depends=r-base-core (>= $(shell R --version | head -n1 | perl -ne 'print / +([0-9]\.[0-9]+\.[0-9])/')~)" >> debian/$(package).substvars +include /usr/share/R/debian/r-cran.mk + +# Remove once upstream adds a NAMESPACE file + +pre-build:: + [ -f NAMESPACE ] || echo 'exportPattern("^[[:alpha:]]+")' > NAMESPACE + +override_dh_strip: + dh_strip -Pr-other-mott-happy.hbrem --dbg-package=r-other-mott-happy.hbrem-dbg + + +override_dh_clean: + $(RM) NAMESPACE + dh_clean diff -Nru r-other-mott-happy-2.1/debian/upstream r-other-mott-happy-2.4/debian/upstream --- r-other-mott-happy-2.1/debian/upstream 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/debian/upstream 2013-08-31 15:52:00.000000000 +0000 @@ -0,0 +1,15 @@ +Contact: Richard Mott +Homepage: http://www.well.ox.ac.uk/happy/happyR +Name: happy.hbrem +Reference: + author: Mott, Richard and Talbot, Christopher J. and Turri, Maria G. and Collins, Allan C. and Flint, Jonathan + title: A method for fine mapping quantitative trait loci in outbred animal stocks + journal: Proc. Natl. Acad. Sci. USA + volume: 97 + number: 23 + pages: 12649-12654 + year: 2000 + PMID: 11050180 + DOI: 10.1073/pnas.230304397 + URL: http://www.pnas.org/content/97/23/12649.abstract + eprint: http://www.pnas.org/content/97/23/12649.full.pdf+html diff -Nru r-other-mott-happy-2.1/debian/upstream-metadata.yaml r-other-mott-happy-2.4/debian/upstream-metadata.yaml --- r-other-mott-happy-2.1/debian/upstream-metadata.yaml 2011-12-17 11:58:58.000000000 +0000 +++ r-other-mott-happy-2.4/debian/upstream-metadata.yaml 1970-01-01 00:00:00.000000000 +0000 @@ -1,15 +0,0 @@ -Contact: Richard Mott -DOI: 10.1073/pnas.230304397 -Homepage: http://www.well.ox.ac.uk/happy/happyR -Name: HAPPY -PMID: 11050180 -Reference: - author: Mott, Richard and Talbot, Christopher J. and Turri, Maria G. and Collins, Allan C. and Flint, Jonathan - title: A method for fine mapping quantitative trait loci in outbred animal stocks - journal: Proc Natl Acad Sci U S A - volume: 97 - number: 23 - pages: 12649-12654 - year: 2000 - URL: http://www.pnas.org/content/97/23/12649.abstract - eprint: http://www.pnas.org/content/97/23/12649.full.pdf+html diff -Nru r-other-mott-happy-2.1/man/AAA-happy.Rd r-other-mott-happy-2.4/man/AAA-happy.Rd --- r-other-mott-happy-2.1/man/AAA-happy.Rd 2008-02-01 17:02:49.000000000 +0000 +++ r-other-mott-happy-2.4/man/AAA-happy.Rd 2013-12-19 14:26:44.000000000 +0000 @@ -207,16 +207,15 @@ } -} + \references{ Mott R, Talbot CJ, Turri MG, Collins AC, Flint J. A method for fine mapping quantitative trait loci in outbred animal stocks. Proc Natl Acad Sci U S A. 2000 Nov 7;97(23):12649-54. } \usage{ -happy( datafile, allelesfile, generations=200, standardise=FALSE, - phase="unknown", file.format="happy", missing.code="NA", -do.dp=TRUE, min.dist=1.0e-5, mapfile=NULL, ancestryfile=NULL, haploid=FALSE ) +happy( datafile, allelesfile, generations=200, phase="unknown", + file.format="happy", missing.code="NA", do.dp=TRUE, min.dist=1.0e-5, mapfile=NULL, ancestryfile=NULL, haploid=FALSE ) happy.matrices( h ) happy.save( h, file ) @@ -225,13 +224,9 @@ \arguments{ \item{datafile}{ name of the text file containing the genotype and phenotype data in HAPPY format} - \item{allelesfile} { name of the text file containing the + \item{allelesfile}{ name of the text file containing the allele/strain data in HAPPY format} \item{generations}{ the number of breeding generations in the HS} - \item{standardise}{ if TRUE then the phenotype is transformed to have - mean 0 and variance 1. This does not affect the identification of QTL - but can make the interpretation of effects easier.} - \item{phase}{ If phase=="unknown" then the phase of the genotypes is unknown and no attempt is made to infer it. If phase="estimate" then it is estimated using parental genotype data @@ -306,11 +301,10 @@ phenotypes} \item{handle}{a numeric index used internally by the C-code. Do not change.} - \item{matrices}{ a list of matrices used in model fitting (only - created after a call to happy.matrices()) + \item{matrices}{ a list of matrices used in model fitting (only created after a call to happy.matrices()).} \item{use.pedigrees}{ boolean variable indicating whether pedigree information was used to help determine the phase of the genotypes} \item{phase.known}{ boolean variable indicating whether or not the phase of the genotypes is assumed to be known} -} + happy.save() will save a happy object to a file so that it can be re-used in a later session with the load() command. @@ -328,8 +322,8 @@ performed once, the data can be persisted to disk and then re-used (e.g. to analyse multiple phenotypes) at a later date. - } - + +} \seealso{ hfit(), mergefit(), happyplot() } \examples{ \dontrun{h <- happy('HS.data', 'HS.alleles', generations=200)} diff -Nru r-other-mott-happy-2.1/man/gauss.Rd r-other-mott-happy-2.4/man/gauss.Rd --- r-other-mott-happy-2.1/man/gauss.Rd 2005-02-25 13:57:19.000000000 +0000 +++ r-other-mott-happy-2.4/man/gauss.Rd 2013-12-19 14:26:44.000000000 +0000 @@ -37,6 +37,8 @@ The model-fitting is implemented in the function gfit() by an iterative process, rather like a simplified version of EM. Is is slower than hfit(), and generally gives similar results as far as overall QTL detection is concered,m but gives more accurate parameter estimates. The log-likelihood for the data is + + \deqn{ L = \sum_{i} \log ( \sum_j p_{ij} \frac{\exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})}{\sqrt{2\pi \sigma^2}}) } @@ -44,29 +46,30 @@ Differentiating wrt to the parameters gives -\deqn{ \frac{\partial L}{\partial \sigma^2} = \sum_i \frac{\sum_j p_{ij} (y_i-\beta_j)^2 \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2}) }{ 2\sigma^4 \sum_j p_{ij} \exp(-\frac{(y_i-\beta_j)^2 }{2\sigma^2})} - \frac{N}{2\sigma^2} } -\deqn{ \frac{\partial L}{\partial \beta_j } = - \sum_i \frac{ p_{ij} \frac{(y_i-\beta_j) }{\sigma^2} \exp( -\frac{(y_i-\beta_j)^2}{2\sigma^2})}{ \sum_j e_{ij}} } -\deqn{ = \frac{1}{\sigma^2} \left( - \sum_i \frac{y_i e_{ij} }{\sum_j e_{ij}} + \beta_j \frac{\sum_i e_{ij} }{\sum_j e_{ij}} \right) } - -write -\deqn{ w_{ij} = \frac{p_{ij} \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2}) }{ \sum_j \frac{p_{ij} \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})}}} - -then the mle satisfies +XXXX -\deqn{ \hat{\sigma^2} = \frac{1}{N} \sum_i \frac{\sum_j p_{ij}(y_i-\beta_j)^2 \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})} {\sum_j p_{ij}\exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})}} +\deqn{ \frac{\partial L}{\partial \sigma^2} = \sum_i \frac{\sum_j p_{ij} (y_i-\beta_j)^2 \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2}) }{ 2\sigma^4 \sum_j p_{ij} \exp(-\frac{(y_i-\beta_j)^2 }{2\sigma^2})} - \frac{N}{2\sigma^2} } +\deqn{ \frac{\partial L}{\partial \beta_j } = - \sum_i \frac{ p_{ij} \frac{(y_i-\beta_j) }{\sigma^2} \exp( -\frac{(y_i-\beta_j)^2}{2\sigma^2})}{ \sum_j e_{ij}} } +\deqn{ = \frac{1}{\sigma^2} \left( - \sum_i \frac{y_i e_{ij} }{\sum_j e_{ij}} + \beta_j \frac{\sum_i e_{ij} }{\sum_j e_{ij}} \right) } \deqn{ \hat{\sigma^2} = \frac{1}{N} \sum_i \sum_j \hat{w}_{ij}(y_i-\hat{\beta}_j)^2 } +write +\deqn{ w_{ij} = \frac{p_{ij} \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2}) }{ \sum_j p_{ij} \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})}} + +then the mle satisfies \deqn{ \hat{\beta_j} = \frac{\sum_i \hat{e}_{ij} y_i}{\sum_i \hat{e}_{ij}} } +\deqn{ \hat{\sigma^2} = \frac{1}{N} \sum_i \frac{\sum_j p_{ij}(y_i-\beta_j)^2 \exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})} {\sum_j p_{ij}\exp(-\frac{(y_i-\beta_j)^2}{2\sigma^2})}} + and the log-likelihood is \deqn{ \hat{L} = \sum_i(\log \sum_j \hat{e}_{ij}) - \frac{N \log(2\pi\hat{\sigma}^2)}{2} } @@ -93,7 +96,7 @@ \item{eps}{ the terminatation accuracy in the model fitting : the log likelihood must change by less than eps in successive iterations} \item{df}{ the degress of freedom to use. If NULL then this is computed as the rank of the data} \item{n}{ the number of observations with non-missing phenotypes} - \item{sigma2} {the variance of the errors in the data} + \item{sigma2}{the variance of the errors in the data} \item{params}{ a list with two components, beta = the group means and sigma = the standard deviation} \item{p}{ vector of paramters. For internal use only} } diff -Nru r-other-mott-happy-2.1/man/happy-internal.Rd r-other-mott-happy-2.4/man/happy-internal.Rd --- r-other-mott-happy-2.1/man/happy-internal.Rd 2008-02-01 16:53:05.000000000 +0000 +++ r-other-mott-happy-2.4/man/happy-internal.Rd 2013-12-19 15:11:57.000000000 +0000 @@ -14,7 +14,15 @@ \alias{hprob2} \alias{hnonrecomb} \alias{h.sum.prob2} - +\alias{delete.happy.cobject} +\alias{happy.load.data} +\alias{hbrem} +\alias{hbrem.locus} +\alias{hbrem.region} +\alias{hbrem.true} +\alias{hbrem.merge.locus} +\alias{hbrem.perm.locus} +\alias{save.happy.internal} \usage{ matrixSquared( matrix1, matrix2 ) @@ -27,6 +35,9 @@ strain.effects( h, fit, family='gaussian' ) glmfit( formula=NA, family='gaussian' ) sdp( strains, alleles ) +hbrem.perm.locus(m, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) +hbrem.true( RX, HaploidInd, Ndip, Nstrain, Nind, Npost=2000, Nbin, Ry ) +hbrem.merge.locus(m, sdp, g, model, Ry, cc, HaploidInd, Ndip, Nstrain, Nind, Npost, Nbin) } \author{Richard Mott} \keyword{internal} diff -Nru r-other-mott-happy-2.1/man/hcache.Rd r-other-mott-happy-2.4/man/hcache.Rd --- r-other-mott-happy-2.1/man/hcache.Rd 2008-02-01 16:58:48.000000000 +0000 +++ r-other-mott-happy-2.4/man/hcache.Rd 2013-12-19 16:10:14.000000000 +0000 @@ -18,12 +18,13 @@ into memory. \code{save.happy()} saves a single happy object as a delayed data package. \code{the.chromosomes()} is a conveniemce funtion that generates a character vector of chromosome names.} -} + \usage{ save.genome( gdir, sdir, prefix, chrs=NULL, file.format="ped", -mapfile=NULL,ancestryfile=NULL, generations=50, phase="unknown", haploid=FALSE ) -genome <- load.genome( sdir, use.X=TRUE, chr=the.chromosomes(use.X=use.X) ) -marker.list <- load.markers( genome, markers ) +mapfile=NULL,ancestryfile=NULL, generations=50, phase="unknown", haploid=FALSE, mc.cores=1 ) +the.chromosomes(autosomes=19, use.X=FALSE) +load.genome( sdir, use.X=TRUE, chrs = the.chromosomes(use.X=use.X), n.chr=NA, models=c("additive", "full", "genotype")) +load.markers( genome, markers, model="additive", include.models=FALSE ) save.happy( h, pkg, dir, model="additive" ) } \arguments{ @@ -37,9 +38,14 @@ \code{save.genome()}. An attempt is made to find files in \code{gdir} named like \code{chrN.prefix.*} where N is the chromosome number (1...20, X, Y), as defined in \code{chrs}.} - \item{chrs}{ List of chromosome numbers to be processed.} - \item{use.X}{ logical to determine whether to use X-chromsome data, in + \item{chrs}{ number of autosomes.} + \item{n.chr}{ Alternative way of specifying the number of + chromosomes. Must be an integer or NA.} + \item{autosomes}{Sequence of autosomal chromosome identifiers} + \item{use.X}{Logical to determine whether to use X-chromsome data, in load.genome().} + \item{models}{ list of strings specifies the types of data to load. "additive" and "full" are design matrices corresponding to the additive and full models, "genotypes" are the raw genotypes.} +\item{include.models}{ Boolean indicating the type of return object - if TRUE then an additional list specifiying if the model is additive, full or genotype is returned.} \item{file.format}{Defines the input genotype file format, either "ped" (Ped file format) or "happy" ( HAPPY .data file format).} \item{mapfile}{ Name of a text file containing the physicla (base @@ -79,6 +85,7 @@ \item{dir}{ Name of directory to create a delayed data package for a single happy object} \item{model}{ One of "additive", "full", "genotype"}. + \item{mc.cores}{Split computation across this number of cores} } \value{ \code{save.genome()} returns NULL. diff -Nru r-other-mott-happy-2.1/man/hdesign.Rd r-other-mott-happy-2.4/man/hdesign.Rd --- r-other-mott-happy-2.1/man/hdesign.Rd 2005-09-20 13:52:39.000000000 +0000 +++ r-other-mott-happy-2.4/man/hdesign.Rd 2013-12-19 14:26:44.000000000 +0000 @@ -29,8 +29,8 @@ \value{ hdesign() returns a design matrix \eqn{d_{ij}}{d[i,j]}, in which the \eqn{i}th row corresponds to the subject \eqn{i}, and the \eqn{j}th column to the corresponding strain or combination of strains or merged strains. - \hprob() returns a matrix \eqn(p_{ix}{p[i,x]}, in which the \eqn{i} th row corresponds to the subject \eqn{i}, and the \eqn{x=s*S+t} th column contains the probability that the ancestral strains are \eqn{s,t} where \eqn{S} is the total number of strains. - \hgenotype() returns a \eqn{Nx2} matrix \eqn(g_{ix}{g[i,x]} in which + hprob() returns a matrix \eqn{p_{ix}{p[i,x]}}, in which the \eqn{i} the row corresponds to the subject \eqn{i}, and the \eqn{x=s*S+t} th column contains the probability that the ancestral strains are \eqn{s,t} where \eqn{S} is the total number of strains. + hgenotype() returns a \eqn{Nx2} matrix \eqn{g_{ix}{g[i,x]}} in which the \eqn{i} th row corresponds to the subject \eqn{i}, and column 1 contains the first allele and column 2 the second allele at the marker specified, or (if \code{collapse=TRUE}) a vector of genotypes with the diff -Nru r-other-mott-happy-2.1/man/hfit.Rd r-other-mott-happy-2.4/man/hfit.Rd --- r-other-mott-happy-2.1/man/hfit.Rd 2007-09-10 09:10:45.000000000 +0000 +++ r-other-mott-happy-2.4/man/hfit.Rd 2013-12-19 14:26:44.000000000 +0000 @@ -110,6 +110,7 @@ \item{permdata}{ a list containing the results of the permutation analysis, or NULL if permute=0. The list contains the following elements: + \itemize{ \item{N}{ The number of permutations } \item{permutation.dist}{ A vector containing sorted ANOVA logP values from the N permutations. These values can be used to estimate the shape of @@ -127,6 +128,7 @@ permutation logP at the marker interval that exceed the logP for that interval. } + } } diff -Nru r-other-mott-happy-2.1/src/Rhappy.c r-other-mott-happy-2.4/src/Rhappy.c --- r-other-mott-happy-2.1/src/Rhappy.c 2008-02-01 16:33:57.000000000 +0000 +++ r-other-mott-happy-2.4/src/Rhappy.c 2013-12-19 16:21:35.000000000 +0000 @@ -11,27 +11,28 @@ #include #include"happy.h" -#include"cl.h" +#include"readline.h" #include"cmp.h" #include"stats.h" static QTL_DATA *qtldata[100]; static int nqtldata = 0; +int entrycmp( const void *a, const void *b ); SEXP happy( SEXP datafile, SEXP allelesfile, SEXP generations, SEXP phase, SEXP file_format, SEXP missing_code, SEXP do_dp, SEXP min_dist, SEXP haploid, SEXP ancestryfile ) { QTL_DATA *q = NULL; ALLELES *a = NULL; FILE *dfp=NULL, *afp=NULL, *anfp=NULL; - const char *afilename, *anfilename; - const char *dfilename; + const char *afilename=NULL, *anfilename=NULL; + const char *dfilename=NULL; int gen = 0; double g; int verbose = 0; int use_parents = 0; int ped_format = 0; int phaseKnown = 0; - SEXP strains,markers, chromosome, subjects, phenotypes, map, family; + SEXP strains,markers, chromosome, subjects, phenotypes, map, family, mother, father; SEXP ans = R_NilValue; SEXP names, handle, class; char *MissingCode; @@ -64,7 +65,7 @@ if ( ! isNumeric(generations) || length(generations) != 1 ) error( "generations is not numeric"); g = REAL(generations)[0]; - gen = (int)g; + gen = (int)g; if ( ! isString(phase) || length(phase) != 1 ) error( "phase is not a string"); @@ -79,17 +80,17 @@ if ( strlen( CHAR(STRING_ELT(missing_code,0)) ) > 0 ) { MissingCode = (char*)CHAR(STRING_ELT(missing_code,0)); } else { - strcpy(MissingCode,ND_ALLELE); + MissingCode = strdup(ND_ALLELE); } if ( ! isNumeric(do_dp) || length(do_dp) != 1 ) error( "do_dp is not numeric(1)"); - Do_dp = (int)REAL(do_dp)[0]; + Do_dp = INTEGER(do_dp)[0]; if ( ! isNumeric(haploid) || length(haploid) != 1 ) error( "haploid is not numeric(1)"); - Haploid = (int)REAL(haploid)[0]; - + Haploid = INTEGER(haploid)[0]; + if ( ! isNumeric(min_dist) || length(min_dist) != 1 ) error( "min_dist is not numeric(1)"); else if ( isNumeric(min_dist) ) @@ -116,11 +117,9 @@ phaseKnown = 1; } - if ( use_parents ) Rprintf( "using parental genotypes to help determine phase\n"); a = input_allele_frequencies( afp, gen, MissingCode, MinDist, verbose ); - Rprintf( "a->markers %d\n", a->markers ); q = read_qtl_data( dfp, (char*)dfilename, a, verbose, use_parents, ped_format, MissingCode ); q->an = read_subject_ancestries( anfp, (char*)anfilename, verbose ); q->phase_known = phaseKnown; @@ -129,47 +128,42 @@ if ( q->an ) check_and_apply_ancestry( q ); - Rprintf( "dfile %s afile %s gen %d\n", dfilename, afilename, gen ); - + Rprintf( "datafile %s allelesfile %s gen %d\n", dfilename, afilename, gen ); - if ( Do_dp ) + if ( Do_dp ) { if ( q->haploid ) { /* heterozygosity(q ); */ create_haploid_summed_dp_matrices( q ); } else create_summed_dp_matrices( q ); + } PROTECT(strains=allocVector(STRSXP,q->S)); for(i=0;iS;i++) { SET_STRING_ELT(strains,i, mkChar(a->strain_name[i])); } - PROTECT(markers=allocVector(STRSXP,q->M)); for(i=0;iM;i++) { SET_STRING_ELT(markers,i, mkChar(a->af[i].marker_name)); } - PROTECT(chromosome=allocVector(STRSXP,q->M)); for(i=0;iM;i++) { SET_STRING_ELT(chromosome,i, mkChar(a->af[i].chromosome)); } - PROTECT(map=allocVector(REALSXP,q->M)); for(i=0;iM;i++) { REAL(map)[i] = a->af[i].position; } - PROTECT(subjects=allocVector(STRSXP,q->N)); for(i=0;iN;i++) { SET_STRING_ELT(subjects,i, mkChar(q->name[i])); } - PROTECT(family=allocVector(STRSXP,q->N)); for(i=0;iN;i++) { if ( q->family && q->family[i] ) @@ -178,23 +172,44 @@ SET_STRING_ELT(family,i, mkChar("")); } + PROTECT(mother=allocVector(INTSXP,q->N)); + if ( q->mother ) { + for(i=0;iN;i++) { + INTEGER(mother)[i] = q->mother[i]+1; + } + } + else { + for(i=0;iN;i++) { + INTEGER(mother)[i] = -1; + } + } + + PROTECT(father=allocVector(INTSXP,q->N)); + if ( q->father ){ + for(i=0;iN;i++) { + INTEGER(father)[i] = q->father[i]+1; + } + } + else { + for(i=0;iN;i++) { + INTEGER(father)[i] = -1; + } + } PROTECT(phenotypes=allocVector(REALSXP,q->N)); for(i=0;iN;i++) { REAL(phenotypes)[i] = q->observed[i]; } - - PROTECT( ans = allocVector( VECSXP, 8 ) ); - PROTECT(names=allocVector(STRSXP, 8 ) ); + PROTECT( ans = allocVector( VECSXP, 10 ) ); + PROTECT(names=allocVector(STRSXP, 10 ) ); PROTECT(handle = allocVector(INTSXP,1)); qtldata[nqtldata] = q; INTEGER(handle)[0] = nqtldata++; - - SET_VECTOR_ELT(ans,0,strains); + /* SET_VECTOR_ELT(ans,0,strains); SET_VECTOR_ELT(names,0,mkChar("strains")); SET_VECTOR_ELT(ans,1,markers); SET_VECTOR_ELT(names,1,mkChar("markers")); @@ -210,13 +225,31 @@ SET_VECTOR_ELT(names,6,mkChar("chromosome")); SET_VECTOR_ELT(ans,7,family); SET_VECTOR_ELT(names,7,mkChar("family")); + */ + SET_VECTOR_ELT(ans,0,strains); + SET_STRING_ELT(names,0,mkChar("strains")); + SET_VECTOR_ELT(ans,1,markers); + SET_STRING_ELT(names,1,mkChar("markers")); + SET_VECTOR_ELT(ans,2,map); + SET_STRING_ELT(names,2,mkChar("map")); + SET_VECTOR_ELT(ans,3,subjects); + SET_STRING_ELT(names,3,mkChar("subjects")); + SET_VECTOR_ELT(ans,4,phenotypes); + SET_STRING_ELT(names,4,mkChar("phenotypes")); + SET_VECTOR_ELT(ans,5,handle); + SET_STRING_ELT(names,5,mkChar("handle")); + SET_VECTOR_ELT(ans,6,chromosome); + SET_STRING_ELT(names,6,mkChar("chromosome")); + SET_VECTOR_ELT(ans,7,family); + SET_STRING_ELT(names,7,mkChar("family")); + SET_VECTOR_ELT(ans,8,mother); + SET_STRING_ELT(names,8,mkChar("mother")); + SET_VECTOR_ELT(ans,9,father); + SET_STRING_ELT(names,9,mkChar("father")); setAttrib( ans, R_NamesSymbol, names ); - - - UNPROTECT(10); - + UNPROTECT(12); /* set the class */ @@ -286,14 +319,15 @@ QTL_DATA *read_qtl_data( FILE *fp, char *name, ALLELES *a, int verbose, int use_parents, int ped_format, char *missingCode ) { QTL_DATA *q = (QTL_DATA*)calloc(1,sizeof(QTL_DATA)); - int max_N = 3000; + int max_N = 10000; int m; - char *buffer = (char*)calloc(100000,sizeof(char)); - char **sbuf = (char**)calloc(2*a->markers,sizeof(char*)); + int bufsize = 100+a->markers*20; + char *buffer = (char*)calloc(bufsize,sizeof(char)); char **mother = NULL; char **father= NULL; double NaN = nan("char-sequence"); - + int *pcount; + int nparents=0; q->alleles = a; q->filename = (char*)strdup(name); q->N = 0; @@ -316,7 +350,6 @@ Rprintf( "Reading phenotype and genotype data from data file %s\n", name ); while ( skip_comments( fp, buffer ) != EOF ) { char *str1, *str2; - double x; m = 0; int ok = 0; if ( q->N >= max_N ) { @@ -338,14 +371,13 @@ if ( use_parents || ped_format ) { char *family = strtok( buffer, " " ); char *id = strtok( NULL, " " ); - char *mum = strtok( NULL, " " ); char *dad = strtok( NULL, " " ); + char *mum = strtok( NULL, " " ); char *sex = strtok( NULL, " " ); char *pheno = strtok( NULL, " " ); - if ( family && id && mum && dad && sex && pheno ) { - char tmpstr[1000]; + if ( family && id && dad && mum && sex && pheno ) { char *endptr; q->family[q->N] = (char*)strdup(family); father[q->N] = (char*)strdup(dad); @@ -458,18 +490,22 @@ Rprintf( "Number of individuals: %-5d\n", q->N ); Rprintf( "Number of markers: %-5d\n", q->M ); Rprintf( "Number of strains: %-5d\n", q->S ); + Rprintf( "Use Parents: %s\n", q->use_parents ? "yes" : "no" ); +#ifdef _USE_HASH_ + + Rprintf( "status %d\n", use_parents || ped_format ); if ( use_parents || ped_format ) { - int i; - ENTRY e; - ENTRY *f; + int i, nparents=0; + PARENT_KEY e; + PARENT_KEY *f; int both = 0; q->mother = (int*)calloc(q->N, sizeof(int)); q->father = (int*)calloc(q->N, sizeof(int)); - - + pcount = (int*)calloc(q->N, sizeof(int)); + if ( hcreate(q->N) == 0 ) error("Could not create hash table"); @@ -493,20 +529,86 @@ q->father[i] = (int)((long)(f->data)); else q->father[i] = -1; - if ( q->mother[i] > -1 && q->father[i] > -1 ) + if ( q->mother[i] > -1 && q->father[i] > -1 ) { both++; + pcount[q->mother[i]]++; + pcount[q->father[i]]++; + } } - + + for( i=0;iN;i++) + nparents += pcount[i]>0 hdestroy(); free(mother); free(father); + free(pcount); + + Rprintf( "Number of subjects with two parents: %-5d\n", both ); + Rprintf( "Number of parents in nuclear families: %-5d\n", nparents ); + + } +#else + + if ( use_parents || ped_format ) { + int i; + int both = 0; + PARENT_KEY *sorted = (PARENT_KEY*)calloc(q->N, sizeof(PARENT_KEY)); + PARENT_KEY *m, *f; + + q->mother = (int*)calloc(q->N, sizeof(int)); + q->father = (int*)calloc(q->N, sizeof(int)); + pcount = (int*)calloc(q->N, sizeof(int)); + + for(i=0;iN;i++) { + sorted[i].key = q->name[i]; + sorted[i].id = i; + } + + qsort( sorted, q->N, sizeof(PARENT_KEY), entrycmp ); + + + for(i=0;iN;i++) { + PARENT_KEY mum, dad; + mum.key = mother[i]; + dad.key = father[i]; + m = (PARENT_KEY*)bsearch( &mum, sorted, q->N, sizeof(PARENT_KEY), entrycmp ); + f = (PARENT_KEY*)bsearch( &dad, sorted, q->N, sizeof(PARENT_KEY), entrycmp ); + if ( m != NULL ) + q->mother[i] = m->id; + else + q->mother[i] = -1; + + if ( f != NULL ) + q->father[i] = f->id; + else + q->father[i] = -1; + if ( q->mother[i] > -1 && q->father[i] > -1 ) { + both++; + pcount[q->mother[i]]++; + pcount[q->father[i]]++; + } + + } + + for( i=0;iN;i++) + nparents += pcount[i]>0; + + free(mother); + free(father); + free(sorted); + free(pcount); Rprintf( "Number of subjects with two parents: %-5d\n", both ); + Rprintf( "Number of parents in nuclear families: %-5d\n", nparents ); + } +#endif + fit_null_qtl_model( q ); + free(buffer); return q; } @@ -534,7 +636,7 @@ for(k=0;kalleles->strains;k++) fit->trait[k] = fit->trait_error[k] = 0.0; - printf("null model mean %e var %e\n", fit->mean, fit->sigma ); + Rprintf("null model mean %e var %e\n", fit->mean, fit->sigma ); return fit->sigma; } @@ -547,7 +649,7 @@ char line[256]; int line_no = 0; int i; - int imax = 10000; + // int imax = 10000; Rprintf( "Reading subject ancestries from %s\n", filename ); skip_comments(fp, line); @@ -563,8 +665,9 @@ char *str=strtok(line," "); an->strain_name = (char**)calloc(strains,sizeof(char*)); for(k=0;kstrain_name[k] = (char*)strdup(str); + } else { Rprintf( "ERROR not enough strain names %d/%d\n", k, strains ); error("fatal HAPPY error"); @@ -760,14 +863,11 @@ char line[10000]; char marker_name[256]; - char strain_name[256]; char allele_name[256]; int markers=-1; int strains=-1; - int alleles=-1; int line_no = 0; int m, s, a; - FILE *errorf = stderr; ALLELES *A=NULL; skip_comments( fp, line ); line_no++; @@ -784,8 +884,9 @@ char *str=strtok(line," "); A->strain_name = (char**)calloc(strains,sizeof(char*)); for(k=0;kstrain_name[k] = (char*)strdup(str); + } else { Rprintf( "ERROR not enough strain names %d/%d\n", k, strains ); error("fatal HAPPY error"); @@ -971,16 +1072,12 @@ /* compute, for each individual, the prior probabilities of the different strains at the flanking markers */ - int i, j, left, right; + int i; int m = locus; - int strains = qtl->S; double total; - double Pr_ss = qtl->alleles->Pr_ss[locus]; - double Pr_st = qtl->alleles->Pr_st[locus]; double *LeftMargin = (double*)calloc(qtl->S,sizeof(double)); double *RightMargin = (double*)calloc(qtl->S,sizeof(double)); double s1 = 1.0/qtl->S; - double s12 = s1*s1; for(i=0;iN;i++) { int s, t; @@ -1076,8 +1173,8 @@ int start, stop, incr; CHROM_PAIR *genotypes = &(qtl->genos[individual]); ALLELES *A = qtl->alleles; - int m, s, t, par; - int p1, p2, offset; + int m, s, t; + int offset; double total; int markers = genotypes->markers; int strains = A->strains; @@ -1183,7 +1280,7 @@ else a = A->af[m].pr_AtoS; - int m1, m2, p1, p2; + int m1=0, m2=0, p1=0, p2=0; if ( use_parents ) { m1 = mgenotypes->chrom1[m]; m2 = mgenotypes->chrom2[m]; @@ -1302,7 +1399,7 @@ { QTL_DATA *q = NULL; - int id; + int id=0; *locus = -1; @@ -1424,19 +1521,21 @@ return Matrix; } -SEXP happyprobs2 ( SEXP handle, SEXP marker ) { - /* returns a matrix giving for each row the probability that the corresponding individual is descened form a pair of strains. happyprobs2 differs from happyprobs only in the row of probabilities return is complete - ie due to symmentry it contains the off-diagonal elements twice */ +SEXP happyprobs2 ( SEXP handle, SEXP marker, SEXP symmetrize ) { + /* returns a matrix giving for each row the probability that the corresponding individual is descened form a pair of strains. happyprobs2 differs from happyprobs only in the row of probabilities return is complete - ie due to symmetry it contains the off-diagonal elements twice. this is useful if the probabilities have been computed using parental information to estimate the phase of the haplotypes. */ int locus = -1; - QTL_DATA *q = validateParams( handle, marker, &locus, 0 ); + QTL_DATA *q = validateParams( handle, marker, &locus, 1 ); SEXP Matrix = R_NilValue; + int Symmetrize; + if ( ! isNumeric(symmetrize) || length(symmetrize) != 1 ) + error( "symmetrize is not numeric(1)"); + Symmetrize = (int)REAL(symmetrize)[0]; if ( locus >= 0 && q->dp_matrices != NULL) { ALLELE_FREQ *af = &(q->alleles->af[locus]); QTL_PRIOR ***p; - int S2 = q->S*q->S; - int i; @@ -1445,18 +1544,36 @@ p = allocate_qtl_priors( q ); compute_qtl_priors( q, p, locus, af->prior ); - PROTECT( Matrix = allocMatrix( REALSXP, q->N, S2) ); + if ( Symmetrize ) { + int S2 = q->S*(q->S+1)/2; + PROTECT( Matrix = allocMatrix( REALSXP, q->N, S2) ); - for(i=0; iN; i++) { - int j, k,m=0; - for(j=0;jS;j++) { - for(k=0;kS;k++,m++) { - REAL(Matrix)[i+q->N*m] = p[i][j][k].prior; + for(i=0; iN; i++) { + int j, k,m=0; + for(j=0;jS;j++) { + for(k=0;kN*m] = p[i][j][k].prior + p[i][k][j].prior; + } + REAL(Matrix)[i+q->N*m] = p[i][j][j].prior; + m++; } } + UNPROTECT(1); } - UNPROTECT(1); + else { + int S2 = q->S*q->S; + PROTECT( Matrix = allocMatrix( REALSXP, q->N, S2) ); + for(i=0; iN; i++) { + int j, k,m=0; + for(j=0;jS;j++) { + for(k=0;kS;k++,m++) { + REAL(Matrix)[i+q->N*m] = p[i][j][k].prior; + } + } + } + UNPROTECT(1); + } for(i=0;iN;i++) { int s1; for(s1=0;s1S;s1++) @@ -1464,7 +1581,7 @@ free(p[i]); } free(p); - + } return Matrix; } @@ -1510,13 +1627,9 @@ int locus = -1; QTL_DATA *q = validateParams( handle, marker, &locus, 0 ); SEXP Genotype = R_NilValue; - SEXP G1 = R_NilValue; - SEXP G2 = R_NilValue; if ( locus >= 0 ) { ALLELE_FREQ *af = &(q->alleles->af[locus]); - CHROM_PAIR *genos = q->genos; - int i; PROTECT( Genotype = allocMatrix( STRSXP, q->N, 2 ) ); @@ -1541,10 +1654,8 @@ SEXP happydesign( SEXP handle, SEXP marker, SEXP model ) { - double **design=NULL; SEXP Design = R_NilValue; char *mod=NULL; - SEXP clss, h; int locus = -1; QTL_DATA *q = validateParams( handle, marker, &locus, 1 ); @@ -1595,10 +1706,10 @@ for(i=0; iN; i++) { int j, k, n=0; for(j=0;jS;j++) - REAL(Design)[i+q->N*n++] = 2*p[i][j][j].prior; + REAL(Design)[i+q->N*n++] = p[i][j][j].prior; for(j=0;jS;j++) for(k=0;kN*n] = 2*(p[i][j][k].prior + p[i][k][j].prior); + REAL(Design)[i+q->N*n] = p[i][j][k].prior + p[i][k][j].prior; } } UNPROTECT(1); @@ -1611,7 +1722,7 @@ int j, k, n=0; for(j=0;jS;j++) for(k=0;kN*n] = 2*p[i][j][k].prior; + REAL(Design)[i+q->N*n] = p[i][j][k].prior; } } UNPROTECT(1); @@ -1660,10 +1771,8 @@ /* HAPLOID GENOMES */ SEXP haploid_happydesign( SEXP handle, SEXP marker ) { - double **design=NULL; SEXP Design = R_NilValue; - char *mod=NULL; - SEXP clss, h; + // char *mod=NULL; int locus = -1; QTL_DATA *q = validateParams( handle, marker, &locus, 1 ); @@ -1688,7 +1797,7 @@ } for(i=0; iN; i++) { - int j, k; + int j; for(j=0;jS;j++) { REAL(Design)[i+q->N*j] = p[i][j].prior; } @@ -1737,8 +1846,8 @@ int start, stop, incr; CHROM_PAIR *genotypes = &(qtl->genos[individual]); ALLELES *A = qtl->alleles; - int m, s, t, par; - int p1, p2, offset; + int m, s, t; + int offset; double total; int markers = genotypes->markers; int strains = A->strains; @@ -1780,7 +1889,6 @@ int g1 = genotypes->chrom1[m]; double pr_ss = Pr_ss[m+offset]; double pr_st = Pr_st[m+offset]; - int m1, m2, p1, p2; double **a; if ( qtl->an ) a = qtl->an->pr_AtoS[individual][m]; @@ -1830,13 +1938,9 @@ /* compute, for each individual, the prior probabilities of the different strains at the flanking markers */ - int i, j, left, right; + int i; int m = locus; - int strains = qtl->S; double total; - double Pr_ss = qtl->alleles->Pr_ss[locus]; - double Pr_st = qtl->alleles->Pr_st[locus]; - double s1 = 1.0/qtl->S; double P[LINKAGE_STATES]; ALLELES *A = qtl->alleles; double d = (A->af[locus+1].position - A->af[locus].position)/100.0; @@ -1874,7 +1978,7 @@ QTL_PRIOR **allocate_haploid_qtl_priors( QTL_DATA *qtl ) { QTL_PRIOR **qp = (QTL_PRIOR**)calloc(qtl->N,sizeof(QTL_PRIOR*)); - int i, k; + int i; for(i=0;iN;i++) { qp[i] = (QTL_PRIOR*)calloc(qtl->S,sizeof(QTL_PRIOR)); @@ -1883,3 +1987,11 @@ return qp; } +int entrycmp( const void *a, const void *b ) { + const PARENT_KEY *A = (PARENT_KEY*)a; + const PARENT_KEY *B = (PARENT_KEY*)b; + + return strcmp( (const char*)A->key, (const char*)B->key ); +} + + diff -Nru r-other-mott-happy-2.1/src/cl.c r-other-mott-happy-2.4/src/cl.c --- r-other-mott-happy-2.1/src/cl.c 2007-05-09 08:40:11.000000000 +0000 +++ r-other-mott-happy-2.4/src/cl.c 1970-01-01 00:00:00.000000000 +0000 @@ -1,928 +0,0 @@ -/* Last edited: Sep 10 18:26 1996 (rmott) */ -/* COMMAND_LINE PARSING AND FILE UTILITY ROUTINES */ -/* - (C) Richard Mott, Genome Analysis Lab, ICRF - */ -#define _GNU_SOURCE - -#include -#include -#include -#include -#include -#include -#include -#include"cl.h" - -#ifdef ALPHA -extern void exit(int); -#endif - -/*command line option to suppress all comments*/ -int comment = 1; - -FILE * -openfile( char *filename, char *mode ) - - /* attempts to open filename, using mode to determine function - (read/write/append), and stops the program on failure */ -{ - FILE *fp; - char *make_legal(); - - if ( ( fp = fopen( make_legal(filename), mode ) ) == NULL ) - { - fprintf( stderr, "\n ERROR: file %s will not open for function %s\n\n", filename, mode ); - exit(1); - } - else - return fp; -} - -int -getfloat( char *format, float *variable, int argc, char **argv ) /* gets a float from command line */ - -/* format can be a string like "-a=%f", or "-a". In the former case is -looks for entries line -a=10.3 AND -a 10.3, in the latter just for -a -10.3 */ - -{ float t; char *s; - char Format[256]; - - sprintf(Format,"%g",*variable); - append_usage( format, "float", Format, 0 ); - - if ( (s = next_arg( format, argc, argv ) ) && sscanf( s, "%f", &t ) == 1 ) - { - *variable = t; - return 1; - } - - - s = format; - while(*s && *s != '=') - s++; - if ( *s == '=' ) - sprintf(Format,"%s", format); - else - sprintf(Format,"%s%s", format, "=%f"); - - while ( --argc > 0 ) - { - if ( sscanf( argv[argc], Format, &t ) == 1 ) - { - *variable = t; - return 1; - } - } - - return 0; -} - -char * -next_arg( char *format, int argc, char **argv ) - -/* scan the command lines for argument matching stub, and return the -NEXT arg - -eg if format == "-max=%d", then argv == "-max 10" returns "10" */ - -{ - char *s = cl_stub( format ); - - argc--; - while( --argc > 0 ) - { - if ( ! strcmp( s, argv[argc] ) ) - return argv[argc+1]; - } - return NULL; -} - -char * -cl_stub( char *format ) /* extract stub from format, eg "-files=%s" -> "-files" */ -{ - static char stub[256]; - char *s=stub; - - while ( *format && *format != '=' ) - *s++ = *format++; - *s = 0; - - return stub; -} - - -int -getint( char *format, int *variable, int argc, char **argv ) /* gets an int from the command line */ - -{ - int t; - char *s; - char Format[256]; - - sprintf(Format,"%d",*variable); - append_usage( format, "integer", Format, 0 ); - - if ( (s = next_arg( format, argc, argv ) ) && sscanf( s, "%d", &t ) == 1 ) - { - *variable = t; - return 1; - } - - s = format; - while(*s && *s != '=') - s++; - if ( *s == '=' ) - sprintf(Format,"%s", format); - else - sprintf(Format,"%s%s", format, "=%d"); - - while ( --argc > 0 ) - { - if ( sscanf( argv[argc], Format, &t ) == 1 ) - { - *variable = t; - return 1; - } - } - - return 0; -} - - -int -clcheck( char *format, int argc, char **argv ) /* checks if a string is on the command line*/ -{ - append_usage( format, "switch", "", 0 ); - - while ( --argc > 0 ) - { - if ( strcmp( format, argv[argc] ) == 0 ) - return 1; - } - - return 0; -} - -int -getbool( char *format, int *bool, int argc, char **argv ) - - /* checks for boolean switch, eg - - -gap=yes -gap=y -gap=1 -gap=true -gap=t -gap TRUE - -gap=no -gap=n -gap=0 -gap=false -gap=f FALSE - - format should be the string "-gap" - - */ - -{ - int arg; - - if ( *bool ) - append_usage( format, "switch", "true", 0 ); - else - append_usage( format, "switch", "false", 0 ); - - if ( getint( format, bool, argc, argv) ) - return 1; - else - { - char reply[256]; - memset(reply,0,256); - if ( getarg( format, reply, argc, argv ) ) - { - if ( ! strcasecmp( reply, "yes" ) || - ! strcasecmp( reply, "y" ) || - ! strcasecmp( reply, "true" ) || - ! strcasecmp( reply, "t" ) ) - { - *bool = 1; - return 1; - } - else if ( ! strcasecmp( reply, "no" ) || - ! strcasecmp( reply, "n" ) || - ! strcasecmp( reply, "false" ) || - ! strcasecmp( reply, "f" ) ) - { - *bool = 0; - return 1; - } - } - } - - if ( clcheck( format, argc, argv ) ) - return 1; - else - return 0; -} - -int -getboolean( char *format, int *bool, int argc, char **argv ) - -/* similar to getbool(), except that function returns 1 and sets *bool = 1 if clcheck is true also . Therefore the switch - --switch - -is equivalent to - --switch=y - --noswitch is equivalent to - --switch=n -*/ - -{ - char noformat[256]; - int Argc = argc; - int n; - - if ( *bool ) - append_usage( format, "switch", "true", 1 ); - else - append_usage( format, "switch", "false", 1 ); - - if ( *format == '-' ) - sprintf(noformat,"-no%s", &format[1] ); - else - sprintf(noformat,"-no%s", format ); - - while ( --Argc > 0 ) - { - if ( ! strcmp( format, argv[Argc] ) ) - { - *bool = 1; - return 1; - } - else if ( ! strcmp( noformat,argv[Argc]) ) - { - *bool = 0; - return 1; - } - } - -/* NB there is a slight bug here in that a command line of the form - - program -notest -test=y - - will yield FALSE - - because it searches for -notest before -test=y - - THerefore you should stick to one convention within a command-line - -*/ - - return getbool( format, bool, argc, argv ); -} - - -int -getarg( char *format, char *arg, int argc, char ** argv ) /* gets a string argument*/ -{ - char *s; - char Format[256]; - - append_usage( format, "text", arg, 0 ); - - s = format; - while(*s && *s != '=') - s++; - if ( *s == '=' ) - sprintf(Format,"%s", format); - else - sprintf(Format,"%s%s", format, "=%s"); - - if ( s = next_arg( format, argc, argv ) ) - { - strcpy( arg, s ); - return 1; - } - while ( --argc > 0 ) - { - if ( sscanf( argv[argc], Format, arg ) > 0 ) return 1; - } - return 0; -} - -int get_restricted_arg( char *format, char *arg, int argc, char ** argv, int stringc, char **strings, int *value ) { - getarg( format, arg, argc, argv ); - return legal_string( arg, strings, stringc, value ); -} - -FILE * -nextfile( char *filename, int argc, char **argv ) - - /* searches the command line for arguments not beginning with "-", and tries - to open the first one it finds for reading. The corresponding argument in the command-line list is then blanked - out so that on the next call a different file is opened. The name of the - file opened is returned in filename. On failure the function returns NULL */ -{ - FILE *fp; - int i; - - *filename = 0; - - while ( --argc > 0 ) - { - if ( argv[argc][0] != 0 && argv[argc][0] != '-' ) - { - if (fp = fopen( filename, "r" )) - { - strcpy( filename, argv[argc] ); - for(i=0;i<=strlen(argv[argc]);i++) - argv[argc][i] = 0; - return fp; - } - } - } - - return NULL; -} - -char * -extension ( char *filename, char *ext) - - /* changes the extension of filename to ext. - - e.g. - extension( "file.dat", "out" ); - - produces file.out - - extension( "file", "out"); - - produces file.out - - extension( "file", ".out"); - - produces file.out - - extension( "file.out", ""); - - produces file - - Note: filename MUST be long enough to cope ! - - Returns the modified string */ -{ - int n; - - if ( ! ext ) - return NULL; - - if ( ext[0] == '.' ) - ext++; - - n = strlen(filename); - - while ( filename[n] != '.' && n > 0 ) - n--; - - if ( filename[n] != '.' ) - { - n = strlen(filename); - filename[n] = '.'; - } - - n++; - - strcpy( &filename[n], ext ); - - if ( filename[n=strlen(filename)-1] == '.' ) - filename[n] = 0; - - return filename; -} - -char * -my_basename (char *filename) - - /* strips any directory from filename - - e.g. - my_basename("/gea/users/rmott/.cshrc"); - - produces .cshrc - - my_basename("/usr"); - - produces usr - - */ - -{ - int n, m; - - n = strlen(filename); - - while( filename[n] != '/' && n > 0 ) - n--; - - if ( filename[n] == '/' ) - n++; - - m = 0; - while ( filename[m] ) - filename[m++] = filename[n++]; - - return filename; -} - -char * -dirname( char *pathname ) - -/* strips off the filename from the path, leaving the directory */ -{ - char *s = &pathname[strlen(pathname)-1]; - - while( *s && s > pathname && *s != '/' ) - s--; - - if ( s == pathname && *s != '/' ) - strcpy(pathname, "./"); - else if ( s == pathname && *s == '/' ) - strcpy(pathname,"/"); - else - *s = 0; - - return pathname; -} - -char * -directory( char *filename, char *dir ) - - /* changes the directory part of filename to dir. - - e.g. - directory( "/home/gea/fred/file.dat", "/usr/local/" ); - - produces /usr/local/file.dat - - directory("/usr/lib", "var"); - - produces var/lib, i.e. a directory slash is inserted if necessary. - - */ - -{ - char name[512]; - - my_basename( filename ); - - if ( dir && *dir ) - { - if ( dir[strlen(dir)-1] != '/' ) - sprintf( name, "%s/%s", dir, filename ); - else - sprintf( name, "%s%s", dir, filename ); - - strcpy( filename, name ); - } - return filename; -} - -char * -rootname( char *filename ) - - /* trims off the directory and the extension from filename */ - -{ - return extension(my_basename(filename),""); -} - - -/* use compiler switch -DVAX to use the form of make_legal which converts . to _ */ - -#ifdef VAX - -char * -make_legal( char *filename ) - - /* converts files like file.in.out to file_in.out */ - -{ - char *s = filename, *t; - int n=0; - - while( *s ) - { - if ( *s == '.' ) - { - n++; - t = s; - } - s++; - } - - if ( n > 1 ) - { - t--; - while ( t != filename ) - { - if ( *t == '.' ) - *t = '_'; - t--; - } - } - return filename; -} - - -#else - -char * -make_legal( char *filename ) - -{ - return filename; -} - -#endif - -FILE * -argfile( char *format, char *mode, int argc, char **argv, char *filename ) - - /* parses the command-line using format, attempts to get a filename and open - it. If format does not match any argument then NULL is returned. If format - matchs an argument but the file will not open then an error message is printed to stderr */ - -{ - FILE *fp=NULL; - - *filename = 0; - if ( getarg( format, filename, argc, argv ) ) - { - if ( ( fp = fopen( filename, mode ) ) == NULL ) - fprintf( stderr, "\n ERROR: file %s will not open for function %s\n\n", filename, mode ); - } - if ( !strcmp( mode, "w" ) ) - append_usage( format, "Writeable File", "", 1 ); - else if ( !strcmp( mode, "r" ) ) - append_usage( format, "Readable File", "", 1 ); - else - append_usage( format, "File", "", 1 ); - - return fp; -} - - -FILE * -argfile_force( char *format, char *mode, int argc, char **argv, char *filename ) - - /* parses the command-line using format, attempts to get a filename and open - it. If format does not match any argument or if format matchs an argument - but the file will not open then the program stops*/ - -{ - FILE *fp; - - if ( getarg( format, filename, argc, argv ) ) - { - if ( ( fp = fopen( filename, mode ) ) == NULL ) - { - fprintf( stderr, "\n ERROR: file %s will not open for function %s\n\n", filename, mode ); - exit(1); - } - else - return fp; - } - else - { - fprintf( stderr, "\n ERROR: argument %s not on command line\n\n", format ); - exit(1); - } -} - - -#include -#include -#include - - -int -file_time( char *filename ) - - /* returns the time that filename was last modified, or 0 if filename - cannot be opened for reading or if the call to stat fails */ - -{ - - FILE *fp; - struct stat buf; - - if ( ! (fp = fopen(filename,"r")) ) - return 0; - else - { - fclose(fp); - if ( ! stat( filename, &buf ) ) - return (int)buf.st_mtime; - else - return 0; - } -} - -char * -file_date( char *filename ) - - /* similar to file_time() execpt that a pointer to the date (in English) */ - -{ - FILE *fp; - struct stat buf; - char date[256], *t; - - strcpy( date, "?"); - - if ( (fp = fopen(filename,"r")) ) - { - fclose(fp); - if ( ! stat( filename, &buf ) ) - { - sprintf( date, "%s", ctime(&buf.st_mtime) ); - t = date; - while ( *t ) - { - if ( *t == '\n') - *t = 0; - t++; - } - - } - } - t = date; - return t; -} - -typedef struct use_list -{ - char *format; - char *text; - char *def; - struct use_list *next; -} -USE_LIST; - -static USE_LIST *use_begin=NULL, *use_end=NULL; - -void -append_usage( char *format, char *text, char *def, int overide ) - -{ - USE_LIST *use; - - if ( ! use_begin ) - use_begin = use_end = (USE_LIST *)calloc( 1, sizeof(USE_LIST) ); - else - { - use = use_begin; - while( use ) - { - if ( ! strcmp( use->format, format ) ) - { - if ( overide ) - break; - else - return; - } - use = use->next; - } - - if ( ! use ) - use_end = use_end->next = (USE_LIST *)calloc( 1, sizeof(USE_LIST) ); - else - use_end = use; - } - -/* printf("append %s %s %s\n", format, text, def ); */ - use_end->format = (char*)strdup( format ); - use_end->text = (char*)strdup( text ); - use_end->def = (char*)strdup( def ); -} - -void -gethelp( int argc, char **argv ) -{ - - int want_help = clcheck("-help", argc, argv ); - int errors = check_usage( stderr, argc, argv ); - if ( want_help || errors ) - print_usage( argc, argv, 1 ); -} - -void -print_usage( int argc, char **argv, int stop ) - -{ - USE_LIST *u; - - u = use_begin; - - fprintf( stderr, "\nusage: %s\n", argv[0] ); - while( u ) - { - fprintf( stderr, " %-15s %-20s [ %s ]\n", u->format, u->text, u->def ); - u = u->next; - } - fprintf(stderr, "\n\n"); - if ( stop ) - exit(1); -} - -int -check_usage( FILE *fp, int argc, char **argv ) -{ - USE_LIST *u; - int errors=0; - char *arg, *noarg; - - while( --argc > 0 ) - { - arg = argv[argc]; - if ( *arg == '-' ) - { - if ( strlen(arg) > 3 && (arg[1] == 'n') && (arg[2] == 'o') ) - noarg = &arg[3]; - else - noarg = NULL; - - u = use_begin; - while( u && strncmp( argv[argc], u->format, strlen(argv[argc]) ) && ( noarg == NULL || strncmp( noarg, &u->format[1], strlen(noarg) ) ) ) - u = u->next; - - if ( ! u && ! isdigit(arg[1]) ) - { - if ( fp ) - fprintf(fp, "WARNING: unknown argument %s\n", argv[argc] ); - errors++; - } - } - } - - return errors; -} - -#define COM "-comfile=" -#define MAX_COM_DEPTH 10 -static int com_depth=0; - -int -add_commands_from_file( int argc, char **argv, int *nargc, char ***nargv ) - - /* looks for a command file argument and parses the command file, adding at to the command line */ - -{ - FILE *fp; - char comfile[256], line[256], *s, **targv = argv; - int n = argc, targc=argc; - - *nargc = targc = argc; - *nargv = targv = argv; - - if ( ++com_depth < MAX_COM_DEPTH && (fp = argfile( "-comfile=%s", "r", argc, argv, comfile ) )) - { - /* printf( "! %2d reading commands from file %s\n", com_depth, comfile );*/ - - while( skip_comments( fp, line ) != EOF ) - targc++; - - rewind(fp); - - targv = (char**)calloc( targc, sizeof(char*) ); - - while( n-- ) - { - if ( argv[n] && strncmp( COM, argv[n] , strlen(COM) ) ) - targv[n] = argv[n]; - else targv[n] = (char*)strdup(" "); - } - targc = argc; - while( skip_comments( fp, s=line ) != EOF ) - { - while( isspace(*s) ) - *s++; - targv[targc++] = (char*)strdup(s); - } - add_commands_from_file( targc, targv, nargc, nargv ); - com_depth--; - return 1; - } - else - { - com_depth--; - return 0; - } -} - -int -legal_string( char *string, char **strings, int size, int *value ) - -/* checks if string is a member os strings, and sets value to the index in the array strings -returns 1 on success and 0 on failure */ -{ - int i; - - if ( string ) - for(i=0;i n ) - return 0; - else { - char *t, *s; - while( m >= 0 ) { - if ( text[n] != suffix[m] ) - return 0; - n--; m--; - } - } - return 1; -} diff -Nru r-other-mott-happy-2.1/src/cl.h r-other-mott-happy-2.4/src/cl.h --- r-other-mott-happy-2.1/src/cl.h 2007-05-09 08:41:14.000000000 +0000 +++ r-other-mott-happy-2.4/src/cl.h 1970-01-01 00:00:00.000000000 +0000 @@ -1,155 +0,0 @@ -/* Last edited: Feb 13 13:49 1996 (rmott) */ -/* CL.H is the header file for the command-line parsing functions in CL.C */ -/* (C) Richard Mott, ICRF */ - - -FILE *openfile( char *filename, char *mode ); - - /* attempts to open filename, using mode to determine function - (read/write/append), and stops the program on failure */ - -FILE *openfile_in_searchpath(char *base, char *mode, char *searchpath, char *fullname ); -FILE *openfile_in_envpath(char *base, char *mode, char *env, char *fullname ); - -int skip_comments( FILE *file, char *string ); -int getfloat( char *format, float *variable, int argc, char **argv ); /* gets a float from command line */ -int getint( char *format, int *variable, int argc, char **argv ); /* gets an int from the command line */ -int clcheck( char *format, int argc, char **argv ); /* checks if a string is on the command line*/ -int getbool( char *format, int *bool, int argc, char **argv ); - - /* checks for boolean switch, eg - - -gap=yes -gap=y -gap=1 -gap=true -gap=t -gap TRUE - -gap=no -gap=n -gap=0 -gap=false -gap=f FALSE - - format should be the string "-gap" - - */ -int getboolean( char *format, int *bool, int argc, char **argv ); /* similar to getbool(), except that function returns 1 and sets *bool = 1 if clcheck is true also . Therefore the switch - --switch - -is equivalent to - --switch=y - -*/ -int getarg( char *format, char *arg, int argc, char **argv ); /* gets a string argument*/ - -FILE *nextfile( char *filename, int argc, char **argv ); - - /* searches the command line for arguments not beginning with "-", and tries - to open the first one it finds for reading. The corresponding argument in the command-line list is then blanked - out so that on the next call a different file is opened. The name of the - file opened is returned in filename. On failure the function returns NULL */ - -char *extension ( char *filename, char *ext); - - /* changes the extension of filename to ext. - - e.g. - extension( "file.dat", "out" ); - - produces file.out - - extension( "file", "out"); - - produces file.out - - extension( "file", ".out"); - - produces file.out - - extension( "file.out", ""); - - produces file - - Note: filename MUST be long enough to cope ! - - Returns the modified string */ - - -char *my_basename (char *filename); - - /* strips any directory from filename - - e.g. - basename("/gea/users/rmott/.cshrc"); - - produces .cshrc - - basename("/usr"); - - produces usr - - */ - -char *dirname( char *pathname ); - -/* strips off the filename from the path, leaving the directory */ - -char *directory( char *filename, char *dir ); - - /* changes the directory part of filename to dir. - - e.g. - directory( "/home/gea/fred/file.dat", "/usr/local/" ); - - produces /usr/local/file.dat - - directory("/usr/lib", "var"); - - produces var/lib, i.e. a directory slash is inserted if necessary. - - */ - -char *rootname( char *filename ); - - /* trims off the directory and the extension from filename */ - -char *make_legal( char *filename ); - -FILE *argfile( char *format, char *mode, int argc, char **argv, char *filename ); - - /* parses the command-line using format, attempts to get a filename and open - it. If format does not match any argument then NULL is returned. If format - matchs an argument but the file will not open then the program stops*/ - -FILE *argfile_force( char *format, char *mode, int argc, char **argv, char *filename ); - - /* parses the command-line using format, attempts to get a filename and open - it. If format does not match any argument or if format matchs an argument - but the file will not open then the program stops*/ - -int file_time( char *filename ); - - /* returns the time that filename was last modified, or 0 if filename - cannot be opened for reading or if the call to stat fails */ - -char *file_date( char *filename ); - - /* similar to file_time() execpt that a pointer to the date (in English) */ - -int legal_string( char *string, char **strings, int size, int *value ); - -/* checks if string is a member os strings, and sets value to the index in the array strings -returns 1 on success and 0 on failure */ - -char *next_arg( char *format, int argc, char **argv ); -char *cl_stub( char *format ); - -int add_commands_from_file( int argc, char **argv, int *nargc, char ***nargv ); -void print_usage( int argc, char **argv, int stop ); -void append_usage( char *format, char *text, char *def, int overide ); -void gethelp( int argc, char **argv ); -int check_usage( FILE *fp, int argc, char **argv ); - - -char -**split_on_separator( char *string, char separator, int *items); - -int get_restricted_arg( char *format, char *arg, int argc, char ** argv, int stringc, char **strings, int *value ); - -char *Now(void); -int EndsWith( char *text, char *suffix ); - diff -Nru r-other-mott-happy-2.1/src/cmp.c r-other-mott-happy-2.4/src/cmp.c --- r-other-mott-happy-2.1/src/cmp.c 2007-05-09 08:36:14.000000000 +0000 +++ r-other-mott-happy-2.4/src/cmp.c 2013-12-19 16:21:35.000000000 +0000 @@ -63,8 +63,9 @@ int n; while ( la && lb ) { - if (n = ( ((int)a[la--]) - ((int)b[lb--]) ) ) + if (n = ( ((int)a[la--]) - ((int)b[lb--]) ) ) { return n; + } } return la-lb; } diff -Nru r-other-mott-happy-2.1/src/happy.h r-other-mott-happy-2.4/src/happy.h --- r-other-mott-happy-2.1/src/happy.h 2007-05-18 10:29:58.000000000 +0000 +++ r-other-mott-happy-2.4/src/happy.h 2013-12-19 16:21:35.000000000 +0000 @@ -131,6 +131,13 @@ } KV; +typedef struct { + char *key; + int id;; +} +PARENT_KEY; + + /* function prototypes */ CHROM_PAIR *new_chrom_pair( int markers ); diff -Nru r-other-mott-happy-2.1/src/hbcore.c r-other-mott-happy-2.4/src/hbcore.c --- r-other-mott-happy-2.1/src/hbcore.c 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/src/hbcore.c 2013-12-19 16:21:35.000000000 +0000 @@ -0,0 +1,299 @@ +#include +#include +#include + +#include "hbcore.h" + + + +XMAT* Xdip(double **Xmat, int nrow, int ncol) { + + int i,j; + double new, dcheck, dncol, dnrow, maxH, Hbar; + double mu_avNi, sd_avNi, var_avNi; + + double *Hvec=NULL, *av_Ni=NULL; + double **newX=NULL, **cumX=NULL; + + XMAT *full=NULL; + + + // memory allocation + + full = (XMAT*)calloc(1,sizeof(XMAT)); + + newX = (double**)calloc(nrow,sizeof(double*)); + cumX = (double**)calloc(nrow,sizeof(double*)); + for (i=0; i= 1.00001) ) { + Rprintf("individual %i : dcheck = %e ERROR HMM probs do not sum to 1\n",i,dcheck); + } + + Hvec[i] = 0.0; + for (j=0; j= 1.00001) ) { + Rprintf("individual %i : dcheck = %e ERROR HMM probs do not sum to 1\n",i,dcheck); + } + + Hvec[i] = 0.0; + for (j=0; j=0; j--) { // load the shuffle table (after 8 warm-ups) + k=(*idum)/IQ1; + *idum=IA1*(*idum-k*IQ1)-k*IR1; + if (*idum < 0) *idum += IM1; + if (j < NTAB) iv[j] = *idum; + } + iy=iv[0]; + } + + k=(*idum)/IQ1; // start here when not initialising + *idum=IA1*(*idum-k*IQ1)-k*IR1; // compute idum = (IA1*idum) % IM1 without overflows, by Schrage's \ +method + if (*idum < 0) *idum += IM1; + k=idum2/IQ2; + idum2=IA2*(idum2-k*IQ2)-k*IR2; // compute idum2 = (IA2*idum) % IM2 likewise + if (idum2 < 0) idum2 += IM2; + j=iy/NDIV; // will be in the range 0..NTAB-1 + iy=iv[j]-idum2; // here idum is shuffled, idum and idum2 are combined to generate o\ +utput + iv[j] = *idum; + if (iy < 1) iy += IMM1; + if ((temp=AM*iy) > RNMX) return RNMX; + else return temp; // because users don't expect endpoint values + +} + diff -Nru r-other-mott-happy-2.1/src/hbcore.h r-other-mott-happy-2.4/src/hbcore.h --- r-other-mott-happy-2.1/src/hbcore.h 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/src/hbcore.h 2013-12-19 16:21:35.000000000 +0000 @@ -0,0 +1,24 @@ +#include +#include // has exit() function +#include +#include +#include + + +typedef struct { + double **X, **cumX; + double *Hvec, *av_Ni; + double Hbar, muNi, sdNi; +} XMAT; + + + +XMAT* Xdip(double **Xmat, int nrow, int ncol); +XMAT* Xhap(double **Xmat, int nrow, int ncol); + + +double NRroundit(double d, int dig); + +float ran1(long *idum); +float ran2(long *idum); +void NRsort(int nr, double *arr); diff -Nru r-other-mott-happy-2.1/src/hbrem.c r-other-mott-happy-2.4/src/hbrem.c --- r-other-mott-happy-2.1/src/hbrem.c 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/src/hbrem.c 2013-12-19 16:21:35.000000000 +0000 @@ -0,0 +1,354 @@ + + +#define _ISOC99_SOURCE +#define _GNU_SOURCE + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include"happy.h" +/*#include"cl.h" +#include"cmp.h" +#include"stats.h" +*/ + +#include "hbcore.h" +#include "hbrem.h" + + + +SEXP hbrem( SEXP RX, SEXP HaploidInd, SEXP Ndip, SEXP Nstrain, SEXP Nind, SEXP Npost, SEXP Nbin, SEXP Ry ) { + + + SEXP Rstats = R_NilValue; + SEXP Reffectmeans = R_NilValue; + SEXP ReffectSDs = R_NilValue; + SEXP ReffectNis = R_NilValue; + SEXP Rstrainmeans = R_NilValue; + SEXP list = R_NilValue; + + int haploid_ind,nbin,ncol,nrow,nsim,nstrains; + + double *y=NULL; + double *pntRstats=NULL, *pntReffectmeans=NULL, *pntReffectSDs=NULL; + double *pntReffectNis=NULL, *pntRX=NULL, *pntRstrainmeans=NULL; + + + PROTECT(RX = AS_NUMERIC(RX)); + pntRX = (double*)NUMERIC_POINTER(RX); + + PROTECT(HaploidInd = AS_INTEGER(HaploidInd)); + haploid_ind = (int)INTEGER_POINTER(HaploidInd)[0]; + + PROTECT(Ndip = AS_INTEGER(Ndip)); + ncol = (int)INTEGER_POINTER(Ndip)[0]; + + PROTECT(Nstrain = AS_INTEGER(Nstrain)); + nstrains = (int)INTEGER_POINTER(Nstrain)[0]; + + PROTECT(Nind = AS_INTEGER(Nind)); + nrow = (int)INTEGER_POINTER(Nind)[0]; + + PROTECT(Npost = AS_INTEGER(Npost)); + nsim = (int)INTEGER_POINTER(Npost)[0]; + + PROTECT(Nbin = AS_INTEGER(Nbin)); + nbin = (int)INTEGER_POINTER(Nbin)[0]; + + /* Rprintf( "nbin %d npost %d nrow %d ncol %d haploid_ind %d\n", nbin, nsim, nrow, ncol, haploid_ind); */ + PROTECT(Ry = AS_NUMERIC(Ry)); + y = (double*)NUMERIC_POINTER(Ry); + + PROTECT(Rstats = NEW_NUMERIC(45)); // allocates storage space + pntRstats = (double*)NUMERIC_POINTER(Rstats); + + PROTECT(Reffectmeans = NEW_NUMERIC(ncol)); // allocates storage space + pntReffectmeans = (double*)NUMERIC_POINTER(Reffectmeans); + + PROTECT(ReffectSDs = NEW_NUMERIC(ncol)); // allocates storage space + pntReffectSDs = (double*)NUMERIC_POINTER(ReffectSDs); + + PROTECT(ReffectNis = NEW_NUMERIC(ncol)); // allocates storage space + pntReffectNis = (double*)NUMERIC_POINTER(ReffectNis); + + PROTECT(Rstrainmeans = NEW_NUMERIC(nstrains)); + pntRstrainmeans = (double*)NUMERIC_POINTER(Rstrainmeans); + + + int i,j,inc; + long idum=0; + float init; + + double **Xmat=NULL; + double *Tcomb=NULL; + + XMAT *xmat=NULL; + POST *SLsamp=NULL; + SSTA *SLstats=NULL; + + + Xmat = (double**)calloc(nrow,sizeof(double*)); + for (j=0; j +#include +#include +#include #include -#include"cl.h" +#include"readline.h" /* functions for reading in lines (C) Richard Mott,, ICRF */ -void -uncomment( char *string ) /* truncates string at the first ! */ +void uncomment( char *string ) /* truncates string at the first ! */ { while ( *string != '!' && *string != '#' && *string != 0 ) string++; *string = 0; } -int -read_line( FILE *file, char *string ) { +int read_line( FILE *file, char *string ) { int c; int i=0; - if ( file ) { - while(c=getc(file)) { + if ( file != NULL ) { + while((c=getc(file))) { if (!i && c==EOF) return EOF; if (i && c==EOF ) return i; if (c=='\n') return i; @@ -30,11 +30,10 @@ return EOF; } -int -next_line( FILE *file ) { +int next_line( FILE *file ) { int c; - if ( file ) { - while(c=getc(file)) { + if ( file != NULL ) { + while((c=getc(file))) { if (feof(file)) return 0; if (c=='\n') return 1; } @@ -42,8 +41,7 @@ return EOF; } -int -not_blank( char *string ) /* checks whether string is full of white space */ +int not_blank( char *string ) /* checks whether string is full of white space */ { while ( *string != 0 ) { if ( ! isspace(*string) ) return 1; @@ -56,8 +54,7 @@ /* reads in successive lines, truncating comments and skipping blank lines */ -int -skip_comments( FILE *file, char *string ) { +int skip_comments( FILE *file, char *string ) { int n = EOF; *string = 0; @@ -70,3 +67,23 @@ } return n; } + +int +legal_string( char *string, char **strings, int size, int *value ) + +/* checks if string is a member os strings, and sets value to the index in the array strings +returns 1 on success and 0 on failure */ +{ + int i; + + if ( string ) + for(i=0;i -#include #include +#include +#include +#include +#include +#include +#include #include"stats.h" #include"cmp.h" @@ -24,7 +32,7 @@ double lin_regression( double *x, double *y, int from, int to, double *intercept, double *slope, double *sigma, double *t_slope, double *stderr_slope, double *stderr_intercept ) { double s_x, s_y, ss_x, ss_y, ss_xy; - int i, k, n; + int k; double N=to-from+1, R; s_x = s_y = ss_x = ss_y = ss_xy = 0.0; @@ -59,7 +67,7 @@ int len = stop-start+1; double *rank = (double*)calloc( len, sizeof(double) ); double **ptr = (double**)calloc( len, sizeof(double*) ); - int n, m; + int n; for(n=0;n 0 ) - p++; - else - q++; - } - p /= len; - q /= len; - - pp = pq = qp = qq = 0.0; - for(n=1;n 0 ) { - if ( diff[n-1]>0 ) - pp++; - else - pq++; - } - else { - if ( diff[n-1]>0 ) - qp++; - else - qq++; - } - } - - len1 = len-1; - - e = q*q*len1; z = ( qq - e ); chi = z*z/e; - e = p*q*len1; z = ( pq - e ); chi += z*z/e; - e = q*p*len1; z = ( qp - e ); chi += z*z/e; - e = p*p*len1; z = ( pp - e ); chi += z*z/e; - printf("chisq = %.3f\n", chi ); - } - - - free(diff); - - - return dw; -} - - double normal_tail( double z ) { diff -Nru r-other-mott-happy-2.1/src/subrout.c r-other-mott-happy-2.4/src/subrout.c --- r-other-mott-happy-2.1/src/subrout.c 1970-01-01 00:00:00.000000000 +0000 +++ r-other-mott-happy-2.4/src/subrout.c 2013-12-19 16:21:35.000000000 +0000 @@ -0,0 +1,1793 @@ +#include +#include +#include + +#include"hbcore.h" +#include"hbrem.h" +#include"cmp.h" + +POST* single_locus_jointpostX(XMAT *xmat,double *y,int nsim,int ncol,int nrow,int nbin,long *idum) { + + FILE *out=NULL; + char check; + char sampoutname[25]="SLXsamp_R.dat"; + + int i,sim; + double nu,kT; + + int **NiX=NULL; + + double *postkT=NULL, *postvar=NULL, *postmu=NULL, *nullvar=NULL, *nullmu=NULL; + double *Lnull=NULL, *Lqtl=NULL; + + double **postT=NULL, **yTbarX=NULL; + + SIMX *selX=NULL; + GRKT *kTdist=NULL; + POST *single=NULL; + + + single = (POST*)calloc(1,sizeof(POST)); + + postkT = (double*)calloc(nsim,sizeof(double)); + postvar = (double*)calloc(nsim,sizeof(double)); + postmu = (double*)calloc(nsim,sizeof(double)); + + postT = (double**)calloc(nsim,sizeof(double*)); + yTbarX = (double**)calloc(nsim,sizeof(double*)); + for (i=0; i= nbin ) { + dn = (double)(*trueX).Ni[ indvec[i][j] ]; + + sum2 = ( sum2 + dn ); + sum1 = ( sum1 + dn*(*SLstats).muT[ indvec[i][j] ] ); + } + } + Tcomb[i] = sum1/sum2; + } + + + for (i=0; i 0.0 ) { + dn = (*xmat).av_Ni[ indvec[i][j] ]; + + sum2 = ( sum2 + dn ); + sum1 = ( sum1 + dn*(*SLstats).muT[ indvec[i][j] ] ); + } + } + Tcomb[i] = sum1/sum2; + } + + + for (i=0; i= nbin ) { + muT[j] = muT[j] + (*SLsamp).T[i][j]; + } + } + + mu_nullmu = mu_nullmu + (*SLsamp).null_mu[i]; + mu_nullvar = mu_nullvar + (*SLsamp).null_var[i]; + + avlik_null = avlik_null + (*SLsamp).lik_null[i]; + avlik_qtl = avlik_qtl + (*SLsamp).lik_qtl[i]; + } + mukT = mukT/dnsim; + mumu = mumu/dnsim; + muvar = muvar/dnsim; + for (j=0; j= nbin ) { + muT[j] = muT[j]/dnsim; + } + } + mu_nullmu = mu_nullmu/dnsim; + mu_nullvar = mu_nullvar/dnsim; + avlik_null = avlik_null/dnsim; + avlik_qtl = avlik_qtl/dnsim; + + modevar = (((*SLsamp).N - 3.0)/((*SLsamp).N + 1.0))*muvar; + mode_nullvar = (((*SLsamp).N - 3.0)/((*SLsamp).N + 1.0))*mu_nullvar; + + varkT = 0.0; + for (i=0; i= nbin ) { + sdT[j] = sdT[j] + ((*SLsamp).T[i][j] - muT[j])*((*SLsamp).T[i][j] - muT[j]); + } + } + } + varkT = varkT/dnsim; + for (j=0; j= nbin ) { + sdT[j] = sdT[j]/dnsim; + sdT[j] = sqrt(sdT[j]); + } + } + + ga = ((mukT*mukT*(1.0 - mukT)/varkT) - mukT); + gb = ga*((1.0 - mukT)/mukT); + modekT = (ga - 1.0)/(ga + gb - 2.0); + if (modekT < 0.0) { + modekT = 0.0; + } else if (modekT > 1.0) { + modekT = 1.0; + } + + + // plug-ins for pD + + Lnull = null_lik(trueX,y,mode_nullvar,mu_nullmu,nrow,nbin); + Lqtl = qtl_lik(trueX,y,modekT,modevar,mumu,muT,nrow,nbin); + + pD_qtl = -2.0*(avlik_qtl - Lqtl); + pD_null = -2.0*(avlik_null - Lnull); + + DIC_qtl = ( pD_qtl - 2.0*avlik_qtl ); + DIC_null = ( pD_null - 2.0*avlik_null ); + + DIC_diff = ( DIC_null - DIC_qtl ); + + // plug-ins for BIC + + Lfocqtl = qtl_Lfoc(trueX,y,modekT,modevar,mumu,nrow,ncol,nbin); + + BIC_qtl = ( (3.0*log((*SLsamp).N)) - (2.0*Lfocqtl) ); + BIC_null = ( (2.0*log((*SLsamp).N)) - (2.0*Lnull) ); + + BF = exp(-(BIC_null - BIC_qtl)/2.0); + logBF = -log10(BF); + + + count1 = 0; + for (i=0; i count1 ) { + count1 = kThst[i]; + } + } + + + // sort posterior sample + // NRsort doesn't use 0 index - entries start at 1 + + qsort( stkT+1, nsim, sizeof(double), dcmp ); + qsort( stmu+1, nsim, sizeof(double), dcmp ); + qsort( stvar+1, nsim, sizeof(double), dcmp ); + +// NRsort(nsim,stkT); +// NRsort(nsim,stmu); +// NRsort(nsim,stvar); + + + // calculate medians + + index = (int)(0.5*((float)nsim)); + cred = stkT[index] + stkT[index + 1]; + medkT = cred/2.0; + cred = stmu[index] + stmu[index + 1]; + medmu = cred/2.0; + cred = stvar[index] + stvar[index + 1]; + medvar = cred/2.0; + + + // calculate HPD credible intervals + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(99*nsim/100)) < nsim ) { + bkT = stkT[i]; + tkT = stkT[i+(99*nsim/100)]; + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTb99 = bkT; + kTt99 = tkT; + } + bmu = stmu[i]; + tmu = stmu[i+(99*nsim/100)]; + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub99 = bmu; + mut99 = tmu; + } + bvar = stvar[i]; + tvar = stvar[i+(99*nsim/100)]; + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb99 = bvar; + vart99 = tvar; + } + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(95*nsim/100)) < nsim ) { + bkT = stkT[i]; + tkT = stkT[i+(95*nsim/100)]; + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt95 = tkT; + kTb95 = bkT; + } + bmu = stmu[i]; + tmu = stmu[i+(95*nsim/100)]; + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub95 = bmu; + mut95 = tmu; + } + bvar = stvar[i]; + tvar = stvar[i+(95*nsim/100)]; + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb95 = bvar; + vart95 = tvar; + } + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(75*nsim/100)) < nsim ) { + bkT = stkT[i]; + tkT = stkT[i+(75*nsim/100)]; + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt75 = tkT; + kTb75 = bkT; + } + bmu = stmu[i]; + tmu = stmu[i+(75*nsim/100)]; + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub75 = bmu; + mut75 = tmu; + } + bvar = stvar[i]; + tvar = stvar[i+(75*nsim/100)]; + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb75 = bvar; + vart75 = tvar; + } + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(50*nsim/100)) < nsim ) { + bkT = stkT[i]; + tkT = stkT[i+(50*nsim/100)]; + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt50 = tkT; + kTb50 = bkT; + } + bmu = stmu[i]; + tmu = stmu[i+(50*nsim/100)]; + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub50 = bmu; + mut50 = tmu; + } + bvar = stvar[i]; + tvar = stvar[i+(50*nsim/100)]; + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb50 = bvar; + vart50 = tvar; + } + i = i + 1; + } + + + + HPDlkT[3] = kTb99; + HPDukT[3] = kTt99; + HPDlkT[2] = kTb95; + HPDukT[2] = kTt95; + HPDlkT[1] = kTb75; + HPDukT[1] = kTt75; + HPDlkT[0] = kTb50; + HPDukT[0] = kTt50; + + HPDlmu[3] = mub99; + HPDumu[3] = mut99; + HPDlmu[2] = mub95; + HPDumu[2] = mut95; + HPDlmu[1] = mub75; + HPDumu[1] = mut75; + HPDlmu[0] = mub50; + HPDumu[0] = mut50; + + HPDlvar[3] = varb99; + HPDuvar[3] = vart99; + HPDlvar[2] = varb95; + HPDuvar[2] = vart95; + HPDlvar[1] = varb75; + HPDuvar[1] = vart75; + HPDlvar[0] = varb50; + HPDuvar[0] = vart50; + + + free(stkT); + free(stmu); + free(stvar); + + (*stats).kThist = kThst; + (*stats).HPDlkT = HPDlkT; + (*stats).HPDukT = HPDukT; + (*stats).modekT = modekT; + (*stats).mukT = mukT; + (*stats).medkT = medkT; + (*stats).ga = ga; + (*stats).gb = gb; + + (*stats).HPDlmu = HPDlmu; + (*stats).HPDumu = HPDumu; + (*stats).mumu = mumu; + (*stats).medmu = medmu; + + (*stats).HPDlvar = HPDlvar; + (*stats).HPDuvar = HPDuvar; + (*stats).modevar = modevar; + (*stats).muvar = muvar; + (*stats).medvar = medvar; + + (*stats).muT = muT; + (*stats).sdT = sdT; + + (*stats).pD_null = pD_null; + (*stats).pD_qtl = pD_qtl; + (*stats).DIC_null = DIC_null; + (*stats).DIC_qtl = DIC_qtl; + (*stats).DIC_diff = DIC_diff; + (*stats).BIC_null = BIC_null; + (*stats).BIC_qtl = BIC_qtl; + (*stats).BF = BF; + (*stats).logBF = logBF; + + + return stats; + +} + + +SSTA* single_locus_sumstatsX(XMAT *xmat, POST *SLsamp, double *y, int nsim, int ncol, int nrow, int nbin) { + + int i,j,index,count1,nprc,sum; + + double cred,dnbin,dnrow,dnsim,prc; + + double mukT,medkT,mumu,medmu,muvar,medvar,modevar,modekT; + double BF,BIC_qtl,BIC_null,logBF,DIC_null,DIC_qtl,DIC_diff; + double avlik_qtl,avlik_null,Lfocqtl,Lqtl,Lnull,pD_null,pD_qtl; + double ga,gb,varkT,mu_nullmu,mu_nullvar,mode_nullvar,av_ysq,av_ybar; + + double bkT,bmu,bvar,tkT,tmu,tvar,kTdiff,mudiff,vardiff; + double kTb99,kTt99,kTb95,kTt95,kTb75,kTt75,kTb50,kTt50; + double mub99,mut99,mub95,mut95,mub75,mut75,mub50,mut50; + double varb99,vart99,varb95,vart95,varb75,vart75,varb50,vart50; + + int *kThst=NULL; + + double *stkT=NULL,*HPDlkT=NULL,*HPDukT=NULL; + double *stmu=NULL,*HPDlmu=NULL,*HPDumu=NULL; + double *stvar=NULL,*HPDlvar=NULL,*HPDuvar=NULL; + + double *muT=NULL,*sdT=NULL,*nT=NULL,*av_yT=NULL; + + SSTA *stats=NULL; + + + nprc = 201; + prc = 200.0; + + + stats = (SSTA*)calloc(1,sizeof(SSTA)); + + HPDlkT = (double*)calloc(4,sizeof(double)); + HPDukT = (double*)calloc(4,sizeof(double)); + + HPDlmu = (double*)calloc(4,sizeof(double)); + HPDumu = (double*)calloc(4,sizeof(double)); + + HPDlvar = (double*)calloc(4,sizeof(double)); + HPDuvar = (double*)calloc(4,sizeof(double)); + + kThst = (int*)calloc(nprc,sizeof(int)); + + stkT = (double*)calloc((nsim+1),sizeof(double)); + stmu = (double*)calloc((nsim+1),sizeof(double)); + stvar = (double*)calloc((nsim+1),sizeof(double)); + + muT = (double*)calloc(ncol,sizeof(double)); + sdT = (double*)calloc(ncol,sizeof(double)); + nT = (double*)calloc(ncol, sizeof(double)); + av_yT = (double*)calloc(ncol,sizeof(double)); + + + + dnsim = (double)nsim; + nsim = (int)nsim; + + dnrow = (double)nrow; + nrow = (int)nrow; + + dnbin = (double)nbin; + nbin = (int)nbin; + + + mukT = 0.0; + mumu = 0.0; + muvar = 0.0; + mu_nullmu = 0.0; + mu_nullvar = 0.0; + avlik_null = 0.0; + avlik_qtl = 0.0; + + for (i=0; i= nbin ) { + muT[j] = muT[j] + (*SLsamp).T[i][j]; + nT[j] = nT[j] + 1.0; + } + } + + mu_nullmu = mu_nullmu + (*SLsamp).null_mu[i]; + mu_nullvar = mu_nullvar + (*SLsamp).null_var[i]; + + avlik_null = avlik_null + (*SLsamp).lik_null[i]; + avlik_qtl = avlik_qtl + (*SLsamp).lik_qtl[i]; + } + + mukT = mukT/dnsim; + mumu = mumu/dnsim; + muvar = muvar/dnsim; + for(j=0; j 0.0 ) { + muT[j] = muT[j]/nT[j]; + } + } + mu_nullmu = mu_nullmu/dnsim; + mu_nullvar = mu_nullvar/dnsim; + avlik_null = avlik_null/dnsim; + avlik_qtl = avlik_qtl/dnsim; + + + modevar = ((dnrow - 3.0)/(dnrow + 1.0))*muvar; + mode_nullvar = ((dnrow - 3.0)/(dnrow + 1.0))*mu_nullvar; + + varkT = 0.0; + for (i=0; i= nbin ) { + sdT[j] = sdT[j] + ((*SLsamp).T[i][j] - muT[j])*((*SLsamp).T[i][j] - muT[j]); + } + } + } + varkT = varkT/dnsim; + for (j=0; j 0.0 ) { + sdT[j] = sdT[j]/nT[j]; + sdT[j] = sqrt(sdT[j]); + } + } + + ga = ((mukT*mukT*(1.0 - mukT)/varkT) - mukT); + gb = ga*((1.0 - mukT)/mukT); + modekT = (ga - 1.0)/(ga + gb - 2.0); + if (modekT < 0.0) { + modekT = 0.0; + } else if (modekT > 1.0) { + modekT = 1.0; + } + + + // plug-ins for pD + + av_ysq = 0.0; + av_ybar = 0.0; + for (i=0; i 0.0 ) { + av_yT[j] = av_yT[j]/(*xmat).av_Ni[j]; + } + } + + Lnull = null_plug((*xmat).av_Ni,av_ysq,av_ybar,mode_nullvar,mu_nullmu,ncol,nrow); + Lqtl = qtl_plug(av_yT,(*xmat).av_Ni,av_ysq,modekT,modevar,mumu,muT,ncol,nrow); + + pD_qtl = -2.0*(avlik_qtl - Lqtl); + pD_null = -2.0*(avlik_null - Lnull); + + DIC_qtl = ( pD_qtl - 2.0*avlik_qtl ); + DIC_null = ( pD_null - 2.0*avlik_null ); + + DIC_diff = ( DIC_null - DIC_qtl ); + + // plug-ins for BIC + + Lfocqtl = qtl_LfocX(av_yT,(*xmat).av_Ni,av_ysq,av_ybar,modekT,modevar,mumu,ncol,nrow); + + BIC_qtl = ( (3.0*log(dnrow)) - (2.0*Lfocqtl) ); + BIC_null = ( (2.0*log(dnrow)) - (2.0*Lnull) ); + + BF = exp(-(BIC_null - BIC_qtl)/2.0); + logBF = -log10(BF); + + + count1 = 0; + for (i=0; i count1 ) { + count1 = kThst[i]; + } + } + + + // sort posterior sample + // NRsort doesn't use 0 index - entries start at 1 + + qsort( stkT+1, nsim, sizeof(double), dcmp ); + qsort( stmu+1, nsim, sizeof(double), dcmp ); + qsort( stvar+1, nsim, sizeof(double), dcmp ); + +// NRsort(nsim,stkT); +// NRsort(nsim,stmu); +// NRsort(nsim,stvar); + + + // calculate medians + + index = (int)(0.5*((float)nsim)); + nsim = (int)nsim; + + cred = stkT[index] + stkT[index + 1]; + medkT = cred/2.0; + cred = stmu[index] + stmu[index + 1]; + medmu = cred/2.0; + cred = stvar[index] + stvar[index + 1]; + medvar = cred/2.0; + + + // calculate HPD credible intervals + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(99*nsim/100)) < nsim ) { + + bkT = stkT[i]; + tkT = stkT[i+(99*nsim/100)]; + + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTb99 = bkT; + kTt99 = tkT; + } + + bmu = stmu[i]; + tmu = stmu[i+(99*nsim/100)]; + + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub99 = bmu; + mut99 = tmu; + } + + bvar = stvar[i]; + tvar = stvar[i+(99*nsim/100)]; + + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb99 = bvar; + vart99 = tvar; + } + + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(95*nsim/100)) < nsim ) { + + bkT = stkT[i]; + tkT = stkT[i+(95*nsim/100)]; + + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt95 = tkT; + kTb95 = bkT; + } + + bmu = stmu[i]; + tmu = stmu[i+(95*nsim/100)]; + + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub95 = bmu; + mut95 = tmu; + } + + bvar = stvar[i]; + tvar = stvar[i+(95*nsim/100)]; + + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb95 = bvar; + vart95 = tvar; + } + + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(75*nsim/100)) < nsim ) { + + bkT = stkT[i]; + tkT = stkT[i+(75*nsim/100)]; + + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt75 = tkT; + kTb75 = bkT; + } + + bmu = stmu[i]; + tmu = stmu[i+(75*nsim/100)]; + + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub75 = bmu; + mut75 = tmu; + } + + bvar = stvar[i]; + tvar = stvar[i+(75*nsim/100)]; + + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb75 = bvar; + vart75 = tvar; + } + + i = i + 1; + } + + + kTdiff = 1.0; + mudiff = ( stmu[nsim] - stmu[1] ); + vardiff = ( stvar[nsim] - stvar[1] ); + + i = 1; + while ( (i+(50*nsim/100)) < nsim ) { + + bkT = stkT[i]; + tkT = stkT[i+(50*nsim/100)]; + + if ( (tkT - bkT) < kTdiff ) { + kTdiff = (tkT - bkT); + kTt50 = tkT; + kTb50 = bkT; + } + + bmu = stmu[i]; + tmu = stmu[i+(50*nsim/100)]; + + if ( (tmu - bmu) < mudiff ) { + mudiff = (tmu - bmu); + mub50 = bmu; + mut50 = tmu; + } + + bvar = stvar[i]; + tvar = stvar[i+(50*nsim/100)]; + + if ( (tvar - bvar) < vardiff ) { + vardiff = (tvar - bvar); + varb50 = bvar; + vart50 = tvar; + } + + i = i + 1; + } + + + HPDlkT[3] = kTb99; + HPDukT[3] = kTt99; + HPDlkT[2] = kTb95; + HPDukT[2] = kTt95; + HPDlkT[1] = kTb75; + HPDukT[1] = kTt75; + HPDlkT[0] = kTb50; + HPDukT[0] = kTt50; + + HPDlmu[3] = mub99; + HPDumu[3] = mut99; + HPDlmu[2] = mub95; + HPDumu[2] = mut95; + HPDlmu[1] = mub75; + HPDumu[1] = mut75; + HPDlmu[0] = mub50; + HPDumu[0] = mut50; + + HPDlvar[3] = varb99; + HPDuvar[3] = vart99; + HPDlvar[2] = varb95; + HPDuvar[2] = vart95; + HPDlvar[1] = varb75; + HPDuvar[1] = vart75; + HPDlvar[0] = varb50; + HPDuvar[0] = vart50; + + + free(stkT); + free(stmu); + free(stvar); + free(nT); + free(av_yT); + + + (*stats).kThist = kThst; + (*stats).HPDlkT = HPDlkT; + (*stats).HPDukT = HPDukT; + (*stats).modekT = modekT; + (*stats).mukT = mukT; + (*stats).medkT = medkT; + (*stats).ga = ga; + (*stats).gb = gb; + + (*stats).HPDlmu = HPDlmu; + (*stats).HPDumu = HPDumu; + (*stats).mumu = mumu; + (*stats).medmu = medmu; + + (*stats).HPDlvar = HPDlvar; + (*stats).HPDuvar = HPDuvar; + (*stats).modevar = modevar; + (*stats).muvar = muvar; + (*stats).medvar = medvar; + + (*stats).muT = muT; + (*stats).sdT = sdT; + + (*stats).pD_null = pD_null; + (*stats).pD_qtl = pD_qtl; + (*stats).DIC_null = DIC_null; + (*stats).DIC_qtl = DIC_qtl; + (*stats).DIC_diff = DIC_diff; + (*stats).BIC_null = BIC_null; + (*stats).BIC_qtl = BIC_qtl; + (*stats).BF = BF; + (*stats).logBF = logBF; + + + return stats; + +} + + +double qtl_LfocX(double *av_yT, double *avNi, double av_ysq, double av_ybar, double kT, double var, double mu, int ncol, int nrow) { + + // integrated version of qtl_plug, compared to null_plug + + int i; + double deno,dncol,dnrow,lik,prod1,sum1,sum2,sum3; + + + dnrow = (double)nrow; + nrow = (int)nrow; + + dncol = (double)ncol; + ncol = (int)ncol; + + + sum2 = 0.0; + prod1 = 0.0; + for (i=0; i 0.0 ) { + deno = ( 1.0 - kT + kT*avNi[i] ); + prod1 = prod1 + log(deno); + sum2 = sum2 + ((avNi[i]*avNi[i]*(av_yT[i] - mu)*(av_yT[i] - mu))/deno); + } + } + prod1 = prod1/2.0; + + sum1 = ( (dnrow*mu*(mu - 2.0*av_ybar)) + av_ysq ); + sum3 = (sum1 - (kT*sum2))/(2.0*var*(1.0 - kT)); + + lik = ( (-dnrow/2.0)*log(2.0*M_PI) + ((dncol - dnrow)/2.0)*log(1.0 - kT) - (dnrow/2.0)*log(var) - prod1 - sum3 ); + + return lik; + +} + + +double qtl_plug(double *av_yT, double *avNi, double av_ysq, double kT, double var, double mu, double *T, int ncol, int nrow) { + + int i; + double dnrow,lik,sum2; + + + dnrow = (double)nrow; + nrow = (int)nrow; + + + sum2 = 0.0; + for (i=0; i 0.0 ) { + sum2 = sum2 + ( (avNi[i])*(mu + T[i])*(mu + T[i] - 2.0*av_yT[i]) ); + } + } + sum2 = sum2 + av_ysq; + + lik = ( ((-dnrow/2.0)*log(2.0*M_PI)) - ((dnrow/2.0)*log(1.0 - kT)) - ((dnrow/2.0)*log(var)) - (sum2/(2.0*var*(1.0 - kT))) ); + + return lik; + +} + + +double qtl_lik(SIMX *trueX, double *y, double kT, double var, double mu, double *T, int nrow, int nbin) { + + int i,m; + double Nsamp,lik,sum2; + + Nsamp = 0.0; + sum2 = 0.0; + for (i=0; i= nbin ) { + sum2 = sum2 + (y[i] - mu - T[m])*(y[i] - mu - T[m]); + Nsamp = Nsamp + 1.0; + } + } + + lik = ( ((-Nsamp/2.0)*log(2.0*M_PI)) - ((Nsamp/2.0)*log(1.0 - kT)) - ((Nsamp/2.0)*log(var)) - (sum2/(2.0*var*(1.0 - kT))) ); + + + return lik; + +} + + +double qtl_Lfoc(SIMX *trueX, double *y, double kT, double var, double mu, int nrow, int ncol, int nbin) { + + int i,m; + double deno,dn,Nsamp,lik,prod1,sum1,sum2,sum3,tsamp; + + double *yTbar=NULL; + + yTbar = (double*)calloc(ncol,sizeof(double)); + + Nsamp = 0.0; + sum1 = 0.0; + for (i=0; i= nbin ) { + sum1 = sum1 + ((y[i] - mu)*(y[i] - mu)); + yTbar[m] = yTbar[m] + y[i]; + Nsamp = Nsamp + 1.0; + } + } + + tsamp = 0.0; + prod1 = 0.0; + sum2 = 0.0; + for (i=0; i= nbin ) { + dn = (double)(*trueX).Ni[i]; + deno = ( 1.0 - kT + kT*dn ); + tsamp = tsamp + 1.0; + + prod1 = prod1 + log(deno); + yTbar[i] = yTbar[i]/dn; + + sum2 = sum2 + ((dn*dn*(yTbar[i] - mu)*(yTbar[i] - mu))/deno); + } + } + prod1 = prod1/2.0; + + sum3 = (sum1 - (kT*sum2))/(2.0*var*(1.0 - kT)); + + lik = ( (-Nsamp/2.0)*log(2.0*M_PI) + ((tsamp - Nsamp)/2.0)*log(1.0 - kT) - (Nsamp/2.0)*log(var) - prod1 - sum3 ); + + free(yTbar); + + return lik; +} + + +double null_plug(double *avNi, double av_ysq, double av_ybar, double var, double mu, int ncol, int nrow) { + + int i; + double dnrow,lik,sum1; + + dnrow = (double)nrow; + nrow = (int)nrow; + + sum1 = ( (dnrow*mu*(mu - 2.0*av_ybar)) + av_ysq ); + + lik = ( ((-dnrow/2.0)*log(2.0*M_PI)) - ((dnrow/2.0)*log(var)) - (sum1/(2.0*var)) ); + + + return lik; + +} + + +double null_lik(SIMX *trueX, double *y, double var, double mu, int nrow, int nbin) { + + int i,j,m; + double lik,sum1,Nsamp; + + Nsamp = 0.0; + sum1 = 0.0; + for (i=0; i= nbin ) { + sum1 = sum1 + (y[i] - mu)*(y[i] - mu); + Nsamp = Nsamp + 1.0; + } + } + + lik = ( ((-Nsamp/2.0)*log(2.0*M_PI)) - ((Nsamp/2.0)*log(var)) - (sum1/(2.0*var)) ); + + + return lik; + +} + + +SIMX* drawX(XMAT *xmat, int ncol, int nrow, long *idum) { + + int i,j,m; + double check,dncol,max,muNi,ran,sdNi,varNi; + + int *Xvec=NULL, *Ni=NULL; + double *Xprob=NULL; + + SIMX *trueX=NULL; + + + trueX = (SIMX*)calloc(1,sizeof(SIMX)); + + Xvec = (int*)calloc(nrow,sizeof(int)); + Xprob = (double*)calloc(nrow,sizeof(double)); + + Ni = (int*)calloc(ncol,sizeof(int)); + + + dncol = (double)ncol; + ncol = (int)ncol; + + + for (i=0; i (*xmat).cumX[i][m-1] ) { + m++; + } + + if (m > ncol) { + Rprintf("hbrem drawX ERROR m = %i, ran = %f\n",m,ran); + return(NULL); + } + + Xvec[i] = m; + Xprob[i] = (*xmat).X[i][m-1]; + + Ni[Xvec[i] - 1] = Ni[Xvec[i] - 1] + 1; + } + + + muNi = 0.0; + for (i=0; i= nbin ) { + yTbar[m] = yTbar[m] + y[i]; + sum1 = sum1 + (y[i]*y[i]); + sum6 = sum6 + y[i]; + Nsamp = Nsamp + 1.0; + } + } + sum6 = sum6/Nsamp; + + Ngen = 0.0; + for (j=0; j= nbin ) { + yTbar[j] = yTbar[j]/((double)(*trueX).Ni[j]); + Ngen = Ngen + 1.0; + } + } + + + kT = 0.0; + inc = 1.0/dprc; + stmax = -1000000.0; + + for (i=0; i= nbin ) { + + dn = (double)(*trueX).Ni[j]; + deno = ( 1.0 - kT + (kT*dn) ); + + sum2 = sum2 + (dn/deno); + prod1 = prod1 + log(deno); + sum3 = sum3 + (dn*dn*(yTbar[j])*(yTbar[j]))/deno; + sum4 = sum4 + (dn*(yTbar[j]))/deno; + } + } + + + inter = ( sum1 - (kT*sum3) - ((1.0-kT)*(sum4*sum4/sum2)) ); + + + sum2 = -log(sum2)/2.0; + prod1 = -prod1/2.0; + + pow1 = ( Ngen - 1.0 )/2.0; + pow2 = ( Nsamp - 1.0 )/2.0; + + kT_pdf[i] = ( (pow1*log(1-kT)) + sum2 + prod1 - (pow2*log(inter)) ); + + if (kT_pdf[i] > stmax) { + stmax = kT_pdf[i]; + } + + kT = kT + inc; + + } + + sum5 = 0.0; + stmin = (stmax - 703.0); + + for (i=0; i (*kTdist).cdf[i] ) { + i++; + } + if ( i > (prc + 1) ) { + Rprintf("error in draw of kT\n"); + return(0.0); + } + + kT = ((double)i)/dprc; + + + return kT; + +} + + +double draw_knownvar(GRKT *kTdist, int *Ni, int ncol, double kT, double nu, int nbin) { + + int j; + + double chi,scale,var; + double deno,dn,sum2,sum3,sum4; + + + if ( kT == 1.0 ) { + + var = 0.0; + + } else { + + chi = rchisq(nu); + + sum2 = 0.0; + sum3 = 0.0; + sum4 = 0.0; + for (j=0; j= nbin ) { + + dn = (double)Ni[j]; + deno = ( 1.0 - kT + kT*dn ); + + sum2 = sum2 + (dn/deno); + sum3 = sum3 + (dn*dn*((*kTdist).yTbar[j])*((*kTdist).yTbar[j]))/deno; + sum4 = sum4 + (dn*((*kTdist).yTbar[j]))/deno; + } + } + + + // scale = rho squared multiplied by nu=degrees of freedom + scale = ( ((*kTdist).sum_ysq)/(1.0 - kT) - (kT/(1.0 - kT))*sum3 - (sum4*sum4/sum2) ); + var = scale/chi; + + } + + return var; + +} + + +double draw_nullvar(SIMX *selX, double *y, int nrow, int nbin) { + + int j,m; + + double chi,nu,scale,var; + double dnrow,sum2,sum3,Nsamp; + + + dnrow = (double)nrow; + nrow = (int)nrow; + + Nsamp = 0.0; + sum2 = 0.0; + sum3 = 0.0; + for (j=0; j= nbin ) { + sum2 = sum2 + (y[j]*y[j]); + sum3 = sum3 + y[j]; + Nsamp = Nsamp + 1.0; + } + } + sum3 = sum3/Nsamp; + + + nu = (Nsamp - 1.0); + chi = rchisq(nu); + + // scale = rho squared multiplied by nu=degrees of freedom + scale = ( sum2 - (Nsamp*sum3*sum3) ); + var = scale/chi; + + + return var; + +} + + +double draw_knownmu(GRKT *kTdist, int *Ni, int ncol, double kT, double var, int nbin) { + + int j; + + double mu,nmean,nsigmasq,sigma; + double deno,dn,sum2,sum4; + + + sum2 = 0.0; + sum4 = 0.0; + for (j=0; j= nbin ) { + + dn = (double)Ni[j]; + deno = ( 1.0 - kT + kT*dn ); + + sum2 = sum2 + (dn/deno); + sum4 = sum4 + (dn*((*kTdist).yTbar[j]))/deno; + } + } + + nmean = sum4/sum2; + nsigmasq = var/sum2; + + sigma = sqrt(nsigmasq); + mu = rnorm(0.0, sigma); + mu = mu + nmean; + + return mu; + +} + + +double draw_nullmu(SIMX *selX, double *y, int nrow, double var, int nbin) { + + int j,m; + + double mu,nmean,nsigmasq,sigma; + double dnrow,sum2,Nsamp; + + + dnrow = (double)nrow; + nrow = (int)nrow; + + + sum2 = 0.0; + Nsamp = 0.0; + for (j=0; j= nbin ) { + sum2 = sum2 + y[j]; + Nsamp = Nsamp + 1.0; + } + } + sum2 = sum2/Nsamp; + + + nmean = sum2; + nsigmasq = var/Nsamp; + + sigma = sqrt(nsigmasq); + mu = rnorm(0.0, sigma); + mu = mu + nmean; + + return mu; + +} + + +double draw_knownTi(GRKT *kTdist, int *Ni, double kT, double var, double mu, int nbin, int index) { + + double dn,deno; + double nmean,nsigmasq,sigma,ti; + + + ti = 0.0; + if ( (int)Ni[index] >= nbin ) { + + dn = (double)Ni[index]; + deno = ( 1.0 - kT + kT*dn ); + + nmean = (kT*dn*((*kTdist).yTbar[index] - mu))/deno; + nsigmasq = (kT*(1.0 - kT)*var)/deno; + + sigma = sqrt(nsigmasq); + ti = rnorm(0.0, sigma); + ti = ti + nmean; + + } + + return ti; + +} +