前言
今天介绍的这篇文章是关于用深度学习的方法预测G4序列的软件DeepG4，作者将该问题转换成分类问题，文章链接：《DeepG4: A deep learning approach to predict cell-type specific active G-quadruplex regions》
即输入一定的序列特征，利用深度卷积神经网络进行训练，最后输出二分类的结果
关于Keras安装
作者利用R包Keras进行训练，那么这个包常用于深度学习的R包，以tensorflow作为后端，那么该包的安装比较复杂，如果是在windows上安装，需要先安装conda（miniconda也可以，默认安装C盘，安装其他盘请添加环境变量）；如果是Linux系统，同样要先安装conda，并且将用户的conda添加到环境变量里面，接下来在R中安装Keras：
install.packages("remotes")
install.packages("reticulate")
library(reticulate)
remotes::install_github("rstudio/tensorflow",force=T)
library(tensorflow)
# 指定tensorflow的版本和python的版本
install_tensorflow(version = "1.15",conda_python_version = "3.5")
remotes::install_github("rstudio/keras", dependencies = TRUE)
library(keras)
这样就安装好Keras了
在Linux系统中，如果安装好conda，还可以利用conda进行安装：
# 利用conda创建名为 r-tensorflow 的虚拟环境
conda create -n r-tensorflow python=3.6
# 进入r-tensorflow 的虚拟环境
source activate r-tensorflow
# 安装tensorflow和keras
conda install -y tensorflow
conda install -y keras
# 在 r-tensorflow 的虚拟环境中安装R
conda install -y r-base
#安装keras包
conda install -y r-keras
然后进入r-tensorflow 虚拟环境中的R里，library(keras)就可以了
DeepG4源码分析
它的核心代码比较简单，我们以序列特征作为分类指标为例子，作者在这里使用的分类指标是将字符信号转换为数字信号，由于基因组文件都是A,G,C,T四个碱基字符，所以很容易就转换为数字信号：
假设说我有一段AGTCG的序列，转换为数字信号即为一个5×4的array，由于G4的序列特征满足
所以某G4序列形成的array按照0，1排列的顺序，从一定程度上可以反应该序列的特征
以DeepG4的测试数据为例，读入的序列为：
[[1]]
A DNAStringSet instance of length 100
width seq
[1] 201 GTTCGGGCCTCGGTCGCGCCGCCGGGTCTTGCAGACGCGAATGTAAACAGAAACA...TGACTCCTGGAGCGACCTTCACGAGGGAAAGCGCGCCCCCCGGCACCCACCCCT
[2] 201 TTTCTATAGTTTTCTTTTGTTTCTACCTCATGACTAGATGATTCACTGCTTGAAC...GTCAAATCTGTCCATCTTCACTGCCACCCTTCAGTACCAAATGACCAGTCTCTT
[3] 201 GCTTAAAAGCCTGTAAGAAAGATATAATTTGATAGAACTGGCTAGGATTTGTCAG...CGTCAGGGAGGGGGTGGGGCCTCCACGTGGGAGATCTTGCCTGGAGGTGGTGGA
[4] 201 TCCCACACCCGGTAGATGTAAGGGAAAAACTGCATTACCCAGAAGGCACTGCCCC...GTGTGACGTCATCTCCGTGGGCCGGTTTGGCCCTGAAACAGTGTGGGGCCTAGA
[5] 201 AGTAGCTACAGAGTTCCTGCTCCAGCAACCAGGAGCCTTGAGGCAGCACAAGGAC...ACCACAATGTCTGCCAAGAAAGAGGATGAGTCACCAAGACCCACAGGAAAGAGG
... ... ...
[96] 201 CACATGCCTTCCTTGGGGACGTGTTCACACATGTGGCCCTAGCTGTGAGAGACAG...CATCTCAGAACAGCTGAGCTGGAAGTGGGTGAATAATAATAATAATAATAATAA
[97] 201 TGGTGGTCTTTCTCTACCGGGCCTGGTAGCCAAAGACAAAGGTCATAATCACTTG...CTATGTACTCTTCAAAGTGCCACCTCCTGGCTGCAAGCCAACCAACACAAAACC
[98] 201 TGACCGTAGACCTCGTGCACTTCTGCTGCGGTCGGGGCCGGAGTCTGGGCTGGAG...GCGATCCAGAGCCAAGCGCCCCGCCCCTGCCCGGGCGCGCTCCCTCCTTAGCCC
[99] 201 TTAACGTCATCAGTCGGGAGGACGACAGCTACGCACGCGCGGGGCACCTCCTCTG...GCCACGGTGGAGGCAGCGGCGAGAGGGGGCGGGGACAAGGAGAGGGCACGCACG
[100] 201 GTGTCCGGGTGAGAGACCTGGAGGTGGGGCCTAGGTGTCTACCCGGCCAGGTGCG...TAAGGCTCGGGGCCAGTCGTCGTCCATTCCTTCCTAACACCTCCCTATCCTCCC
源码：
library(Biostrings)
library(rtracklayer)
sequences <- readDNAStringSet("test_G4_data.fa")
# 基本的input
X = sequences
X.atac = NULL
Y=NULL
model=NULL
lower.case=F
treshold = 0.5
seq.size = 201
retrain=FALSE
retrain.path=""
log_odds=F
# 定义序列长度
model.regular.size.accepted <- 201
# 对四个碱基编码
tabv = c("T"=4,"G"=3,"C"=2,"A"=1)
# 转换为序列的数字信号 array
X <- DNAToNumerical(X,tabv = tabv,lower.case=lower.case,seq.size = seq.size)
X_inputs <- X
rm(X)
# 载入已经训练好的模型，h5文件
model <- keras::load_model_hdf5('best_model.h5')
model <- keras::keras_model(inputs = model$input,
outputs = keras::get_layer(model, index = log_odds_index)$output)
# 预测
res <- stats::predict(model,X_inputs)
最终的结果为：
结果为一个分数，文章建议该score介于0到1之间为active的G4 regions，而位于该区域外的即为非active的G4 regions。推测该范围是利用正负样本的序列，分别计算它们的score，然后分别看它们的分布，总结出来的经验范围
关于DNAToNumerical函数：
DNAToNumerical <- function(x,tabv = c("T"=4,"G"=3,"C"=2,"A"=1),lower.case=F,seq.size = 201){
x <- Biostrings::as.matrix(x)
listMat <- list()
for(i in 1:length(tabv)){
nuc_index <- tabv[[i]]
nuc_value <- names(tabv[i])
mat <- matrix(0,nrow(x),ncol(x))
mat[x==nuc_value] <- 1
listMat[[nuc_index]] <- mat
}
arrayout <- array(unlist(listMat), dim = c(nrow(listMat[[1]]), ncol(listMat[[1]]), length(listMat)))
return(arrayout)
}
DeepG4支持自己训练数据
DeepG4也支持自己训练模型，并且自己测试：
ATAC <- system.file("extdata", "Peaks_BG4_G4seq_HaCaT_GSE76688_hg19_201b_Accessibility.bw", package = "DeepG4")
ATAC <- import.bw(ATAC)
# Read positive and segative set of sequences
bed.pos <- import.bed(system.file("extdata", "Peaks_BG4_G4seq_HaCaT_GSE76688_hg19_201b.bed", package = "DeepG4"))
bed.neg <- import.bed(system.file("extdata", "Peaks_BG4_G4seq_HaCaT_GSE76688_hg19_201b_Ctrl_gkmSVM.bed", package = "DeepG4"))
# Generate classes
Y <- c(rep(1,length(bed.pos)),rep(0,length(bed.neg)))
BED <- c(bed.pos,bed.neg)
Input_DeepG4 <- DeepG4InputFromBED(BED=BED,ATAC = ATAC,GENOME=BSgenome.Hsapiens.UCSC.hg19)
# 训练模型
training <- DeepG4(X=Input_DeepG4[[1]],X.atac=Input_DeepG4[[2]],Y,retrain=TRUE,retrain.path = "DeepG4_retrained.hdf5")
重新训练样本需要加参数retrain=TRUE和相应保存h5文件的路径retrain.path，我们看其源代码：
# 该model重新训练的model
history <- keras::fit(
model,
X_inputs,
Y,
epochs = 20,
batch_size = 128,
validation_split = 0.2,
verbose= 0)
res <- stats::predict(model,X_inputs)
if(retrain.path == ""){
retrain.path <- paste0("DeepG4_retrained_",Sys.Date(),".hdf5")
}
keras::save_model_hdf5(model,retrain.path)
这样我们就可以使用自己训练的模型进行训练了
DeepG4加入染色质可及性特征进行预测
当然DeepG4也可以将G4区域的内的染色质可及性作为一种指标来训练模型
library(rtracklayer)
library(BSgenome.Hsapiens.UCSC.hg19)
library(DeepG4)
BED <- system.file("extdata", "test_G4_data.bed", package = "DeepG4")
BED <- import.bed(BED)
# 读取染色质可及性的数据
ATAC <- system.file("extdata", "Peaks_BG4_G4seq_HaCaT_GSE76688_hg19_201b_Accessibility.bw", package = "DeepG4")
ATAC <- import.bw(ATAC)
Input_DeepG4 <- DeepG4InputFromBED(BED=BED,ATAC = ATAC,GENOME=BSgenome.Hsapiens.UCSC.hg19)
# 进行预测
predictions <- DeepG4(X=Input_DeepG4[[1]],X.atac = Input_DeepG4[[2]])
