前言
刚刚踏入生信圈的你，有没有因为各种数据格式而困惑过？比如基因组Fasta序列，又或者二代测序的Fastq数据。这些数据有哪些特点？又有哪些日常操作是生信小白必须要掌握的？今天，就请跟随小编一起，了解一下这些数据格式的前世今生吧！
关卡·1
已知某细菌的protein序列(ref.protein.fasta)，现给定一组未知身份的protein序列(query.protein.fasta)，想要研究给定序列与该细菌的哪些protein序列相同或相似？
那么，我们需要具备哪些技能才可以顺利通关呢？
技能一：Fasta
1)查看方式
cat 或者less
2)数据格式
Fasta格式的数据，包括两个要素，一个是以大于号“>”开头的序列名称；另一个是序列主体，序列可以是一行，也可以是多行(习惯上按照固定宽度分隔)。
技能二：BLAST
在学习BLAST之前，我们先来介绍一下序列比对。序列比对一般分为全局比对(global alignment)和局部比对(local alignment)。全局比对是指将序列作为整体来进行比对，适用于整体相似度较高的序列间的比对。局部比对是指将序列之间所有相似性的局部片段都进行比对，适用于整体相似性较低、或者相似度未知的序列之间的比对。图1展示了全局比对和局部比对的含义。
图1 | 全局比对和局部比对示意图
了解序列比对的概念之后，就要郑重介绍一下BLAST啦。
Basic Local Alignment Search Tool(BLAST)是基于BLAST算法[1]的序列比对工具，由NCBI研发并维护，因其快速可靠、功能全面、使用方便而被广泛应用。表1总结了BLAST的程序及应用范围。BLAST分析可以通过网页https://blast.ncbi.nlm.nih.gov/Blast.cgi进行，当然，这并不是生信小白最好的选择。我们推荐使用命令行的形式来分析哦！
表1 | BLAST程序及应用范围
BLAST命令行分析
step1: 创建数据库
makeblastdb -in ref.protein.fasta -dbtype prot -out seqDatabase
makeblastdb: 程序名称，将给定的氨基酸序列创建数据库
-in: 输入文件，即需要创建数据库的氨基酸序列
-dbtype: 数据库类型，核酸对应nucl，氨基酸对应prot
-out: 所创建的数据库的前缀名称
step2: 选择合适的程序进行比对
blastp -query query.protein.fasta -db seqDatabase -outfmt 6 -num_threads 2 -out results.txt
blastp: 程序名称，氨基酸与氨基酸序列的比对
-query: 指定查询序列
-db: 指定查询序列与哪一个数据库进行比对
-outfmt: 输出格式
-num_threads:线程数
-out: 输出结果文件
输出结果见图二，包含有12列：
图2 | 比对结果
# 1. Query id：查询序列ID标识；
# 2. Subject id：比对上的目标序列ID标识；
# 3. identity：序列比对的一致性百分比；
# 4. alignment length：能够比对上的区域长度；
# 5. mismatches：比对区域的碱基错配数；
# 6. gap openings：比对区域的gap数目；
# 7. q_start：比对区域在查询序列(Query id)上的起始位点；
# 8. q_end：比对区域在查询序列(Query id)上的终止位点；
# 9. s_start：比对区域在目标序列(Subject id)上的起始位点；
# 10. s_end：比对区域在目标序列(Subject id)上的终止位点；
# 11. e-value: 对score值可靠性的评估，e-value值越小，score值越可靠；
# 12. bit score: 比对结果的bit score值，分越高，表示相似程度越大。
如果你掌握了以上技能，就可以轻松通过第一关啦！
关卡·2
已知某细菌的核酸序列(genome.fasta)，现给定一组fastq序列(query.R1.fastq和query.R2.fastq)，计算给定序列与该细菌核酸序列的回比率是多少？
那么，我们需要具备哪些技能才可以顺利通过第二关呢？
技能一：Fastq
Fastq格式一般是二代测序数据的存储方式，可以保存碱基及其质量得分信息。
1)查看方式
head -n4 query.R1.fastq
2)数据格式
第一行：以@开头，之后为序列的标识符以及描述信息
第二行：序列信息，以ATCGN表示，N代表未知碱基
第三行：+，之后可以再次加上序列的标识及描述信息(可选)
第四行：碱基质量得分信息，与第二行的序列相对应，长度必须与第二行相同
技能二：BWA
随着测序技术不断的进步，各种比对工具应运而生。这里我们介绍一款适用于短序列比对的软件：BWA[2]。
BWA命令行分析
step1:为参考基因组创建索引
bwa index genome.fasta
step2: 将短序列回比到参考基因组中
bwa mem genome.fasta query.R1.fastq query.R2.fastq > result.sam
技能三：SAM和BAM
小伙伴们可能已经发现了，通过比对软件分析后，我们会得到以sam结尾的文件，此时sam文件的解读就显得尤为重要啦！SAM的全称为The sequence alignment/map format，是用来存储测序reads与参考基因组比对结果信息的一种文件格式。SAM格式的诞生，规范了不同比对软件产生的结果，更有利于下游的分析。不过SAM文件一般较大，所以最好通过samtools[3]将其转换成二进制的BAM文件。
1)查看方式
SAM文件可以直接用less查看；
SAM格式转换成BAM格式：
samtools view -bS result.sam > result.bam
BAM文件可以通过samtools view result.bam |less -S 来查看
2)数据格式
SAM格式由两部分构成，分别是标头信息(header)和比对信息(record)。
Header主要存放基础信息，比如SAM格式的版本、比对序列的排列顺序、使用的比对程序等，均以@开头。如果小伙伴需要获取更详细的信息，可以访问https://www.samformat.info/sam-format-header哦！
Record每一行记录的是一条read与参考基因组的比对信息，以tab键分割，下面展示的是前11列的信息。
# 1. 测序read名称；
# 2. flag值，即位标识，每一个数字代表一种比对情况。flag值可以是多个数值相加的总和；我们可以直接在网站https://broadinstitute.github.io/picard/explain-flags.html上输入flag值，就可以迅速得到read的比对情况；
# 3. 参考基因组的序列名称；
# 4. 比对上的位置，从1 开始计数，如果没有比对上，则为0；
# 5. Mapping质量；
# 6. 简要比对信息表达式，使用数字加字母表示比对结果。数字表示碱基个数，字母表示比对情况；
# 7. 与之配对的read所比对上的参考序列的名称，如果没有mate序列，用*表示，如果存在mate序列，且比对到同一参考序列，用=表示；
# 8. mate 序列在参考序列上的位置，如果没有，则此处为0；
# 9. 插入片段的长度，一般R1端的为正，R2端的为负；
# 10. 序列信息；
# 11. 碱基质量，记录方法与Fastq相同。
通过两个关卡后，恭喜小伙伴们已经初步掌握了Fasta和Fastq文件格式及常用的处理方法。让我们相约下期，继续在生信的海洋里遨游吧！
参考文献
1. Altschul, S.F., Gish, W., Miller, W., Myers, E.W. & Lipman, D.J. (1990). Basic local alignment search tool. J. Mol. Biol. 215:403-410.
2. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform.[J]. Bioinformatics, 2010, 25(5):1754-1760.
3. Li H, Handsaker B, Wysoker A, et al. The Sequence Alignment/Map (SAM) Format and SAMtools[J]. Transplantation Proceedings, 2009, 19(1 Pt 2):1653-4.
END
本文系欧易生物原创
转载请注明本文转自欧易生物
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
