R包ggseqlogo|R包ggseqlogo |绘制序列分析图
简介
在生物信息分析中,经常会做序列分析图(sequence logo),这里的序列指的是核苷酸(DNA/RNA链中)或氨基酸(在蛋白质序列中)。sequence logo图是用来可视化一段序列某个位点的保守性,据根提供的序列组展示位点信息。常用于描述序列特征,如DNA中的蛋白质结合位点或蛋白质中的功能单元。
实现以上可视化过程的工具有很多,本文介绍一个使用起来非常简单,不拖泥带水的R包ggseqlogo,只要你根据此包要求的数据格式上传一堆DNA序列或者氨基酸序列,再根据现成的命令流程就能画出logo图。
安装到作图的代码如下:
安装
安装方式有两种
#直接从CRAN中安装
install.packages("ggseqlogo")
#从GitHub中安装
devtools::install.github("omarwagih/ggseqlogo")
数据加载 ggseqlogo提供了测试数据
ggseqlogo_sample
。#加载包
library(ggplot2)
library(ggseqlogo)
#加载数据
data(ggseqlogo_sample)
ggseqlogo_sample
数据集是一个列表,里面包含了三个数据集:- seqs_dna:12种转录因子的结合位点序列
- pfms_dna:四种转录因子的位置频率矩阵
- seqs_aa:一组激动酶底物磷酸化位点序列
#seqs_dna
head(seqs_dna)[1]
## $MA0001.1
##[1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG"
##[6] "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG"
#pfms_dna
head(pfms_dna)[1]
## $MA0018.2
##[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## A001101028
## C11090370
## G1100210011
## T90000812
#seqs_aa
head(seqs_aa)[1]
## $AKT1
##[1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG"
##[4] "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"
外部数据读入 也可以用自己的数据集,支持两种格式,序列和矩阵。
# 长度为7的motif。每一行为一条序列,长度相同,每一列的碱基组成代表对应位置的碱基偏好性。
fasta = "ACGTATG
ATGTATG
ACGTATG
ACATATG
ACGTACG"fasta_input <- read.table(fasta, header=F, row.names=NULL)
fasta_input <- as.vector(fasta_input$V1)# 长度为5的motif矩阵示例,每一列代表一个位置,及碱基在该位置的出现次数。共4行,每一行代表一种碱基
matrix <- "Base12345
A102081
C112123
G40911
T00019
"
matrix_input <- read.table(matrix, header=T, row.names=1)
matrix_input <- as.matrix(matrix_input)
可视化
ggplot()+geom_logo(seqs_dna$MA0001.1)+theme_logo()
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/d218cc635ae34744b4c5f48630a6f09d.png)
文章图片
image ggseqlogo提供了一个直接绘图的函数
ggseqlogo()
,这是一个包装函数。下面命令结果同上面的。ggseqlogo(seqs_dna$MA0001.1)
输入格式 ggseqlogo支持以下几种类型数据输入:
- 序列
- 矩阵
ggseqlogo(pfms_dna$MA0018.2)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/64f272587c7a4aa2a2b7ca99aa278c63.png)
文章图片
image 方法 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)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/e90cc056677c4b56b7de45c5bb5d0305.jpg)
文章图片
image 序列类型 ggseqlogo支持氨基酸、DNA和RNA序列类型,默认情况下ggseqlogo会自动识别数据提供的序列类型,也可以通过
seq_type
选项直接指定序列类型。ggseqlogo(seqs_aa$AKT1, seq_type="aa")
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/cf28e5e335444a8f87981015fc2d7c7d.jpg)
文章图片
image 自定义字母 通过
namespace
选项来定义自己想要的字母类型#用数字来代替碱基
seqs_numeric <- chartr("ATGC", "1234", seqs_dna$MA0001.1)
ggseqlogo(seqs_numeric, method="prob", namespace=1:4)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/893a0ea04a04471cb918d298e00c6740.jpg)
文章图片
image 配色 ggseqlogo可以使用
col_scheme
参数来设置配色方案,具体可参考?list_col_schemes
ggseqlogo(seqs_dna$MA0001.1, col_scheme="base_pairing")
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/444b0df9c0294414a6bc5a136d3a1b36.jpg)
文章图片
image 自定义配色 ggseqlogo提供函数
make_col_scheme
来自定义离散或者连续配色方案离散配色
csl <- make_col_scheme(chars = c("A","T", "C", "G"), groups = c("gr1","gr1", "gr2","gr2"), cols = c("purple","purple","blue","blue"))
ggseqlogo(seqs_dna$MA0001.1,col_scheme=csl)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/ccda831d742d438490b43e963b8e15dd.png)
文章图片
image 连续配色
cs2 <- make_col_scheme(chars = c("A", "T", "C", "G"), values = 1:4)
ggseqlogo(seqs_dna$MA0001.1, col_scheme=cs2)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/39d50351432f4a08b26f0c83a11e96e5.jpg)
文章图片
image 同时绘制多个序列标志
ggseqlogo(seqs_dna, ncol = 4)
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/7f898fce5a2a44d08c97a2070ac7e6f5.jpg)
文章图片
image 上述命令实际上等同于
ggplot()+geom_logo(seqs_dna)+theme_logo()+
facet_wrap(~seq_group,ncol = 4,scales = "free_x")
自定义高度 通过创建矩阵可以生成每个标志的高度,还可以有负值高度
set.seed(1234)
custom_mat <- matrix(rnorm(20), nrow = 4, dimnames = list(c("A","T","C", "G")))
ggseqlogo(custom_mat,method="custom",seq_type="dna")+
ylab("my custom height")
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/b14ae48f5d1d4e1a8af309c3ff91a7d8.jpg)
文章图片
image 字体 可以通过
font
参数来设置字体,具体可参考?list_fonts
fonts <- list_fonts(F)
p_list <- lapply(fonts, function(f){
ggseqlogo(seqs_dna$MA0001.1,font=f)+ggtitle(f)
})
do.call(gridExtra::grid.arrange,c(p_list, ncol=4))
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/a534700aab504c118479d6d5c2c02eff.jpg)
文章图片
image 注释 注释的话跟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()
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/448ee6c380ff4a17b7df0bb89ac83e8d.jpg)
文章图片
image 图形组合 将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")
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/7b2e79617325402ba71089d9b36fdc58.jpg)
文章图片
image 严涛:浙江大学,作物遗传育种在读博士。爱鼓捣各种可视化软件,爱折腾。点击阅读原文跳转其博客。
【R包ggseqlogo|R包ggseqlogo |绘制序列分析图】更多数据可视化可参考在线工具:http://www.ehbio.com/ImageGP/
![R包ggseqlogo|R包ggseqlogo |绘制序列分析图](https://img.it610.com/image/info10/6a00c050a09a4071a978a89a98d80eda.jpg)
文章图片
image
推荐阅读
- 喂,你结婚我给你随了个红包
- CET4听力微技能一
- 放下心中的偶像包袱吧
- 社保代缴公司服务费包含哪些
- Beego打包部署到Linux
- 世界之大,包罗万象--|世界之大,包罗万象-- 读《我不过低配的人生》
- 用npm发布一个包的教程并编写一个vue的插件发布
- 积极探索|积极探索 绽放生命 ???——心心相印计划:青少年心理工作研讨小组全国大型公益行动第二次活动包头市青山区分校圆满成功
- 那个喝大了的女人在群里发了一晚上的红包
- HttpClient对外部网络的操作