R/test_segregation.R

In onemap: Construction of Genetic Maps in Experimental Crosses

#######################################################################
##                                                                     ##
## Package: onemap                                                     ##
##                                                                     ##
## File: test_segregation.R                                            ##
## Contains: test_segregation_of_a_marker,                             ##
## test_segregation, plot.onemap_segreg_test,                          ##
## print.onemap_segreg_test, Bonferroni_alpha, select_segreg           ##
##                                                                     ##
## Written by Antonio Augusto Franco Garcia with minor modifications   ##
## by Marcelo Mollinari and Cristiane Taniguti                         ##
## copyright (c) 2015 Antonio Augusto Franco Garcia                    ##
##                                                                     ##
## First version: 2015/04/18                                           ##
## License: GNU General Public License version 3 or later              ##
##                                                                     ##
#######################################################################

globalVariables(c("Marker", "p.value"))

##' test_segregation_of_a_marker
##'
##' Applies the chi-square test to check if markers are following the
##' expected segregation pattern, i. e., 1:1:1:1 (A), 1:2:1 (B), 3:1 (C) and 1:1 (D)
##' according to OneMap's notation. It does not use Yate's correction.
##'
##' First, the function selects the correct segregation pattern, then it
##' defines the H0 hypothesis, and then tests it, together with percentage of
##' missing data.
##' 
##' @importFrom stats chisq.test cor.test dist na.omit qchisq rf setNames
##'
##' @param x an object of class \code{onemap}, with data and additional information.
##' @param marker the marker which will be tested for its segregation.
##' @param simulate.p.value a logical indicating whether to compute p-values by Monte Carlo simulation.
##'
##' @return a list with the H0 hypothesis being tested, the chi-square statistics,
##' the associated p-values, and the \% of individuals genotyped.
##'
##' @examples
##' 
##' data(onemap_example_bc) # Loads a fake backcross dataset installed with onemap
##' test_segregation_of_a_marker(onemap_example_bc,1)
##'
##' data(onemap_example_out) # Loads a fake outcross dataset installed with onemap
##' test_segregation_of_a_marker(onemap_example_out,1)
##' 
##' @export
test_segregation_of_a_marker <- function(x, marker, simulate.p.value=FALSE) {
  # Break if marker has any genotype information
  if(all(x$geno[,marker] ==0)) {
    warning(cat("Marker", marker, "do not have genotype information for any sample. We suggest to remove it from the dataset\n"))
    return(list(Hypothesis=NA, qui.quad=NA, p.val=NA, perc.genot=NA))
  } else {
    ## Segregation pattern for each marker type
    p.a <- rep(1/4, 4); p.b <- c(1/4, 1/2, 1/4); p.c <- c(3/4, 1/4); p.d <- rep(1/2, 2)
    ## Counting each category, removing missing data (coded as 0)
    count <- table(x$geno[,marker], exclude=0)
    ## Do the chisq test, using the appropriate expected segregation
    ## grepl() allows finding the marker type (it has the letter in the argument)
    ## Some markers may not have data for some classes, so fill them with 0
    if (grepl("A.H.B",x$segr.type[marker])) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      if (is.element(3,x$geno[,marker])) c3 <- count[names(count)==3] else c3 <- 0
      qui <- chisq.test(c(c1,c2,c3), p=p.b, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:2:1"
    }
    else if (grepl("C.A",x$segr.type[marker])) {
      if (is.element(2,x$geno[,marker])) c1 <- count[names(count)==2] else c1 <- 0
      if (is.element(1,x$geno[,marker])) c2 <- count[names(count)==1] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=rev(p.c), correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "3:1"
    }
    else if (grepl("D.B",x$segr.type[marker])) {
      if (is.element(2,x$geno[,marker])) c1 <- count[names(count)==2] else c1 <- 0
      if (is.element(1,x$geno[,marker])) c2 <- count[names(count)==1] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=rev(p.c), correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "3:1"
    }
    else if (grepl("A.H",x$segr.type[marker])) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=p.d, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:1"
    }
    else if (grepl("A.B",x$segr.type[marker])) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(3,x$geno[,marker])) c2 <- count[names(count)==3] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=p.d, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:1"
    }
    else if (grepl("A",x$segr.type[marker]) & inherits(x,"outcross")) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      if (is.element(3,x$geno[,marker])) c3 <- count[names(count)==3] else c3 <- 0
      if (is.element(4,x$geno[,marker])) c4 <- count[names(count)==4] else c4 <- 0
      qui <- chisq.test(as.vector(c(c1,c2,c3,c4)), p=p.a, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:1:1:1"
    }
    else if (grepl("B",x$segr.type[marker]) & inherits(x,"outcross")) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      if (is.element(3,x$geno[,marker])) c3 <- count[names(count)==3] else c3 <- 0
      qui <- chisq.test(as.vector(c(c1,c2,c3)), p=p.b, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:2:1"
    }
    else if (grepl("C",x$segr.type[marker]) & inherits(x,"outcross")) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=p.c, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "3:1"
    }
    else if (grepl("D",x$segr.type[marker]) & inherits(x,"outcross")) {
      if (is.element(1,x$geno[,marker])) c1 <- count[names(count)==1] else c1 <- 0
      if (is.element(2,x$geno[,marker])) c2 <- count[names(count)==2] else c2 <- 0
      qui <- chisq.test(as.vector(c(c1,c2)), p=p.d, correct = FALSE, simulate.p.value = simulate.p.value)
      H0 <- "1:1"
    }
    #impossible to test: dominant and co-dominant mixed in the same marker
    #however, it will not work with "qui <- NA"; using NULL instead
    else if (grepl("M.X",x$segr.type[marker])) {
      qui <- NULL
      qui$statistic <- NA
      qui$p.value <- NA
      H0 <- "Check data!"
    }
    
    return(list(Hypothesis=H0, qui.quad=qui$statistic, p.val=qui$p.value,
                perc.genot=100*(sum(table(x$geno[,marker], exclude=0))/x$n.ind)))
  }
}
##'

##' test_segregation
##'
##' Using OneMap internal function test_segregation_of_a_marker(),
##' performs the Chi-square test to check if all markers in a dataset are following
##' the expected segregation pattern, i. e., 1:1:1:1 (A), 1:2:1 (B), 3:1 (C) and 1:1 (D)
##' according to OneMap's notation.
##'
##' First, it identifies the correct segregation pattern and corresponding H0 hypothesis,
##' and then tests it.
##'
##' @param x an object of class \code{onemap}, with data and additional information.
##' @param simulate.p.value a logical indicating whether to compute p-values by Monte Carlo simulation.
##' 
##' @return an object of class onemap_segreg_test, which is a list with marker name,
##' H0 hypothesis being tested, the chi-square statistics, the associated p-values
##' and the \% of individuals genotyped. To see the object, it is necessary to print
##' it.
##'
##' @examples
##' \donttest{
##'  data(onemap_example_out) # Loads a fake outcross dataset installed with onemap
##'  Chi <- test_segregation(onemap_example_out) # Performs the chi-square test for all markers
##'  print(Chi) # Shows the results
##' }
##' 
##' @export
test_segregation <- function(x, simulate.p.value = FALSE) {
    if (inherits(x,"onemap")) {
        y <- list(Marker=dimnames(x$geno)[[2]],
                  Results.of.tests=sapply(1:x$n.mar, function(onemap.object, marker)
                      test_segregation_of_a_marker(onemap.object, marker, simulate.p.value = simulate.p.value),
                      onemap.object=x))
                                        # sapply iterates from 1 to x$n.mar; x is fixed (onemap object with data)
        class(y) <- c("onemap_segreg_test")
        invisible(y) #returns y without showing it
    }
    else stop("This is not an OneMap object with raw data")
}
##'

##' Show the results of segregation tests
##'
##' It shows the results of Chisquare tests performed for all markers in a onemap object
##' of cross type outcross, backcross, F2 intercross or recombinant inbred lines.
##'
##' @param x an object of class onemap_segreg_test
##'
##' @param ... currently ignored
##'
##' @return a dataframe with marker name, H0 hypothesis, chi-square statistics,
##' p-values, and % of individuals genotyped.
##'
##' @examples
##' 
##' \donttest{
##'  data(onemap_example_out) # Loads a fake outcross dataset installed with onemap
##'  Chi <- test_segregation(onemap_example_out) # Performs the chi-square test for all markers
##'  print(Chi) # Shows the results
##' }
##' 
##' @method print onemap_segreg_test
##' @export
print.onemap_segreg_test <- function(x,...) {
    Z <- data.frame(Marker=x$Marker,
                    H0=unlist(x$Results.of.tests[1,]),
                    Chi.square=unlist(x$Results.of.tests[2,]),
                    p.value=unlist(x$Results.of.tests[3,]),
                    Perc.genot=round(unlist(x$Results.of.tests[4,]),2))
    colnames(Z) <- c("Marker","H0","Chi-square","p-value","% genot.")
    return(Z)
}
##'

##' Plot p-values for chi-square tests of expected segregation
##'
##' Draw a graphic showing the p-values (re-scaled to -log10(p-values)) associated with the
##' chi-square tests for the expected segregation patterns for all markers in a dataset.
##' It includes a vertical line showing the threshold for declaring statistical significance
##' if Bonferroni's correction is considered, as well as the percentage of markers that
##' will be discarded if this criterion is used.
##'
##' @param x an object of class onemap_segreg_test (produced by onemap's function
##' test_segregation()), i. e., after performing segregation tests
##' @param order a variable to define if p-values will be ordered in the plot
##'
##' @param ... currently ignored
##'
##' @return a ggplot graphic
##'
##' @import ggplot2
##'
##' @examples
##' \donttest{
##'  data(onemap_example_bc) # load OneMap's fake dataset for a backcross population
##'  BC.seg <- test_segregation(onemap_example_bc) # Applies chi-square tests
##'  print(BC.seg) # Shows the results
##'  plot(BC.seg) # Plot the graph, ordering the p-values
##'  plot(BC.seg, order=FALSE) # Plot the graph showing the results keeping the order in the dataset
##'
##'  data(onemap_example_out) # load OneMap's fake dataset for an outcrossing population
##'  Out.seg <- test_segregation(onemap_example_out) # Applies chi-square tests
##'  print(Out.seg) # Shows the results
##'  plot(Out.seg) # Plot the graph, ordering the p-values
##'  plot(Out.seg, order=FALSE) # Plot the graph showing the results keeping the order in the dataset
##' }
##' @method plot onemap_segreg_test
##' @export
plot.onemap_segreg_test <- function(x, order=TRUE,...) {
                                        # Create a data frame
    Z <- data.frame(Marker=x$Marker,
                    X.square=unlist(x$Results.of.tests[2,]),
                    p.value=unlist(x$Results.of.tests[3,]))
    Bonf <- -log10(.05/nrow(Z)) #Bonferroni's threshold'
    Z$signif <- factor(ifelse(-log10(Z$p.value)<Bonf,"non sign.","sign."))
    Z$order <- 1:nrow(Z)
                                        # % of distorted
    perc <- 100*(1-(table(Z$signif)[1]/nrow(Z)))
                                        # Keeping markers in their original order (not alphanumeric), or by p-values (default)
    if (order!=TRUE) Z$Marker <- factor(Z$Marker, levels = Z$Marker[order(Z$order)])
    else Z$Marker <- factor(Z$Marker, levels = Z$Marker[order(Z$p.value, decreasing=TRUE)])
                                        # Plotting
    g <- ggplot(data=Z, aes(x=Marker, y=-log10(p.value)))
    g <- g + ylab(expression(-log[10](p-value)))
    g <- g + geom_point(aes(color=signif), stat="identity",size=2.5)
    g <- g + scale_colour_manual(name=paste("Bonferroni\n","(",round(perc,0),"% distorted)",sep=""),
                                 values = c("#46ACC8","#B40F20"))
    g <- g + geom_hline(yintercept = Bonf, colour="#E58601", linetype = "longdash")
    g <- g + coord_flip()
    if (nrow(Z)>30) g <- g + theme(axis.text.y = element_blank())
    g
}
##'


##' Calculates individual significance level to be used to achieve a global alpha (with Bonferroni)
##'
##' It shows the alpha value to be used in each chi-square segregation test, in order to achieve
##' a given global type I error. To do so, it uses Bonferroni's criteria.
##'
##'@importFrom methods is
##'
##' @param x an object of class onemap_segreg_test
##' @param global.alpha the global alpha that
##'
##' @return the alpha value for each test (numeric)
##'
##' @examples
##' \donttest{
##'  data(onemap_example_bc) # Loads a fake backcross dataset installed with onemap
##'  Chi <- test_segregation(onemap_example_bc) # Performs the chi-square test for all markers
##'  print(Chi) # Shows the results of the Chi-square tests
##'  Bonferroni_alpha (Chi) # Shows the individual alpha level to be used
##' }
##' 
##' @export
Bonferroni_alpha <- function(x, global.alpha=0.05) {
    if (!inherits(x,"onemap_segreg_test")) stop("This is not an object of class onemap_segreg_test")
    alpha.Bonf <- global.alpha/length(x$Marker)
    return(alpha.Bonf)
}
##'

##' Show markers with/without segregation distortion
##'
##' A function to shows which marker have segregation distortion if Bonferroni's correction is
##' applied for the Chi-square tests of mendelian segregation.
##'
##' @param x an object of class onemap_segreg_test
##' @param distorted a TRUE/FALSE variable to show distorted or non-distorted markers
##' @param numbers a TRUE/FALSE variable to show the numbers or the names of the markers
##' @param threshold a number between 0 and 1 to specify the threshold (alpha) to be considered in the test. 
##' If NULL, it uses the threshold alpha = 0.05. Bonferroni correction is applied for multiple test correction. 
##'
##' @return a vector with marker names or numbers, according to the option for "distorted" and "numbers"
##'
##' @examples
##' \donttest{
##'  # Loads a fake backcross dataset installed with onemap
##'  data(onemap_example_out)
##'  # Performs the chi-square test for all markers
##'  Chi <- test_segregation(onemap_example_out)
##'  # To show non-distorted markers
##'  select_segreg(Chi)
##'  # To show markers with segregation distortion
##'  select_segreg(Chi, distorted=TRUE)
##'  # To show the numbers of the markers with segregation distortion
##'  select_segreg(Chi, distorted=TRUE, numbers=TRUE)
##' }
##' @export
select_segreg <- function(x, distorted=FALSE, numbers=FALSE, threshold = NULL) {
    if (!inherits(x,"onemap_segreg_test")) stop("This is not an object of class onemap_segreg_test")
    Z <- data.frame(Marker=x$Marker,
                    p.value=unlist(x$Results.of.tests[3,]))
    if(is.null(threshold)) thr <- Bonferroni_alpha(x, global.alpha = 0.05)
    else  thr <- Bonferroni_alpha(x, global.alpha = threshold)
    if (distorted==FALSE) Z <- subset(Z, p.value>=thr)
    else Z <- subset(Z, p.value<thr)
    if (numbers==TRUE) return(which(x$Marker %in% as.vector(Z[,1])))
    else return(as.vector(Z[,1]))
}
##'