利用motifStack将homer结果可视化#1

#本人也不熟悉R语言，完全照猫画虎，有错误的地方还请指出，非常感谢！

前面利用homer将我们的基因集进行了启动子的motif分析，最后输出的结果中有以下几类：

由于我没有进行富集分析，所以最终输出的结果中也没有富集的结果，其中de nove的结果和knowmotif的结果分别在两个文件夹中，接下来的示例是用knowmotif的结果来展示。

点开文件夹后出来的是单个的结果展示，分别是homer输出的logo图和ATCG分布的矩阵文件.motif文件，后续用到的就是这个motif文件。

先打开motif文件可以看到输出的具体结果，保守的具体序列和预测出最佳匹配的对象，这里我只保留了对象名SPL5

但是，从其他人的教程中我看到的基本上是pcm文件，所以该如何对motif文件进行转换呢？这里我照搬了一个博主的方案m6A可视化，用MotifStack美化或自定义HOMER的结果 – 王进的个人网站 (jingege.wang)

他在博客中写了如何将motif的结果转换为后续分析的标准格式

按照描述，我先把文件中第一行删除，只留下矩阵数字等

现在开始进行转化格式得操作

首先安装motifstack包，运行前先调用BiocManager这个包，没有得话就先在线安装。

BiocManager::install("motifStack")

安装完成后调用motifStack包。但是出现报错，缺少一个依赖包，因此还需要再安装‘GenomicAlignments’这个包。仍然使用BiocManager来安装。

BiocManager::install("GenomicAlignments")

但是又出现报错

检查了问题，是因为我无法访问这个网站，解决得办法有两种，一种是去把包下载到本地，进行本地安装，另一种是开个vpn，再进行安装。这里我使用的是vpn，再尝试安装就成功了。

后续出现的依赖包缺失都用此方法进行安装。由于我第一次使用R，所以里面有很多包我都没安装，以下是这次运行时所安装的包的合集。

BiocManager::install("GenomicAlignments")
BiocManager::install("BiocIO")
BiocManager::install("annotate")
BiocManager::install("KEGGREST")
BiocManager::install("GO.db")

由于需要把外部motif文件读入到程序中，因此我们首先要知道文件存在的路径，可以手敲，还有一种方便的方法来获取文件路径，file.choose()，运行这条命令后会弹出文件选择框，选择好文件后会列出文件的路径信息。

file.choose()

确定路径信息后，开始读入文件数据。但是需要注意的是，这里仍然需要手动修改一下，因为R识别的路径信息和windows的路径命名符号有差别，需要将\\改为/，并且需要在英文输入法下进行修改。

motiffile <- read.table('/Users/JHK/Desktop/motif-test.txt')

录入成功后在Environment中查看一下数据

homer导出的文件数据结构和jaspar中的结构不一样，具体是homer的四列分别代表ACGT四个碱基在不同位置上的贡献值，而jaspar则是分布在四行，并非四列。这里目前四列的还没有命名，系统给出的默认命名时V1，V2，V3，V4所以，我们先给四列进行命名，分别是ACGT。

names(motiffile) <- c("A","C","G","T")

之后是给四列元素赋值，将对应列的数据分别赋值给ACGT

A <- motiffile[,1]
C <- motiffile[,2]
G <- motiffile[,3]
T <- motiffile[,4]

现在就有了这个motif中四个碱基的具体赋值结果，接下来需要将四个结果排列成jaspar或者说pcm的结构。

data <- rbind(A,C,G,T)

查看排列好的结果。（这里我感觉用excel的转置复制就行了。。。。。）

将结果赋值给pcm

pcm <- data[,1:ncol(data)]

再给每一行的数据命名

rownames(pcm) <- c("A","C","G","T")

将上述表格转化为矩阵文件，命名为motif-test

motif <- new("pcm", mat=as.matrix(pcm), name="motif-test")

最后就是绘制logo图

opar<-par(mfrow=c(4,1))
par(opar)
plot(motif)

完整的代码如下，具体请参考上面给的连接，我只是照猫画虎做了一次

BiocManager::install("motifStack")
library(motifStack)
library(grid)
library(motifStack)
BiocManager::install("GenomicAlignments")
BiocManager::install("BiocIO")
BiocManager::install("annotate")
BiocManager::install("KEGGREST")
BiocManager::install("GO.db")

file.choose()
motiffile <- read.table('/Users/JHK/Desktop/motif-test.txt')
names(motiffile) <- c("A","C","G","T")
A <- motiffile[,1]
C <- motiffile[,2]
G <- motiffile[,3]
T <- motiffile[,4]
data <- rbind(A,C,G,T)
pcm <- data[,1:ncol(data)]
rownames(pcm) <- c("A","C","G","T")
motif <- new("pcm", mat=as.matrix(pcm), name="motif-test")
opar<-par(mfrow=c(4,1))
par(opar)
plot(motif)

输出的结果，可以保存为其他格式

写在后面：

到这里我们就成功画出了一个motif的logo图，但是我具体的需求是对多个motif结果进行比较，除了绘制出logo图以外，还需要对motif进行聚类分析，要有类似于进化树一样的结果。