文章/答案/技术大牛

发布

社区首页 >专栏 >R语言实现sequence logos绘制

R语言实现sequence logos绘制

一粒沙

发布于 2019-09-08 06:24:00

1.6K00

代码可运行

文章被收录于专栏：R语言交流中心R语言交流中心

运行总次数：0

代码可运行

我们前面讲过在python中如何实现测序图标（sequence logos）的绘制。今天给大家介绍一个在R语言中实现DNA，RNA以及氨基酸的logos绘制的R包motifStack。首先我们看下包的安装，主要是通过bioconductor进行安装，具体的代码不再赘述，请参看bioconductor官网。如果能完整展示绘图还需要另外一个ghostscript的软件，其官网（https://www.ghostscript.com/）：

下载和自己系统对应的安装包，然后进行那些下一步………。软件就安好好了。接下来需要我们设置下环境变量，那我们就一个图展示下：

至此前期工作准备完了，然后就是重启电脑。如果还是无法画图那就可以在运行绘图时，前面直接运行如下代码：

Sys.setenv(R_GSCMD=“F:/gs9.27/bin/ gswin64.exe”)

接下来我们直接看此包是如何实现logos的绘制的：

1. DNA sequence logos绘制：

##数据的读入library("motifStack")pcm <-read.table(file.path(find.package("motifStack"),                           "extdata", "bin_SOLEXA.pcm"))pcm <- pcm[,3:ncol(pcm)]rownames(pcm) <-c("A","C","G","T")

##数据模型构建motif <- new("pcm",mat=as.matrix(pcm), name="bin_SOLEXA")#name就是绘图的标题

##将所有的属性进行绘制对比opar<-par(mfrow=c(4,1))plot(motif)#plot the logo with same heightplot(motif, ic.scale=FALSE,ylab="probability")#try a different fontplot(motif, font="mono,Courier")#try a different font and a different colorgroupmotif@color <-colorset(colorScheme='basepairing')plot(motif,font="Times")

2. RNA sequence logos的绘制：

rna <- pcmrownames(rna)[4] <- "U"motif <- new("pcm",mat=as.matrix(rna), name="RNA_motif")plot(motif)

3. protein logos绘制：

protein<-read.table(file.path(find.package("motifStack"),"extdata","cap.txt"))protein<-t(protein[,1:20])

motif<-pcm2pfm(protein)motif<-new("pfm", mat=motif,name="CAP",           color=colorset(alphabet="AA",colorScheme="chemistry"))plot(motif)

3. position-specific affinity matrix logos的绘制：

motif<-matrix(  c(   .846, .631, .593, .000, .000, .000, .434, .410, 1.00, .655, .284, .000,.000, .771, .640, .961,   .625, .679, .773, 1.00, 1.00, .000, .573, .238, .397, 1.00, 1.00, .000,.298, 1.00, 1.00, .996,   1.00, 1.00, 1.00, .228, .000, 1.00, 1.00, .597, .622, .630, .000, 1.00,1.00, .871, .617, 1.00,   .701, .513, .658, .000, .000, .247, .542, 1.00, .718, .686, .000, .000,.000, .595, .437, .970  ),nrow=4, byrow = TRUE)rownames(motif) <- c("A","C", "G", "T")motif<-new("psam", mat=motif,name="affinity logo")
plot(motif)

以上都是对单个样本的绘制，那么多个样本的绘制就会用到下面的函数motifStack ，其可以绘制"stack", "tree", "phylog","radialPhylog"四种形式的logos：

motifs<-importMatrix(dir(file.path(find.package("motifStack"),"extdata"),"pcm$", full.names = TRUE))

## plot stacksmotifStack(motifs,layout="stack", ncex=1.0)

motifStack(motifs, layout="tree")

motifStack(motifs,layout="radialPhylog")

那如果想美化，那么可以参考下面的代码：

library(RColorBrewer)color <- brewer.pal(12,"Set3")## plot logo stack with radial stylemotifStack(motifs,layout="radialPhylog",          circle=0.3, cleaves = 0.2,          clabel.leaves = 0.5,          col.bg=rep(color, each=5), col.bg.alpha=0.3,          col.leaves=rep(color, each=5),          col.inner.label.circle=rep(color, each=5),          inner.label.circle.width=0.05,          col.outer.label.circle=rep(color, each=5),          outer.label.circle.width=0.02,          circle.motif=1.2,          angle=350)

除此以外，此包还提供了更复杂的绘图：

1. Sequence logos 云图需要用到函数motifCloud：

library("MotifDb")matrix.fly <- query(MotifDb,"Dmelanogaster")motifs2 <- as.list(matrix.fly)## use data from FlyFactorSurveymotifs2 <-motifs2[grepl("Dmelanogaster\\-FlyFactorSurvey\\-",                         names(motifs2))]## format the namesnames(motifs2) <-gsub("Dmelanogaster_FlyFactorSurvey_", "",                      gsub("_FBgn\\d+$", "",                           gsub("[^a-zA-Z0-9]","_",                                 gsub("(_\\d+)+$","", names(motifs2)))))motifs2 <-motifs2[unique(names(motifs2))]pfms <- sample(motifs2, 50)## creat a list of object of pfmmotifs2 <- lapply(names(pfms),                  function(.ele,pfms){new("pfm",mat=pfms[[.ele]], name=.ele)}                  ,pfms)## trim the motifsmotifs2 <- lapply(motifs2, trimMotif,t=0.4)groups <-rep(paste("group",1:5,sep=""), each=10)names(groups) <- names(pfms)## assign group colorsgroup.col <- brewer.pal(5,"Set3")names(group.col)<-paste("group",1:5,sep="")## use MotIV to calculate the distances ofmotifsjaspar.scores <-MotIV::readDBScores(file.path(find.package("MotIV"),                                               "extdata",                                              "jaspar2010_PCC_SWU.scores"))d <- MotIV::motifDistances(lapply(pfms,pfm2pwm))hc <- MotIV::motifHclust(d,method="average")## convert the hclust to phylog objectphylog <- hclust2phylog(hc)## reorder the pfms by the order of hclustleaves <- names(phylog$leaves)pfms <- pfms[leaves]## create a list of pfm objectspfms <- lapply(names(pfms),function(.ele, pfms){                               new("pfm",mat=pfms[[.ele]], name=.ele)}               ,pfms)## extract the motif signaturesmotifSig <- motifSignature(pfms, phylog,groupDistance=0.01, min.freq=1)## draw the motifs with a tag-cloud style.motifCloud(motifSig, scale=c(6, .5),          layout="rectangles",          group.col=group.col,          groups=groups,          draw.legend=TRUE)

2. 多组放射树状图绘制：

## get the signatures from object ofmotifSignaturesig <- signatures(motifSig)## set the inner-circle color for eachsignaturegpCol <- sigColor(motifSig)## plot the logo stack with radial style.plotMotifStackWithRadialPhylog(phylog=phylog,pfms=sig,                              circle=0.4,cleaves = 0.3,                              clabel.leaves =0.5,                              col.bg=rep(color,each=5), col.bg.alpha=0.3,                             col.leaves=rep(rev(color), each=5),                             col.inner.label.circle=gpCol,                              inner.label.circle.width=0.03,                              angle=350,circle.motif=1.2,                             motifScale="logarithmic")

3. Circles图的绘制：

motifCircos(phylog=phylog, pfms=pfms,pfms2=sig,           col.tree.bg=rep(color, each=5), col.tree.bg.alpha=0.3,           col.leaves=rep(rev(color), each=5),           col.inner.label.circle=gpCol,           inner.label.circle.width=0.03,           col.outer.label.circle=gpCol,           outer.label.circle.width=0.03,           r.rings=c(0.02, 0.03, 0.04),           col.rings=list(sample(colors(), 50),                           sample(colors(),50),                           sample(colors(),50)),           angle=350, motifScale="logarithmic")

4. 树状图热图的组合绘制：

## plot the logo stack with pile style.motifPiles(phylog=phylog, pfms=pfms,pfms2=sig,           col.tree=rep(color, each=5),           col.leaves=rep(rev(color),each=5),           col.pfms2=gpCol,           r.anno=c(0.02, 0.03, 0.04),           col.anno=list(sample(colors(), 50),                          sample(colors(), 50),                          sample(colors(),50)),           motifScale="logarithmic",           plotIndex=TRUE,           groupDistance=0.01)