R语言实现DNA结构预测

大家对DNA应该都有一定的了解，那么DNA同样不仅仅是具有一级结构的碱基序列，而且还具有二级结构（双螺旋），三级结构（超螺旋）的特征。今天给大家介绍一个来预测DNA结构的R包DNAshapeR，其从基因组测序数据中以超高速、高通量的方式预测DNA形状特征。该软件包以核苷酸序列或基因组间隔作为输入，并生成各种图形表示，以供进一步分析。DNA预测使用滑动五聚体窗口，其中512个不同五聚体中的每一个都有独特的结构特征，从而在每个核苷酸位置（周向）定义了小沟宽（MGW），滚动，螺旋桨扭曲（ProT）和螺旋扭曲（HelT）的向量（周等人，2013）。MGW和ProT定义碱基对参数，而Roll和HelT代表碱基对步长参数。首先我们看下需要安装的包：

BiocManager::install("DNAshapeR")
BiocManager::install("BSgenome.Scerevisiae.UCSC.sacCer3")#CpGmethylation此包需要外网
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")#referencegenome此包需要外网
BiocManager::install("AnnotationHub")
BiocManager::install("Biostrings")

接下来我们直接通过实例来看下包的操作：

1. 序列的结构预测

library(DNAshapeR)
fn <- system.file("extdata","CGRsample.fa", package = "DNAshapeR")
pred <- getShape(fn)

plotShape(pred$MGW)

plotShape(pred$ProT)

plotShape(pred$Roll)

plotShape(pred$HelT)

2. 对于参考基因组上的间隔进行预测

library(BSgenome.Scerevisiae.UCSC.sacCer3)
 
# Create a query GRanges object
gr <- GRanges(seqnames =c("chrI"),
           strand = c("+", "-", "+"),
           ranges = IRanges(start = c(100, 200, 300), width = 100))
getFasta(gr, Scerevisiae, width = 100,filename = "tmp.fa")
fn <- "tmp.fa"
pred <- getShape(fn)
heatShape(pred$ProT, 20)#其中20指的是色块的合并需要的大小即热图的单元格数量

3. 利用公共的平台去检索对应的基因序列做结构预测,当然，我们参考的这个平台的包需要连接外网。

library(BSgenome.Hsapiens.UCSC.hg19)
library(AnnotationHub)
 
ah <- AnnotationHub()

ah <- subset(ah, species=="Homosapiens")
ah <- query(ah, c("H3K4me3","Gm12878", "Roadmap"))
getFasta(ah[[1]], Hsapiens, width = 150,filename = "tmp.fa")

fn <- "tmp.fa"
pred <- getShape(fn)

4. 预测CpG甲基化背景下的DNA形态特征，可以制备一个FASTA序列文件，其中符号Mg: M表示前导链上甲基化CpG的胞嘧啶，符号g表示后随链上甲基化CpG的胞嘧啶。例如，

>seq1
GTGTCACMgCGTCTATACG

library(DNAshapeR)
fn_methy <- system.file("extdata","MethylSample.fa", package = "DNAshapeR")
pred_methy <- getShape(fn_methy,methylate = TRUE)

5. 基因组浏览器的类似功能，展示DNA形状的预测结构

fn2 <- system.file("extdata","SingleSeqsample.fa", package = "DNAshapeR")
pred2 <- getShape(fn2)
trackShape(fn2, pred2) # Only for singlesequence file

#如上图，横轴是序列中的每一个碱基的序列，纵轴是对每个点的预测值。

6. DNA序列的编码

library(Biostrings)
fn3 <- system.file("extdata","PBMsample_short.fa", package = "DNAshapeR")
pred3 <- getShape(fn3)
featureType <- c("1-mer","1-shape")
featureVector <- encodeSeqShape(fn3,pred3, featureType)
head(featureVector)

得到DNA序列的编码，后面具体如何应用那就很广了，在此不再赘述。