read.ped <- function(file, marker.names=NULL)
{
    raw.data <- read.table(file, header=FALSE)

    subject.data <- raw.data[,1:6]
    allele.mat   <- as.matrix(raw.data[,7:ncol(raw.data)])
    num.markers  <- ncol(allele.mat)/2
    geno.mat     <- matrix(ncol=num.markers, nrow=nrow(subject.data))

    for (i in 1:num.markers)
    {
        a <- 2*(i-1) + 1
        geno.mat[,i] <- allele.mat[,a] + allele.mat[,a+1] - 1
    }

    data <- cbind(subject.data, as.data.frame(geno.mat))

    if (is.null(marker.names))
    {
        marker.names <- paste("m", sep="", 1:num.markers)
    }
    else if (num.markers!=length(marker.names))
    {
        stop("Number of marker names supplied to method (",
                length(marker.names),
                ") different from the number of marker columns in ",
                "ped file (", num.markers, ")\n")
    }
    colnames(data) <- c("Family", "SUBJECT.NAME", "Father", "Mother",
            "Sex", "Phenotype", as.character(marker.names))
    data
}

dominance.part <- function(g)
{
    1 - g^2
}

include.dominance <- function(g)
{
    as.factor(g)
}


scan.markers <- function(formula, h0=NULL, data, marker.data, ...)
{
    if (!is.character(formula))
    {
        formula <- deparse(substitute(formula))
    }
    if (!is.null(h0))
    {
        if (!is.character(h0))
        {
            h0 <- deparse(substitute(h0))
        }
    }
    else
    {
        h0 <- sub("~ .*", "~ 1", formula)
    }

    results <- NULL
    fit0 <- lm(as.formula(h0), data=data)
    for (i in 1:nrow(marker.data))
    {
        marker      <- as.character(marker.data$name[i])
        data$MARKER <- data[,marker]
        fit         <- lm(as.formula(formula), data=data)
        an          <- anova(fit0, fit)

        pval    <- an$"Pr(>F)"[2]
        result  <- data.frame(
                chrom   = I(as.character(marker.data$chrom[i])),
                cM      = marker.data$cM[i],
                pval    = pval,
                logP    = -log10(pval))
        results <- rbind(results, result)
    }
    rownames(results) <- as.character(marker.data$name)
    results
}

permute.scan.markers <- function(formula, h0=NULL, data, marker.data, perms=200, ...)
{
    if (!is.character(formula))
    {
        formula <- deparse(substitute(formula))
    }
    phenotype.name <- sub("\\s*~ .*", "", formula)
    phenotype <- data[,phenotype.name]

    max.logPs <- numeric(perms)
    cat("Permuting...\n")
    for (i in 1:perms)
    {
        cat(paste("[",i,"]",sep=""))
        data[,phenotype.name] <- sample(phenotype)
        results <- scan.markers(formula, h0=h0, data=data,
                marker.data=marker.data)
        max.logPs[i] <- max(na.rm=T, results$logP)
    }
    cat("...Done\n")
    max.logPs
}


plot.chrom.scan <- function(data,
        add=FALSE,
        threshold=NULL,
        col="black",
        max.threshold.height=0.25, ...)
{
    if (!add)
    {
        plot(range(data$cM),
                c(0, max(c(data$logP,threshold/max.threshold.height))),
                xlab="cM", ylab="logP", type="n", las=1, ...)
    }
    if (!is.null(threshold))
    {
        abline(h=threshold, lty=2)
    }
    lines(data$cM, data$logP, col=col, type="b")
}