我们前面讲过在python中如何实现测序图标（sequence logos）的绘制。今天给大家介绍一个在R语言中实现DNA，RNA以及氨基酸的logos绘制的R包motifStack。首先我们看下包的安装，主要是通过bioconductor进行安装，具体的代码不再赘述，请参看bioconductor官网。如果能完整展示绘图还需要另外一个ghostscript的软件，其官网（https://www.ghostscript.com/）：
下载和自己系统对应的安装包，然后进行那些下一步………。软件就安好好了。接下来需要我们设置下环境变量，那我们就一个图展示下：
至此前期工作准备完了，然后就是重启电脑。如果还是无法画图那就可以在运行绘图时，前面直接运行如下代码：
代码语言：javascript复制Sys.setenv(R_GSCMD=“F:/gs9.27/bin/ gswin64.exe”)接下来我们直接看此包是如何实现logos的绘制的：
1. DNA sequence logos绘制：
代码语言：javascript复制##数据的读入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的绘制：
代码语言：javascript复制rna <- pcmrownames(rna)[4] <- "U"motif <- new("pcm",mat=as.matrix(rna), name="RNA_motif")plot(motif)
3. protein logos绘制：
代码语言：javascript复制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的绘制：
代码语言：javascript复制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：
代码语言：javascript复制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")
那如果想美化，那么可以参考下面的代码：
代码语言：javascript复制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：
代码语言：javascript复制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. 多组放射树状图绘制：
代码语言：javascript复制## 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图的绘制：
代码语言：javascript复制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. 树状图热图的组合绘制：
代码语言：javascript复制## 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)
当然此包还提供了JS的网页前端展示，在此我们就不再去描述。
