FASTA&FASTQ格式

2023-11-20 本文已影响0人孟令君

参考来源
生信常用数据格式: 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 定义行示例：

›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