## ----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")