ggseqlogo

#文章完全来源于该篇博文：
https://www.plob.org/article/12630.html

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

• 安装

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

• 演示数据集

#加载包
library(ggplot2)
library(ggseqlogo)
#加载数据
data(ggseqlogo_sample)

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

• seqs_dna:12种转录因子的结合位点序列
• pfms_dna:四种转录因子的位置频率矩阵
• seqs_aa:一组激动酶底物磷酸化位点序列

#seqs_dna
head(seqs_dna)[1]
 
## $MA0001.1<br>## [1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG"<br>## [6] "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG"
 
#pfms_dna<br>head(pfms_dna)[1]
 
## $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
 
#seqs_aa
head(seqs_aa)[1]
 
## $AKT1
## [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG"
## [4] "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"

• 可视化

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

• 输入格式
  ggseqlogo支持以下几种类型数据输入：
  • 序列
  • 矩阵

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

ggseqlogo(pfms_dna$MA0018.2)

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

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

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

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

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

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

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

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

• 支持分组和刻面

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

• 注释
  注释的话跟ggplot2是一样的

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

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

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