inst/doc/gwas-manual.R

In mlmm.gwas: Pipeline for GWAS Using MLMM

## ----setup, include = FALSE----------------------------------------------
    knitr::opts_chunk$set(
      collapse = TRUE,
      comment = "#>"
    )

## ----data----------------------------------------------------------------
    library(mlmm.gwas)
    data("mlmm.gwas.AD")
    ls()

## ----matrixFormat, include=FALSE-----------------------------------------
    # Change le formatage par d??faut des matrices
    # On le met ici pour pas qu'il soit visible quand on fait ls() ci dessus
    knit_print.matrix <- function(x, ...){
        knitr::knit_print(as.data.frame(x), ...)
    }

## ----floweringDateAD-----------------------------------------------------
    head(floweringDateAD)

## ----X-------------------------------------------------------------------
    Xa[1:5,1:5]

## ----K-------------------------------------------------------------------
    ##You can get K from X with the following command line
    ##We are not doing it here because the X matrix included in data("mlmm.gwas.AD")
    ##is a subset of a larger matrix.
    #K.add = Xa %*% t(Xa)

    K.add[1:5,1:5]

## ----mlmm----------------------------------------------------------------
    res_mlmm = mlmm_allmodels(floweringDateAD, list(Xa), list(K.add))

## ----manhattan, fig.height=5, fig.width=5--------------------------------
    manhattan.plot( res_mlmm )

## ----eBIC----------------------------------------------------------------
    sel_XX = frommlmm_toebic(list(Xa), res_mlmm)
    res_eBIC = eBIC_allmodels(floweringDateAD, sel_XX, list(K.add), ncol(Xa))
    
    res_eBIC

## ----threshold-----------------------------------------------------------
    res_threshold <- threshold_allmodels(threshold=NULL, res_mlmm)

## ----estimation----------------------------------------------------------
    sel_XXclass = fromeBICtoEstimation(sel_XX, res_eBIC, res_threshold)
    effects = Estimation_allmodels(floweringDateAD, sel_XXclass, list(K.add))
    
    effects

## ----boxplot-------------------------------------------------------------
    genotypes.boxplot(Xa, floweringDateAD, "SNP303", effects)

## ----completeExample, eval=FALSE-----------------------------------------
#      library(mlmm.gwas)
#      data("mlmm.gwas.AD")
#  
#      res_mlmm = mlmm_allmodels(floweringDateAD, list(Xa), list(K.add))
#  
#      manhattan.plot( res_mlmm )
#  
#      sel_XX = frommlmm_toebic(list(Xa), res_mlmm)
#      res_eBIC = eBIC_allmodels(floweringDateAD, sel_XX, list(K.add), ncol(Xa))
#  
#      res_threshold <- threshold_allmodels(threshold=NULL, res_mlmm)
#  
#      sel_XXclass = fromeBICtoEstimation(sel_XX, res_eBIC, res_threshold)
#      effects = Estimation_allmodels(floweringDateAD, sel_XXclass, list(K.add))
#  
#      genotypes.boxplot(Xa, floweringDateAD, "SNP303", effects)

## ----runEntirePipeline---------------------------------------------------
    results = run_entire_gwas_pipeline(floweringDateAD, list(Xa), list(K.add), threshold=NULL)
    names(results)