Motif可视化——从PFM矩阵到sequence logo
Sequence Logo是如何画出来的呢？
计算PFM
还是以MEME example中第一个的Motif结果举个栗子~
根据包含该Motif的各序列中，统计每个位点四种核苷酸出现的次数，并计算频数即得到Position Frequency Matrix（PFM）。
以第一个位点为例，A出现了2次，T出现了12次，因此矩阵的第一列A计数为2，T计数为12，C和G计数为0，其他位点以此类推。
计算每个位点的四个核苷酸频率，得到PPM矩阵，该文件一般是一个N行4列的矩阵，矩阵中描述了每一个碱基在每一个位置上出现的频率。以第一行为例，A出现的频率为2/(2+12)=0.142857；T出现的频率为14/(2+12)=0.857143。假设PFM各个位点之间相互独立，那么例如该Motif序列为：TACTGTATATAHAHMCAG，则该序列的出现的概率为碱基在每个位点上的概率乘积，0.86*0.86*1*1*1*0.93……以此类推。
PPM矩阵文件一般在网页版的结果中可直接下载，点击Submit/Download，选择Probability Matrix的数据格式。
获取PWM
接下来，我们要获取到位置权重矩阵position weight matrix（PWM）。用下方公式中k表示A/C/G/T，j表示位点，b表示背景碱基频率，M表示PPM矩阵中的碱基频率。
背景碱基含量在网页结果的最后一部分中有统计结果，如下图所示。
因此可根据公式将PPM矩阵转化为PWM矩阵，如下图所示。以第一行第一列为例：log2(0.142857/0.29)=-1.02148117。根据PWM矩阵可以计算序列的得分，即将碱基在各个位点的在PWM矩阵中的值进行加和，从而可以判断该序列是一段随机序列还是功能位点。例如上述序列TACTGTATATAHAHMCAG的的得分为1.56+1.56+2.25+1.79+2.25+1.68+……，该得分如果大于0，则认为该序列是一个潜在的功能位点，预测到该Motif；如果得分小于0，则认为该序列是一段随机序列；如果等于0，则认为各有50%的概率。
绘制seqlog
R包ggseqlogo，可基于PPM矩阵进行Motif的可视化。
library(ggseqlogo)
library(ggplot2)
Motif <- t(read.table("ppm.txt"))
rownames(Motif) <- c("A","C","G","T")
##list_col_schemes(v = T)查看配色
p1 <- ggseqlogo(Motif,method="prob",col_scheme="base_pairing")
p2 <- ggseqlogo(Motif,method="bits",col_scheme="nucleotide")
gridExtra::grid.arrange(p1,p2)
参考资料
Position weight matrix - Wikipedia ( https://en.wikipedia.org/wiki/Position_weight_matrix )
作者：snail
审稿：童蒙
编辑：amethyst
往期回顾