R语言绘图包系列：
- R语言绘图包01--优秀的拼图包patchwork
- R语言绘图包02--热图pheatmap
- R语言绘图包03--火山图EnhancedVolcano
- R语言绘图包04--GOplot：富集分析结果可视化
- R语言绘图包05--韦恩图的绘制：ggvenn和VennDiagram
- R语言绘图包06--基因表达相关性绘图corrplot
- R语言绘图包07--多集合可视化UpSetR
- R语言绘图包08--森林图的绘制：forestplot
序列分析图（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")
