可视化motif logo三种方法(seqLogo/ggseqlogo/LogoJS)

# 安装
source ( "http://bioconductor.org/biocLite.R" )
biocLite("seqLogo")

# 导入CTCF的PPM矩阵
library(seqLogo)
ctcf <- read.table('~/Documents/writethings/20200409/ppm.txt', header = F)

# 转为PWM矩阵（转置是因为必须要让行分别是ATCG的PWM score）
trans_ctcf <- makePWM(t(ctcf))
seqLogo(trans_ctcf)

bad = (sapply( body(seqLogo), "==", "grid.newpage()") |
        sapply( body(seqLogo), "==", "par(ask = FALSE)"))
body(seqLogo)[bad] = NULL
tf_names <- c('CTCF', 'SP1')
grid.newpage()
for(i in 0:1){

   ppm <- read.table(paste('~/Documents/writethings/20200409/ppm',i+1 ,'.txt', sep = ''), header = F)
   pwm <- makePWM(t(ppm))

   pushViewport(viewport(x=0.25+0.55*(i%%2),
                         y=0.5+0.55*(i%/%2),
                         width=0.5, height=0.6))
   seqLogo(pwm)
   grid.text(tf_names[i+1], x=0.5, y=1,
             hjust=0.5, vjust=1)
   popViewport()

}

library(ggseqlogo)
ppm <- t(ppm) 
rownames(ppm) <- c("A", "G", "C", "T")
ggseqlogo(ppm)