inst/doc/metaGE-vignette.R

In metaGE: Meta-Analysis for Detecting Genotype x Environment Associations

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = '70%'
)

## ----intall metaGE------------------------------------------------------------
if (!require('metaGE')){
  install.packages('metaGE')
} 
library('metaGE')

## ----setup,warning=FALSE,message=FALSE----------------------------------------

library(tidyverse)
library(DT)

## For graphical displays
library(data.table)
library(corrplot)
library(ggplot2)


## ----listing files of GWAS results--------------------------------------------
## Get the folder containing the association file
RepData <- system.file("extdata", package = "metaGE")

## Get the complete list of association files
File.list <- list.files(RepData ,full.names = TRUE) %>% 
  tibble(Names = .)%>% 
  mutate(ShortNames = Names %>%
           str_remove(pattern = paste0(RepData,"/")) %>%
           str_remove(pattern = "_3DF.txt"))  %>%
  select(ShortNames,Names) %>% 
  deframe

File.list

## ----looking at single file---------------------------------------------------
## Have a look at the first one
fread(File.list[1]) %>% head()  %>% datatable()

## ----metaGE collect-----------------------------------------------------------
###Build the dataset
## First provide the list of variable names 
Names.list <- list(MARKER="Marker_Name",
                   CHR="Chromosome",
                   POS="Marker_Position",
                   FREQ="Maf",
                   EFFECT="SNP_Weight",
                   PVAL="Pvalue",
                   ALLELE0="Allele1",
                   ALLELE1="Allele2")

MinFreq <- 0.07

## Now collect
MetaData <- metaGE.collect(FileNames = File.list, VariableNames = Names.list, MinFreq = MinFreq)
head(MetaData$Data) %>% datatable()


## ----metaGE collect folders---------------------------------------------------
## Get the list of the association files
File.list1 <- File.list[1:2]

## Get the list of the other association files
File.list2 <- File.list[3]

File.listc <- list("List1" = File.list1, "List2" = File.list2)
File.listc

###Build the dataset
## First provide the list of variable names 
Names.listc <- list(Names.list, Names.list)

MinFreq <- 0.07

## Now collect
MetaData <- metaGE.collect(FileNames = File.listc, VariableNames = Names.listc, MinFreq = MinFreq)

head(MetaData$Data)  %>% datatable()

## ----import data--------------------------------------------------------------
# Import the data
data("metaData")

## ----build matcorr------------------------------------------------------------
Threshold <- 0.8
MatCorr <- metaGE.cor(metaData, Threshold = Threshold)

## ----show matcorr,fig.align = 'center',fig.height=8,fig.width=8---------------
corrplot(MatCorr,order = "hclust")

## ----metaGE fit---------------------------------------------------------------
## Fixed Effect
FeDF <- metaGE.fit(metaData, MatCorr, Method = "Fe")
head(FeDF %>% select(CHR, POS, MARKER, Mu, Tau, PVALUE))  %>% datatable()

## Random Effect
ReDF <- metaGE.fit(metaData, MatCorr, Method = "Re")
head(ReDF %>% select(CHR, POS, MARKER, Mu, Tau, PVALUE))  %>% datatable()

## ----metaGE pvalplot,fig.height=3,fig.width=5,out.width='45%'-----------------
Pvalue.list <- list('PVALUE.Fe'= FeDF$PVALUE, 'PVALUE.Re'= ReDF$PVALUE)
plott <- map(names(Pvalue.list),~metaGE.pvalplot(Pvalue.list[[.x]], Main= .x) )

## ----correcting pvalues-------------------------------------------------------
CorrectingPValues <- function(Pvalues){

  ## Get a p0 estimate
  p0 = min(2*sum(Pvalues > 0.5)/length(Pvalues),1-1/length(Pvalues))

  ## Get the corrected p-values
  CorrectedPval <- stats::p.adjust(Pvalues, method = 'BH')*p0

  return(CorrectedPval)
}

## ----FDR control--------------------------------------------------------------
## FDR control
Alpha <- 0.05
Signif <- lapply(Pvalue.list, FUN = function(i){
                  return(CorrectingPValues(i) %>% 
                    `<`(Alpha) %>% 
                    which)})

lapply(X = Signif,FUN = length)

## ----manhattan plot,fig.align='center',fig.height=6,fig.width=10--------------

PvalThresholdFe <-FeDF[Signif$PVALUE.Fe,]$PVALUE%>% max %>% max(.,0)
manhattanFe <- metaGE.manhattan(Data = FeDF,VarName = 'PVALUE', Threshold = PvalThresholdFe,Main = '-log10(Pval) alongside the chromosome Fe method' )
print(manhattanFe)

PvalThresholdRe <- ReDF[Signif$PVALUE.Re,]$PVALUE%>% max %>% max(.,0)
manhattanRe <- metaGE.manhattan(Data = ReDF,VarName = 'PVALUE', Threshold = PvalThresholdRe,Main = '-log10(Pval) alongside the chromosome Re method')
print(manhattanRe)


## ----heatmap,fig.align='center',fig.height=6,fig.width=10---------------------
heatmapFe <- metaGE.heatmap(Data = FeDF[Signif$PVALUE.Fe,],Prefix = "Z.",Main = "Z-scores of Fe significant markers across environments")


heatmapRe <- metaGE.heatmap(Data = ReDF[Signif$PVALUE.Re,],Prefix = "Z.", Main = "Z-scores of Re significant markers across environments")


## ----local score FE-----------------------------------------------------------
x <- 3
FeDF_ls <- metaGE.lscore(Data = FeDF, PvalName = "PVALUE", xi = x)

## ----sigzones Fe--------------------------------------------------------------
FeDF_ls$SigZones %>% datatable()

## ----manhattan Fe lscore,fig.align='center',fig.height=6,fig.width=10---------
manhattanFe_lscore <- metaGE.manhattan(Data = FeDF_ls$Data,
                              VarName = "SCORE",
                              Score = T,
                              SigZones = FeDF_ls$SigZones )
print(manhattanFe_lscore+ggtitle('Score local alongside the chromosome Fe method'))


## ----import environment data--------------------------------------------------
data("envDesc")
envDesc %>% datatable()

## ----incidence matrix---------------------------------------------------------
## Build the incidence matrix 
(Incidence.Temp <- metaGE.incidence(VarName = "Temp",Covariate = envDesc,EnvName = "ShortName", Data = metaData))

## ----metaGE  contrast test----------------------------------------------------
## Build the list of Incidence
Incidence.list <- list(Temp = Incidence.Temp,
                       Diff.Temp = Incidence.Temp)

#Build the list of Contrast
Contrast.list <- list(Temp = NULL,
                      Diff.Temp = matrix(c(1,-1,0,0,1,-1),2,byrow = T)) 


ContrastDF <- metaGE.test(Data = metaData, MatCorr = MatCorr,
                          Incidence = Incidence.list,
                          Contrast = Contrast.list)

## ----FDR control test contrast------------------------------------------------
Alpha <- 0.05
SignifContrast <- apply(ContrastDF %>% select(contains("PVALUE.")),2,
                        FUN = function(i){return(CorrectingPValues(i) %>% 
                                                    `<`(Alpha) %>% 
                                                      which)})

lapply(SignifContrast,length)


## ----contrast heatmap,fig.align='center',fig.height=6,fig.width=10------------
# Specify the groups of environment
heatmap_Temp <- metaGE.heatmap(Data = ContrastDF[SignifContrast$PVALUE.Temp,],EnvGroups = envDesc[,1:2],Prefix = "Z.",Main = "Z-scores of markers with contrasted effect according the temperature")



## ----regression test----------------------------------------------------------
RegressionDF <- metaGE.test(Data=metaData,MatCorr = MatCorr,Covariate = envDesc[,c(1,5,6)],EnvName = "ShortName")

## ----FDR control test regression----------------------------------------------
SignifRegression <- apply(RegressionDF %>% select(contains("PVALUE.")),2,
                        FUN = function(i){return(CorrectingPValues(i) %>% 
                                                    `<`(Alpha) %>% 
                                                      which)})

lapply(SignifRegression,length)


## ----significant marker test regression---------------------------------------
RegressionDF[SignifRegression$PVALUE.Tnight.mean,] %>% select(CHR, POS, MARKER, PVALUE.Tnight.mean) %>% arrange(PVALUE.Tnight.mean) %>% head()  %>% datatable()

## ----metaGE regplot,fig.align='center',fig.height=7,fig.width=10--------------
metaGE.regplot(Data = metaData, Covariate = envDesc,VarName = "Tnight.mean",EnvName = "ShortName", MarkerName = "AX-90548528", aesCol = "Classification")