参考来源
生信常用数据格式: FASTA 格式
20160406 FASTA 与 FASTQ格式详解
FASTA Format for Nucleotide Sequences (nih.gov)
FASTA
1.1 历史
1985年，David J. Lipman 和 William R. Pearson 发表了一篇文章，Rapid and Sensitive Protein Similarity Searches。用C语言开发的程序: FASTP，进行蛋白质相似性搜索。1988年，两人发表Improved tools for biological sequence comparison
，FASTP 功能从原来的搜索蛋白质相似性 (FAST Protein) 扩展到了 (FAST ALL) ，可以完成所有（蛋白质和核苷酸）序列相似性搜索。
1.2 简介
- FASTA是一种基于ASCII 码的文本格式，可以存储核苷酸序列或氨基酸序列。它存储多个序列记录。
- 每一个序列数据以单行描述开始（必须单行），后跟紧跟一行或多行序列数据。也就是说，序列名称和注释位于序列之前。
- 是BLAST组织数据的基本格式，生物的基因组序列数据也是用FASTA格式存储的。
1.3 组成
描述行
- 描述行 (The description line, Defline, Header or Identifier line): 以一个大于号(">")开头，代表唯一的SeqID (sequence identifier)
每个序列的 SeqID 必须是唯一的，且不得包含任何空格。NCBI数据库工作人员在处理提交的材料时，将用Accession number取代 SeqID。
>SeqABCD
SeqID 后必须有序列来源的信息，格式为[modifier=text]。请勿在"="前空格。至少应包括该生物的学名。可用的modifier及其格式的完整形式为:
›SeqABCD [organism=Mus musculus] [strain=C57BL/6]
FASTA 定义行的最后一个可选部分是序列标题，它将用作DEFINITION 字段。标题应包含序列的简要描述。核苷酸和蛋白质的标题有一种首选格式。在处理过程中，数据库工作人员会将所提供的标题更改为适当的格式。
SeqABCD [organism=Mus musculus] [strain=C57BL/6] Mus musculus neuropilin 1 (Nrp1) mRNA, complete cds.
请注意，在所有情况下，FASTA 定义行不得包含任何硬回车。所有信息必须在一行文本中。如果在导入 FASTA 序列时遇到困难，请仔细检查编辑软件是否在 FASTA 定义行中添加了回车。
›SeqABCD [organism=Mus musculus] [strain=C57BL/6] Mus musculus neuropilin 1 (Nrp1) mRNA, complete cds.
请注意，在所有情况下，FASTA 定义行不得包含任何硬回车。所有信息必须在一行文本中。如果在导入 FASTA 序列时遇到困难，请仔细检查编辑软件是否在 FASTA 定义行中添加了回车。
格式正确的核苷酸序列 FASTA 定义行示例：
›Seq1 [organism=Streptomyces lavendulae] [strain=456A] Streptomyces lavendulae strain 456A mitomycin radical oxidase (mcrA) gene, complete cds.
›ABCD [organism=Plasmodium falciparum] [isolate=ABCD] Plasmodium falciparum isolate ABCD merozoite surface protein 2 (msp2) gene, partial cds.
›DNA.new [organism=Homo sapiens] [chromosome=17] [map=17q21] [moltype=mRNA] Homo sapiens breast and ovarian cancer susceptibility protein (BRCA1) mRNA, complete cds.
序列行
- 序列行 (Sequence Line)：一行或多行的核苷酸序列或肽序列，其中碱基对或氨基酸使用单字母代码表示。
定义行之后的一行开始核苷酸序列。与 FASTA 定义行不同，核苷酸序列可以包含回车。建议每行序列不超过 80 个字符。请仅在核苷酸序列中使用 IUPAC 定义的字符。
序列应使用标准的 IUB/IUPAC 氨基酸和核酸代码表示，以下情况除外：
接受小写字母，并将其映射为大写字母；
可用单个连字符或破折号表示长度不确定的空格；
查询序列中的任何数字都应删除或用适当的字母代码代替（例如，N 表示未知核酸残基，X 表示未知氨基酸残基）。
实际操作中，程序处理时自动去掉空格和换行符，把序列读成1行再处理
1.4 编码方式
Nucleic Acid Codes:
A --> adenosine M --> A C (amino)
C --> cytidine S --> G C (strong)
G --> guanine W --> A T (weak)
T --> thymidine B --> G T C
U --> uridine D --> G A T
R --> G A (purine) H --> A C T
Y --> T C (pyrimidine) V --> G C A
K --> G T (keto) N --> A G C T (any)
- gap of indeterminate length
Accepted Amino Acid Codes:
A ALA alanine P PRO proline
B ASX aspartate or asparagine Q GLN glutamine
C CYS cystine R ARG arginine
D ASP aspartate S SER serine
E GLU glutamate T THR threonine
F PHE phenylalanine U selenocysteine
G GLY glycine V VAL valine
H HIS histidine W TRP tryptophan
I ILE isoleucine Y TYR tyrosine
K LYS lysine Z GLX glutamate or glutamine
L LEU leucine X any
M MET methionine * translation stop
N ASN asparagine - gap of indeterminate length
1.5 示例
>sp|P69905|HBA_HUMAN Hemoglobin subunit alpha OS=Homo sapiens GN=HBA1
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR
从UniRef数据库中下载的人类血红蛋白α亚基的序列。
第1行中，P69905是这个序列在UniRef中的编号，HBA_HUMAN是这个序列的简称，Hemoglobin subunit alpha是全称，OS=Homo sapiens是物种，GN=HBA1是指基因的名字为HBA1.
后面的内容就是这个HBA1基因对应的蛋白的序列。
人类血红蛋白a亚基对应的mRNA序列，从NCBI RefSeq数据库中下载。
>gi|13650073|gb|AF349571.1| Homo sapiens hemoglobin alpha-1 globin chain (HBA1) mRNA, complete cds
CCCACAGACTCAGAGAGAACCCACCATGGTGCTGTCTCCTGACGACAAGACCAACGTCAAGGCCGCCTGG
GGTAAGGTCGGCGCGCACGCTGGCGAGTATGGTGCGGAGGCCCTGGAGAGGATGTTCCTGTCCTTCCCCA
CCACCAAGACCTACTTCCCGCACTTCGACCTGAGCCACGGCTCTGCCCAGGTTAAGGGCCACGGCAAGAA
GGTGGCCGACGCGCTGACCAACGCCGTGGCGCACGTGGACGACATGCCCAACGCGCTGTCCGCCCTGAGC
GACCTGCACGCGCACAAGCTTCGGGTGGACCCGGTCAACTTCAAGCTCCTAAGCCACTGCCTGCTGGTGA
CCCTGGCCGCCCACCTCCCCGCCGAGTTCACCCCTGCGGTGCACGCCTCCCTGGACAAGTTCCTGGCTTC
TGTGAGCACCGTGCTGACCTCCAAATACCGTTAAGCTGGAGCCTCGGTGGCCATGCTTCTTGCCCCTTTG
G
对于第1行，所有来自于NCBI的序列都有一个gi号，就是数据库中的流水号，具有唯一性。gb|AF349571.1是genebank编号的信息。后面就是对序列的一个通俗的描述。这里使用的是mRNA，包含完整的CDS（Coding sequence）区域。
为了保证数据的统一性，mRNA的序列还是用T来表示而不是U。因为U只是在RNA中替换了原来的T，所以为了下游的方便分析处理，无论RNA序列还是DNA序列都是使用T而不是U。
相同的序列被不同的人处理之后、甚至是在不同的网站上或者数据库中头信息都不尽相同，以下的情况都是可能存在的。一般用一个空格把头信息分为两个部分：第一部分是序列名字，它和大于号（>）紧接在一起；第二部分是注释信息。
>ENSMUSG00000020122|ENSMUST00000125984
> ENSMUSG00000020122|ENSMUST00000125984
>ENSMUSG00000020122|ENSMUST00000125984|epidermal growth factor receptor
>ENSMUSG00000020122|ENSMUST00000125984|Egfr
>ENSMUSG00000020122|ENSMUST00000125984|11|ENSFM00410000138465
1.6 优势
FASTA格式非常的简单和容易被理解，降低了序列操纵和分析的难度，易于传播和交流。
FASTQ
简介
fastq格式中的q为quality,一般用来存储原始测序数据及其质量，扩展名一般为fastq或者fq。illumina，BGISEQ一般是双末端测序，一对文件，命名为_R1.fq.gz与_R2.fq.gz。
2.1 命名
Illumina测序仪下机FASTQ命名为（NextSeq CN500下机数据为bcl格式，经过bcl2fastq转化后名称类似），例如：
Samplexx_S53_L002_R1_001.fastq.gz
Samplexx：样本名，与上机时在sampleSheet中填写的一致； S53：S后跟的数字与样本在sampleSheet中的顺序一致，从1开始； L002：lane编号； R1：R1表示read1，R2表示read2。R1和R2为paired end reads。 001：001，通常为001；
Undetermined_S0_L001_R1_001.fastq.gz：存储index不匹配的reads
2.2 示例
示例1
@SIM:1:FCX:1:15:6329:1045:GATTACT+GTCTTAAC 1:N:0:ATCCGA
TCGCACTCAACGCCCTGCATATGACAAGACAGAATC
+
<>;##=><9=AAAAAAAAAA9#:<#<;<<<????#=
第一行，Sequence identifier 序列标识以及相关的描述信息
官网给的格式解释如下：
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos>:<UMI> <read>:<is filtered>:<control number>:<index>
@SIM:1:FCX:1:15:6329:1045:GATTACT+GTCTTAAC 1:N:0:ATCCGA
@, seqID开始符
SIM，测序仪的ID号
1，设备的run number
FCX，flow cell的ID号
1，lane号
15，tile号（tile为flow cell上最小单位，测序时每测一个碱基，照相一次）
6329，flow cell中簇位置的X坐标
1045，flow cell中簇位置的Y坐标
GATTACT+GTCTTAAC 1，当sampleSheet存在UMI时该项存在；为Read1的UMI序列+Read2的UMI序列信息
1，1 表示 single read 2 表示 paired end
N，是否过滤，Y表示被过滤，否则为N
0，0表示不存在控制位
ATCCGA，index序列
image.png
第二行，序列信息，例如
TCGCACTCAACGCCCTGCATATGACAAGACAGAATC
第三行，Quality score identifier line (consisting only of a +)
以“+”开头，为节省存储空间什么也不加
旧版的FASTQ文件中会直接重复第一行的信息
第四行，Quality score，测序质量值
描述第二列中每个碱基的可靠程度，用ASCII码表示，我们平时长听到的Q20，Q30即为该字符对应的值，例如
<>;##=><9=AAAAAAAAAA9#:<#<;<<<????#=
示例二
Illumina平台测序的真实数据，包含1条reads的信息。
@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,
第1行主要储存序列测序时的坐标等信息
@ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
@ 开始的标记符号
ST-E00126:128:HJFLHCCXX 测序仪唯一的设备名称
2 lane的编号
1101 tail的坐标
7405 在tail中的X坐标
1133 在tail中的Y坐标
第2行是测序得到的序列信息，一般用ATCGN来表示。
第3行以“+”开始，一般是空的。
第4行储存的是质量信息，与第2行的碱基序列一一对应，其中的每一个符号对应的ASCII值成为phred值，可以简单理解为对应位置碱基的质量值，越大说明测序的质量越好。
FASTQ质量值计算方法
测序仪进行测序时，会自动根据荧光信号的强弱给出一个参考的测序错误概率（error probility，P）
问题与解决
- P值肯定是越小越好。1%储存成0.01是不高效的，因为1个碱基的信息，占用了至少4个字符。
将P取log10之后再乘以-10，得到的结果为Q。
比如，P=1%，那么对应的Q=-10*log10（0.01）=20
- 直接使用Q值表示质量值，数字直接连起来得加分隔符浪费空间。
ASCII码前0到32个为非可见字符，如空格等，所以需要Q值加上质量体系值（33或者64）
ASCII
把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基。
如Q=20，Phred = 20 + 33 = 53，对应的符号是”5”
反推
第一个碱基T对应的碱基质量ASCII码是<，查询ASCII码表中<对应的ASCII值为60，如果当前测序仪使用的质量体系为Phred33，则T对应的碱基质量值Q为27（60-33），可进一步推算出Q = -10log(p_error)中p_error。
Phred
Solexa
Illumina
不同测序仪的不同Phred值对应的ASCII表
质量体系
获取与转换
使用sratools工具中的prefetch或者fastq-dump软件都可以下载fastq文件。
利用fastq-dump文件可以将sra文件直接转换为fastq格式，注意，如果是illumina的双末端测序，需要添加 --split-files选项，如果需要压缩格式，需要添加 --gzip选项。最终会生成SRR8651554_1.fastq.gz，SRR8651554_2.fastq.gz两个文件。
fastq-dump --split --gzip ~/ncbi/public/sra/SRR8651554.sra
目前绝大部分的软件都可以直接处理压缩格式，因此一般的fastq格式都是压缩格式呈现的，扩展名为fq.gz，压缩gzip或者解压缩 gunzip
网友评论