DeepGenomeScan : A Deep Learning Approach for Whole Genome Scan (WGS) and Genome-wide Association Studies (GWAS)

This package implements the genome scan and genome-wide association studies using deep neural networks (i.e, Multi-Layer Perceptron (MLP), Convolutional Neural Network (CNN)). DeepGenomeScan offers heuristic computational framework integrating different neural network architectures (i.e.,Multi-Layer Perceptron (MLP), convolutional neural network(CNN)) and robust resampling methods, as well as the Model-Agnostic interpretation of feature importance for convolutional neural networks. DeepGenomeScan, in other words, deep learning for genome-wide scanning, is a deep learning approach for detecting signatures of natural selection and for performing various omics-based genome-wide association studies, such as GWAS, PWAS, TWAS, MWAS. The design makes the implemention user-friendly. It is compatible with most self-defined machine learning models (the self-defined models shuold be complete, including tunable parameters, fitted model, predicted model, examples can be found in our tutorial). Users can adopt the package's framework to study various ecological and evolutionary questions.

Install packages

library("devtools")

devtools::install_github("xinghuq/DeepGenomeScan")
devtools::install_github("xinghuq/CaretPlus/pkg/caret")

Dependencies and environment requirements

Note: Environment requirements: python should be installed and the python package of Keras and Tensorflow should also be installed and work properly with the system

Checking the python environment

library("rappdirs")
library("reticulate")
reticulate::use_python("/usr/bin/python3")
library(caret) ### for ML calling functions and performance estimation, users should use the modified version at xinghuq/CaretPlus/caret instead of the original version
library(keras)  
library("tensorflow")

checking if Tensorflow works properly
K0=keras::backend()

requireNamespace("KLFDAPC")

 if (!requireNamespace("KLFDAPC", quietly=TRUE))

  devtools::install_github("xinghuq/KLFDAPC")
  
 if (!requireNamespace("DA", quietly=TRUE))
 
  devtools::install_github("xinghuq/DA")

if (!requireNamespace("keras", quietly=TRUE))

 install.packages("keras")
  
 if (!requireNamespace("tensorflow", quietly=TRUE))
 
  install.packages("tensorflow")
 if (!requireNamespace("kerasR", quietly=TRUE))
 
  install.packages("kerasR")

library

library(DeepGenomeScan)
library(caret)### for ML calling functions and performance estimation
library(keras) ### for DL
library("tensorflow")
library("caretEnsemble")
library(kerasR)
library("h2o")
library(NeuralNetTools)

Example

Preparing data

f <- system.file('extdata',package='DeepGenomeScan')
infile <- file.path(f, "sim1.csv")
sim_example=read.csv(infile)
genotype=sim_example[,-c(1:14)]
env=sim_example[,2:11]
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}
genotype_norm=as.data.frame(apply(genotype,2,normalize))

Setting the resampling method

econtrol1 <- trainControl(## 5-fold CV, repeat 5 times
  method = "adaptive_cv",
  number = 5,
  ## repeated ten times
  repeats = 5,search = "random")
set.seed(999)
options(warn=-1)

DeepGenomeScan with "mlph2o" model

h2o_mlp<- DeepGenomeScan(as.matrix(genotype_norm),env$envir1,
                                  method="mlph2o",
                                  metric = "RMSE",## "Accuracy", "RMSE","Rsquared","MAE"
                                  tuneLength = 10, ### 11 tunable parameters 11^2
                                  # tuneGrid=CNNGrid, ### or search 100 combinations of parameters using random tuneLength=100
                                  trControl = econtrol1)

#### There is a model specific varIMP for this model
varImp(h2o_mlp,scale = FALSE)

out <- as.data.frame(h2o::h2o.varimp(h2o_mlp$finalModel), stringsAsFactors = TRUE)
colnames(out)[colnames(out) == "scaled_importance"] <- "Overall"
rownames(out) <- out$names

Calculating the p-values for SNPs and plot the SNP Manhattan plot

DLqvaluesarsine<-function(DL_data,K)
{
  loadings<-DL_data# [,1:as.numeric(K)]
  normdat <- apply(loadings, 2, normalize)
  asindat=apply(normdat,2, function(x) {asin(sqrt(x))})
  resmaha <- covRob(asindat, distance = TRUE, na.action= na.omit, estim="donostah")$dist
  lambda <- median(resmaha)/qchisq(0.5,df=K)
  reschi2test <- pchisq(resmaha/lambda,K,lower.tail=FALSE)
  qval <- qvalue(reschi2test)
  q.values_DL<-qval$qvalues
  padj <- p.adjust(reschi2test,method="bonferroni")
  return(data.frame(p.values=reschi2test, q.values=q.values_DL,padj=padj,mahaD=resmaha))
}



DLsim1=apply(out,2,normalize) #18

Simqvaluear=DLqvaluesarsine(DLsim1,1)

Manhattan plot

ggplot() +
  geom_point(aes(x=which(Loci=="Neutral"), y=-log10(Simqvaluear[-which(Loci!="Neutral"),1])), col = "gray83") +
  geom_point(aes(x=which(Loci!="Neutral"), y=-log10(Simqvaluear[-which(Loci=="Neutral"),1]), colour = Selected_Loci)) +
  xlab("SNPs") + ylab("DNN -log10(p-value)") +ylim(c(0,100))+theme_bw()
plot(out$Overall,  ylab="SNP importance")

Package tutorial

More about how to construct your own model using DeepGenomeScan can be found in the tutorial

Welcome any feedback and pull request.

Version

The current version is 0.5.5 (Sep 2, 2020).

Contact

qinxinghu@gmail.com

Citation

Qin X, Chiang CWK, Gaggiotti OE. 2022. Deciphering signatures of natural selection via deep learning. Briefings in Bioinformatics. bbac354, https://doi.org/10.1093/bib/bbac354.