R语言绘图包9--序列分析图的绘制ggseqlogo

2022-10-17 本文已影响0人 Hayley笔记

R语言绘图包系列：

序列分析图（sequence logo）的序列指的是核苷酸(DNA/RNA链中)或氨基酸(在蛋白质序列中)。sequence logo图是用来可视化一段序列某个位点的保守性，根据提供的序列组展示位点信息。常用于描述序列特征，如DNA中的蛋白质结合位点或蛋白质中的功能单元。

ggseqlogo包提供了一种简便的序列分析图绘制方法。其主函数为ggseqlogo()，该函数可以直接接受已经比对好的DNA或者氨基酸序列，但是有一点需要注意：比对好的序列中可以有gap，但是所有序列的长度必须相同，否则会报错。

1. 安装

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

2. 加载数据集

演示数据集

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

ggseqlogo_sample数据集包含了三个列表：
pfms_dna：四种转录因子的位置频率矩阵
seqs_aa：一组激动酶底物磷酸化位点序列
seqs_dna：12种转录因子的结合位点序列

#pfms_dna
head(pfms_dna)[1]
## $MA0018.2
##   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## A    0    0   11    0    1    0    2    8
## C    1    1    0    9    0    3    7    0
## G    1   10    0    2   10    0    1    1
## T    9    0    0    0    0    8    1    2
#seqs_aa
head(seqs_aa[[1]])
# [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG" "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"
#seqs_dna
seqs_dna$MA0001.1
#  [1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG" "CCAAATATAG" "CCAAATATAG" "CCATAAATGG"
# [14] "CCATAAATGG" "CCATAAATGG" "CCATAAATGG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAATATGG" "CCAAAAATGG" "CCAAAAATGG" "CCATTTATAG" "CCATTTATAG"
# [27] "CCATTTATAG" "CCATTAATAG" "CCATTAATAG" "CCATTAATAG" "CCAATTATAG" "CCATTAATGG" "CCAATAATAG" "CCAATTATGG" "CCAAAATTAG" "CCAAAATTAG" "CCAAAATTAG" "CCACATATAG" "CCATTTTTAG"
# [40] "CCATTTTTAG" "CCATTTTTAG" "CCATAAGTGG" "CTATATATAG" "CCATACATGG" "CTATAAATAG" "CTATAAATAG" "CCTAATATAG" "CCTAATATAG" "CCTAATATAG" "CCAAAACTAG" "CTAAATATAG" "CTAAATATAG"
# [53] "CCATAAACAG" "CTAAAAATAG" "CTAAAAATAG" "CCTAAAATGG" "CCTAAAATGG" "CCAGATTTAG" "CCTTTTATAG" "CCTTTTATAG" "CCAAAAAAAG" "CCAAATAAGG" "CCAGAATTAG" "CCATTTCTGG" "CTAAATTTAG"
# [66] "CCCGATATAG" "CCAAAAATAC" "CCATATTGGG" "CCGTTTTTAG" "CAATTTATAG" "CTAAAAGTAG" "CTAAAAGTAG" "CCATTTTAAG" "CCATTTTAAG" "CCATTTATGT" "CTGTAAATAG" "CCATTTTGGG" "CAATAATTAG"
# [79] "CTCAAAATAG" "CCAATTTAGG" "CCAATTTAGG" "CCTAATACGG" "CTATATTAAG" "CTATAAATGT" "CCTTTTATGT" "CCCTTTTTCG" "CTAAAAGAAG" "CTATATTTGT" "CTATAATTGT" "CTACTTAAAG" "CCGAAAAATG"
# [92] "CCTAAAGTGT" "CCTTTTCTAC" "CCCTCACAAG" "TAATTTCTCG" "GCACACAGAC" "TTAGGATTCG"

可以看到seqs_dna中的MA0001.1对象其实是一个含有97个碱基序列的向量，而且，这97个碱基序列其实都是独立的字符串元素。因此，如果大家利用R随机生成一串字符数相同的字符串向量，也是可以用来绘制seqlogo的。

也可以导入自己的数据
支持两种格式，序列和矩阵。

# 长度为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 <- "Base    1    2    3    4    5
A    10    2    0    8    1
C    1    12    1    2    3
G    4    0    9    1    1
T    0    0    0    1    9
"
matrix_input <- read.table(matrix, header=T, row.names=1)
matrix_input <- as.matrix(matrix_input)

3. 可视化

使用数据中的位置频率矩阵生成seqlogo
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），结果和p1是一样的
ggplot()+geom_logo(seqs_dna$MA0001.1)+theme_logo()

在上面的seqlogo中，横轴是每一个碱基的位点，纵轴默认情形下是Bits (也就是在不设置method时，默认出的是上面那张图)。关于bits和probability的区别，在此不做赘述，大家可以去维基百科查询。总之，记住一点，不管是什么指标，都是用来衡量某个位点某个碱基出现的可能性的大小。显然，碱基的字符越高（宽度是相同的），其出现的可能性越大。

使用矩阵绘图

df <- data.frame(seq=sample(c('A','T','C','G'),1000,replace = T),type=rep(c('gt1','gt2','gt3','gt4'),each=250))
mat <- table(df$seq,df$type)
mat   
#     gt1 gt2 gt3 gt4
#   A  74  60  61  64
#   C  55  59  58  50
#   G  61  61  59  61
#   T  60  70  72  75
ggseqlogo(mat)

序列类型
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")

自定义配色：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)

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

同时绘制多个序列标志

ggseqlogo(seqs_dna, ncol = 4)
##直接绘图（不使用ggseqlogo），结果是一样的
ggplot()+geom_logo(seqs_dna)+theme_logo()+  
facet_wrap(~seq_group,ncol = 4,scales = "free_x")

R语言绘图包9--序列分析图的绘制ggseqlogo

1. 安装

2. 加载数据集

3. 可视化

猜你喜欢

热点阅读