R语言读入FASTA 蛋白质序列文件 & 生成SeqLogo

一、效果展示

1. FASTA文件蛋白质序列文件

FASTA和FASTQ是我们存储核苷酸序列信息（就是DNA序列）或者蛋白质序列信息最常使用的两种文本文件，虽然看起来名字有些古怪，但它们完全是纯文本文件（如同.txt）！名字的发音分别是fast-A和fast-Q。

【注】这里的序列、序列数据，指的其实就是表示DNA或者蛋白质的一条字符串。

FASTA文件主要由两个部分构成：序列头信息（有时包括一些其它的描述信息）和具体的序列数据。头信息独占一行，以大于号（>）开头作为识别标记，其中除了记录该条序列的名字之外，有时候还会接上其它的信息。紧接的下一行是具体的序列内容，直到另一行碰到另一个大于号（>）开头的新序列或者文件末尾。
下面给出一个FASTA文件的例子。

2. 数据可视化

二、使用R语言实现

编译器：RStudio

1、读入FASTA文件：R语言-Biostrings包

如果你要使用R对存储DNA序列的fasta文件进行操作，那么Biostrings是你的不二之选。

1.1 安装

Biostrings包和普通的R包安装方法不一样，用以下命令安装：

source(“http://bioconductor.org/biocLite.R”)
biocLite(“Biostrings”) #安装Biostrings

1.2 使用

seqdata = readBStringSet(filepath, format=“fasta”, nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=TRUE)

参数解释：

• nrec：最多从文件读入多少记录，负数则表示该项被忽略
• skip：一个非负整数，表示在开始读入记录之前要跳过的记录数
• seek.first.rec：TRUE或者FALSE（默认）。如果为TRUE，则一定是以“>”（对于FASTA）或“@”（对于FASYTQ）作为第一行。如果找不到这样的行，就会出现错误。值得一提的是，TRUE还可以用于解析带有FASTA数据的CFF3文件。
• use.names：返回的向量是否应该被命名，对于FASTA，名称取自描述行；对于FASTQ，它们取自记录序列ID。删除名称可以帮助减少内存占用，例如：对于包含数百万个记录的FASTQ文件。

2. 绘制seq logo图：R语言-ggseqlogo包

2.1 安装

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

2.2 读数据

fasta_input <- read.table(fasta, header=F, row.names=NULL)
fasta_input <- as.vector(fasta_input$V1)

2.3 可视化

ggseqlogo(seqs_dna$MA0001.1)

2.4 自定义

ggseqlogo(data, facet = “wrap”, scales = “free_x”, ncol = NULL,
nrow = NULL, …)

参数：

• data：序列的特征向量或命名的序列列表，所有序列必须具有相同的宽度
• facet：刻面类型，可以是 ‘wrap’ 或 ‘grid’
• scales：刻面尺度，请参阅 facet_wrap
• ncol：列数，仅适用于facet=‘wrap’，请参阅facet_wrap
• nrow：同 ncol
• method：高度方法，可以是"bits" or “probability”（默认值：“bits”）
• seq_type：序列类型，可以是“auto”，“aa”，“dna”，“rna”或“other”之一（默认：“auto”，自动猜测序列类型）
• namespace：命名空间，用于自定义命名空间的单个字母的字符向量。可以是字母数字，包括希腊字符。
• font：字体名称。请参阅list_fonts可用字体。
• stack_width：字母堆栈的宽度在0和1之间（默认值：0.95）
• rev_stack_order：如果TRUE，字母堆栈的顺序颠倒，默认值：FALSE
• col_scheme：应用于序列标志的配色方案。请参阅list_col_schemes可用字体。（默认：“自动”，根据颜色方案自动选取seq_type）。还可以传递使用该make_col_scheme函数创建的自定义颜色方案对象
• low_col，high_col：如果使用定量颜色方案，则渐变的低端和高端的颜色（默认：“黑色”和“黄色”）。
• na_col：颜色方案中缺少字母的颜色（默认值：“grey20”）
• ……

3. Demo

#加载包
library(ggplot2)
library(ggseqlogo)
library(Biostrings)

fasta_input = readBStringSet("/Users/wyc/Documents/rcode/graduation_project/ecoli_test.fasta", format="fasta", nrec=-1L, skip=0L, seek.first.rec=FALSE, use.names=FALSE)
print(fasta_input)
fasta_input <- as.vector(seqdata)
print(fasta_input)
ggseqlogo(fasta_input, seq_type="dna")

控制台输出：

数据可视化：

三、使用在线工具Weblogo实现

WebLogo（http://weblogo.threeplusone.com）绘制seqlogo的老牌在线工具。相比于在R上绘制seqlogo图，网页版在线工具更加轻松容易。但同时也存在一定的局限性，就是不适合分析大批量数据。

1. 上传文件或者输入序列

2. 自定参数（可选）

3. 生成logo

四、扩展——Two Sample Logo

http://www.twosamplelogo.org/cgi-bin/tsl/tsl.cgi
Two Sample Logo是一个基于Web的应用程序，用于计算和可视化两组对齐的氨基酸或核苷酸样本之间的差异。计算对齐的序列组中每个位置的每个残基的统计学显着性，其中零假设是根据阳性和阴性样品中的相同分布产生残基。两个样品标志物可用于确定各种活性位点，蛋白质修饰位点周围的统计学上显着的残基，或发现共享相同序列基序的两组序列之间的差异。

该软件支持两种类型的图形表示：（i）对于每个残基符号使用相同大小绘制统计上显着的残差，（ii）使用与两个样本之间的差异成比例的符号的大小来绘制统计上显着的符号。 。残基分为两组：（i）富含阳性样品，和（ii）贫样阳性样品。在所有类型的表示中，可以使用各种颜色方案来绘制符号，并且可以将共享的残差添加到图中以便可视化图案本身。使用二项分布（更准确但更慢的选项）或t检验（不太准确，但明显更快）计算p值。

两个样本徽标是使用Ruby编程语言创建的，基于WebLogo代码。

参考：
从零开始完整学习全基因组测序（WGS）数据分析：第2节 FASTA和FASTQ

把fasta序列读入到R里面去
R包Biostrings的使用( 转载）
Biostrings 处理fasta文件

R包ggseqlogo 绘制seq logo图
ggseqlogo

Seq logo 在线绘制工具——Weblogo

Two Sample Logo

其他资料：
Biostrings.pdf（英文的，我是看不下午去的，码住备用吧）
w3school R语言教程（R语言中文教程）
TwoSampleLogo.pdf