利用ggseqlogo绘制seqlogo图

来源：TaoYan评论5,897

利用ggseqlogo绘制seqlogo图-图片1

简介

sequence logo图用来可视化一段序列某个位点的保守性，据根提供的序列组展示位点信息。这方面有很多在线小工具可以完成，这里使用R包ggseqlogo进行可视化。

安装

安装方式有两种

  1. #直接从CRAN中安装
  2. install.packages("ggseqlogo")
  3. #从GitHub中安装
  4. devtools::install.github("omarwagih/ggseqlogo")

数据加载

ggseqlogo提供了测试数据ggseqlogo_sample。

  1. #加载包
  2. library(ggplot2)
  3. library(ggseqlogo)
  4. #加载数据
  5. data(ggseqlogo_sample)

ggseqlogo_sample数据集是一个列表，里面包含了三个数据集：

  • seqs_dna:12种转录因子的结合位点序列
  • pfms_dna:四种转录因子的位置频率矩阵
  • seqs_aa:一组激动酶底物磷酸化位点序列
  1. #seqs_dna
  2. head(seqs_dna)[1]
  3.  
  4. ## $MA0001.1<br>## [1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG"<br>## [6] "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG"
  5.  
  6. #pfms_dna<br>head(pfms_dna)[1]
  7.  
  8. ## $MA0018.2<br>## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]<br>## A 0 0 11 0 1 0 2 8<br>## C 1 1 0 9 0 3 7 0<br>## G 1 10 0 2 10 0 1 1<br>## T 9 0 0 0 0 8 1 2
  9.  
  10. #seqs_aa
  11. head(seqs_aa)[1]
  12.  
  13. ## $AKT1
  14. ## [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG"
  15. ## [4] "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"

可视化

  1. ggplot()+geom_logo(seqs_dna$MA0001.1)+theme_logo()

利用ggseqlogo绘制seqlogo图-图片2

ggseqlogo提供了一个直接绘图的函数ggseqlogo(),这是一个包装函数。下面命令结果同上面的。

  1. ggseqlogo(seqs_dna$MA0001.1)

输入格式

ggseqlogo支持以下几种类型数据输入：

  • 序列
  • 矩阵

下面是使用数据中的位置频率矩阵生成的seqlogo

  1. ggseqlogo(pfms_dna$MA0018.2)

利用ggseqlogo绘制seqlogo图-图片3

方法

ggseqlogo通过method选项支持两种序列标志生成方法：bits和probability。

  1. p1 <- ggseqlogo(seqs_dna$MA0001.1, method="bits")
  2. p2 <- ggseqlogo(seqs_dna$MA0001.1, method="prob")
  3. gridExtra::grid.arrange(p1,p2)

利用ggseqlogo绘制seqlogo图-图片4

序列类型

ggseqlogo支持氨基酸、DNA和RNA序列类型，默认情况下ggseqlogo会自动识别数据提供的序列类型，也可以通过seq_type选项直接指定序列类型。

  1. ggseqlogo(seqs_aa$AKT1, seq_type="aa")

 

利用ggseqlogo绘制seqlogo图-图片5

自定义字母

通过namespace选项来定义自己想要的字母类型

  1. #用数字来代替碱基
  2. seqs_numeric <- chartr("ATGC", "1234", seqs_dna$MA0001.1)
  3. ggseqlogo(seqs_numeric, method="prob", namespace=1:4)

 

利用ggseqlogo绘制seqlogo图-图片6

配色

ggseqlogo可以使用col_scheme参数来设置配色方案，具体可参考?list_col_schemes

  1. ggseqlogo(seqs_dna$MA0001.1, col_scheme="base_pairing")

利用ggseqlogo绘制seqlogo图-图片7

自定义配色

ggseqlogo提供函数make_col_scheme来自定义离散或者连续配色方案

离散配色

  1. csl <- make_col_scheme(chars = c("A","T", "C", "G"), groups = c("gr1","gr1", "gr2","gr2"), cols = c("purple","purple","blue","blue"))
  2. ggseqlogo(seqs_dna$MA0001.1,col_scheme=csl)

 

利用ggseqlogo绘制seqlogo图-图片8

连续配色

  1. cs2 <- make_col_scheme(chars = c("A", "T", "C", "G"), values = 1:4)
  2. ggseqlogo(seqs_dna$MA0001.1, col_scheme=cs2)

 

利用ggseqlogo绘制seqlogo图-图片9

同时绘制多个序列标志

  1. ggseqlogo(seqs_dna, ncol = 4)

利用ggseqlogo绘制seqlogo图-图片10
上述命令实际上等同于

  1. ggplot()+geom_logo(seqs_dna)+theme_logo()+facet_wrap(~seq_group,ncol = 4,scales = "free_x")

自定义高度

通过创建矩阵可以生成每个标志的高度，还可以有负值高度

  1. set.seed(1234)
  2. custom_mat <- matrix(rnorm(20), nrow = 4, dimnames = list(c("A","T","C", "G")))
  3. ggseqlogo(custom_mat,method="custom",seq_type="dna")+
  4. ylab("my custom height")

利用ggseqlogo绘制seqlogo图-图片11

字体

可以通过font参数来设置字体，具体可参考?list_fonts

  1. fonts <- list_fonts(F)
  2. p_list <- lapply(fonts, function(f){
  3. ggseqlogo(seqs_dna$MA0001.1,font=f)+ggtitle(f)
  4. })
  5. do.call(gridExtra::grid.arrange,c(p_list, ncol=4))

利用ggseqlogo绘制seqlogo图-图片12

注释

注释的话跟ggplot2是一样的

  1. ggplot()+
  2. annotate("rect", xmin = 0.5, xmax = 3.5, ymin = -0.05, ymax = 1.9, alpha=0.1, col="black", fill="yellow")+
  3. geom_logo(seqs_dna$MA0001.1, stack_width = 0.9)+
  4. annotate("segment", x=4, xend = 8, y=1.2, yend = 1.2, size=2)+
  5. annotate("text", x=6, y=1.3, label="Text annotation")+
  6. theme_logo()

利用ggseqlogo绘制seqlogo图-图片13

图形组合

将ggseqlogo生成的图形与ggplot2生成的图形组合在一起。

  1. p1 <- ggseqlogo(seqs_dna$MA0008.1)+theme(axis.text.x = element_blank())
  2. aln <- data.frame(
  3. letter=strsplit("AGATAAGATGATAAAAAGATAAGA", "")[[1]],
  4. species=rep(c("a","b","c"), each=8),
  5. x=rep(1:8,3)
  6. )
  7. aln$mut <- "no"
  8. aln$mut[c(2,15,20,23)]="yes"
  9. p2 <- ggplot(aln, aes(x, species)) +
  10. geom_text(aes(label=letter, color=mut, size=mut)) +
  11. scale_x_continuous(breaks=1:10, expand = c(0.105, 0)) + xlab('') +
  12. scale_color_manual(values=c('black', 'red')) +
  13. scale_size_manual(values=c(5, 6)) +
  14. theme_logo() +
  15. theme(legend.position = 'none', axis.text.x = element_blank())
  16. bp_data <- data.frame(
  17. x=1:8,
  18. conservation=sample(1:100, 8)
  19. )
  20. p3 <- ggplot(bp_data, aes(x, conservation))+
  21. geom_bar(stat = "identity", fill="grey")+
  22. theme_logo()+
  23. scale_x_continuous(breaks = 1:10, expand = c(0.105, 0))+
  24. xlab("")
  25. suppressMessages(require(cowplot))
  26. plot_grid(p1,p2,p3,ncol = 1, align = "v")

利用ggseqlogo绘制seqlogo图-图片14

SessionInfo

  1. sessionInfo()
  2. ## R version 3.4.3 (2017-11-30)
  3. ## Platform: x86_64-pc-linux-gnu (64-bit)
  4. ## Running under: Ubuntu 17.10
  5. ##
  6. ## Matrix products: default
  7. ## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
  8. ## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
  9. ##
  10. ## locale:
  11. ## [1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C
  12. ## [3] LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8
  13. ## [5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8   
  14. ## [7] LC_PAPER=zh_CN.UTF-8 LC_NAME=C
  15. ## [9] LC_ADDRESS=C LC_TELEPHONE=C
  16. ## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C       
  17. ##
  18. ## attached base packages:
  19. ## [1] stats graphics grDevices utils datasets methods base
  20. ##
  21. ## other attached packages:
  22. ## [1] cowplot_0.9.2 ggseqlogo_0.1 ggplot2_2.2.1
  23. ##
  24. ## loaded via a namespace (and not attached):
  25. ## [1] Rcpp_0.12.15 knitr_1.20 magrittr_1.5 munsell_0.4.3
  26. ## [5] colorspace_1.3-2 rlang_0.2.0 stringr_1.3.0 plyr_1.8.4
  27. ## [9] tools_3.4.3 grid_3.4.3 gtable_0.2.0 htmltools_0.3.6
  28. ## [13] yaml_2.1.16 lazyeval_0.2.1 rprojroot_1.3-2 digest_0.6.15
  29. ## [17] tibble_1.4.2 gridExtra_2.3 evaluate_0.10.1 rmarkdown_1.8
  30. ## [21] labeling_0.3 stringi_1.1.6 compiler_3.4.3 pillar_1.1.0
  31. ## [25] scales_0.5.0 backports_1.1.2

 

发表评论

匿名网友

拖动滑块以完成验证
加载失败