R包ggseqlogo |繪制序列分析圖

前言

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測序分析?(重磅綜述：三萬字長文讀懂單細(xì)胞RNA測序分析的最佳實(shí)踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

簡介

在生物信息分析中，經(jīng)常會(huì)做序列分析圖（sequence logo），這里的序列指的是核苷酸(DNA/RNA鏈中)或氨基酸(在蛋白質(zhì)序列中)。sequence logo圖是用來可視化一段序列某個(gè)位點(diǎn)的保守性，據(jù)根提供的序列組展示位點(diǎn)信息。常用于描述序列特征，如DNA中的蛋白質(zhì)結(jié)合位點(diǎn)或蛋白質(zhì)中的功能單元。

實(shí)現(xiàn)以上可視化過程的工具有很多，本文介紹一個(gè)使用起來非常簡單，不拖泥帶水的R包ggseqlogo，只要你根據(jù)此包要求的數(shù)據(jù)格式上傳一堆DNA序列或者氨基酸序列，再根據(jù)現(xiàn)成的命令流程就能畫出logo圖。

安裝到作圖的代碼如下：

安裝

安裝方式有兩種

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

數(shù)據(jù)加載

ggseqlogo提供了測試數(shù)據(jù)ggseqlogo_sample。

#加載包
library(ggplot2)
library(ggseqlogo)
#加載數(shù)據(jù)
data(ggseqlogo_sample)

ggseqlogo_sample數(shù)據(jù)集是一個(gè)列表，里面包含了三個(gè)數(shù)據(jù)集：

• seqs_dna:12種轉(zhuǎn)錄因子的結(jié)合位點(diǎn)序列

• pfms_dna:四種轉(zhuǎn)錄因子的位置頻率矩陣

• seqs_aa:一組激動(dòng)酶底物磷酸化位點(diǎn)序列

#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]
## 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]

## $AKT1
## ? [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG"
## ? [4] "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"

外部數(shù)據(jù)讀入

也可以用自己的數(shù)據(jù)集，支持兩種格式，序列和矩陣。

# 長度為7的motif。每一行為一條序列，長度相同，每一列的堿基組成代表對應(yīng)位置的堿基偏好性。
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矩陣示例，每一列代表一個(gè)位置，及堿基在該位置的出現(xiàn)次數(shù)。共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)

可視化

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

ggseqlogo提供了一個(gè)直接繪圖的函數(shù)ggseqlogo(),這是一個(gè)包裝函數(shù)。下面命令結(jié)果同上面的。

ggseqlogo(seqs_dna$MA0001.1)

輸入格式

ggseqlogo支持以下幾種類型數(shù)據(jù)輸入：

• 序列

• 矩陣

下面是使用數(shù)據(jù)中的位置頻率矩陣生成的seqlogo

ggseqlogo(pfms_dna$MA0018.2)

方法

ggseqlogo通過method選項(xiàng)支持兩種序列標(biāo)志生成方法：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序列類型，默認(rèn)情況下ggseqlogo會(huì)自動(dòng)識別數(shù)據(jù)提供的序列類型，也可以通過seq_type選項(xiàng)直接指定序列類型。

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

自定義字母

通過namespace選項(xiàng)來定義自己想要的字母類型

#用數(shù)字來代替堿基
seqs_numeric <- chartr("ATGC", "1234", seqs_dna$MA0001.1)
ggseqlogo(seqs_numeric, method="prob", namespace=1:4)

配色

ggseqlogo可以使用col_scheme參數(shù)來設(shè)置配色方案，具體可參考?list_col_schemes

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

自定義配色

ggseqlogo提供函數(shù)make_col_scheme來自定義離散或者連續(xù)配色方案

離散配色

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)

連續(xù)配色

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

同時(shí)繪制多個(gè)序列標(biāo)志

ggseqlogo(seqs_dna, ncol = 4)

上述命令實(shí)際上等同于

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

自定義高度

通過創(chuàng)建矩陣可以生成每個(gè)標(biāo)志的高度，還可以有負(fù)值高度

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")

字體

可以通過font參數(shù)來設(shè)置字體，具體可參考?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))

注釋

注釋的話跟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")

嚴(yán)濤：浙江大學(xué)，作物遺傳育種在讀博士。愛鼓搗各種可視化軟件，愛折騰。點(diǎn)擊閱讀原文跳轉(zhuǎn)其博客。

更多數(shù)據(jù)可視化可參考在線工具：http://www.ehbio.com/ImageGP/

（點(diǎn)擊圖片直達(dá)網(wǎng)站手冊）

你可能還想看

• Illumina測序儀比較和各種測序應(yīng)用模式圖，助力了解高通量測序

• 咦！這樣畫基因結(jié)構(gòu)圖夠好看！
  

• MotifStack：多motif序列比較和可視化

• 在線RaxML構(gòu)建系統(tǒng)發(fā)育樹

• OrthoMCL鑒定物種同源基因 （安裝+使用）

往期精品(點(diǎn)擊圖片直達(dá)文字對應(yīng)教程)

后臺回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集