参考: 从零开始完整学习全基因组测序数据分析：第2节 FASTA和FASTQ
作者：碱基矿工

前言

在生物信息分析过程中，各种组学（基因、转录、蛋白等）各种数据库，不同类型的数据可能都有其特有的数据文件与特殊的存储格式。
而FASTA与FASTQ 正是其中之二。它们被用于存储核苷酸序列信息（即DNA序列）或蛋白质序列信息。

虽然被叫做FASTA与FASTQ，但实际上这两个文件均以类似.txt 的纯文本形式储存。

正文

FASTA

FASTA 来源于一款软件，该软件的名称就是“FASTA”，被用于进行序列比对。而名字中最后一个字母A，就是“assignment”。其最初是由William. R. Pearson 和 David. J. Lipman在1988年所编写，目的是用于生物序列数据的处理。

之后，FASTA 就被作为存储有顺序的序列数据的后缀。而序列，本质上就是用于表示核酸（碱基）或蛋白质（氨基酸）的字符串。文件后缀除了.fasta外，还有.fa 或压缩形式如.fa.gz。

必须强调的是序列是有顺序的。

因此在核酸序列的FASTA文件中，我们可以通过指定特定的某个数，便可知道某个DNA碱基在基因组上的具体位置。

FASTA 文件的构成

FASTA 文件主要由两个部分构成：1）序列头信息。2）具体的序列数据。

>AF466308.1 Dictyostelium discoideum ABC transporter AbcB5 (abcB5) gene, complete cds
TTTTAAATTTAAAACAAATACTATTATATTAATATATAAATAAATTTTATAATTACCTTTTTTTTTTGTT
'''
删去了中间的碱基信息
'''
GATATTTTGAATGATTTTCATCCAAAACTATTGAAATAATTTCTAAAGCTTTTTTAAAATGATCTAAAGA
T

上面的代码块内，便是一个标准的FASTA文件，从ncbi 的数据库中下载的原核生物的abcb5 基因序列。
虽然FASTA 文件信息构成很简单，但是往往序列信息的量非常的庞大。比如人类有30亿个基因，即文件至少由30亿个字符构成。而即使压缩起来，文件也有1Gb。
来源：https://www.ncbi.nlm.nih.gov/nuccore/AF466308.1?report=fasta

序列头信息

序列头信息指下面这段的内容，头信息独占一行，以大于号>开头作为识别标记。而紧接的下一行则是序列内容，直到碰到一行为>的新序列，或文件末尾。
头信息一般基因名称及其他一些注释信息。

>AF466308.1 Dictyostelium discoideum ABC transporter AbcB5 (abcB5) gene, complete cds

令人头疼的一个点是，序列头信息并没有被严格地限制。有时候，不同的人，不同的数据库，可能存储序列的FASTA 文件都不尽相同。

如下面五个不同的头信息：

>ENSMUSG00000020122|ENSMUST00000125984
> ENSMUSG00000020122|ENSMUST00000125984
>ENSMUSG00000020122|ENSMUST00000125984|epidermal growth factor receptor
>ENSMUSG00000020122|ENSMUST00000125984|Egfr
>ENSMUSG00000020122|ENSMUST00000125984|11|ENSFM00410000138465

“不成文”的规范

正因为头信息的混乱造成了较差的阅读性，甚至导致错误，因此也有了一些对头信息的不成文的规范：
1）第一部分是序列名字，与>相连。
2）第二部分用空格与序列名字相隔，表示注释信息，可以没有。

>gene_00284728 length=231;type=dna

具体的序列数据

序列数据由碱基（ATCG）或蛋白质氨基酸（各种氨基酸缩写）构成。一般每行60个字母，超过便换下一行。

额外注意

由于FASTA 文件的本质是文本文件，因此里面的内容是无法自动检验的。这也就意味着其中的序列可能存在内容重复的情况。对于某些标准的参考序列来说，这点不必担心，但其他一些序列，在使用前可以根据序列名称，进行简单的查重，以防万一。

FASTQ

FASTA 是经过人为排列（头信息与序列信息）的序列文件。而FASTQ 则储存的是测序仪产生的原始数据。它由测序的图像数据转换而来，也是文本文件。

FASTQ 文件构成

FASTQ 文件一般由四行内容构成。每四行为一个独立的单元，也叫read。

@DJB775P1:248:D0MDGACXX:7:1202:12362:49613
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
+
JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA
@DJB775P1:248:D0MDGACXX:7:1202:12782:49716
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
+
IIIIIIIIIIIIIIIHHHHHHFFFFFFEECCCCBCECCCCCCCCCCCCCCCC

FASTQ 是存储测序数据最普遍，也是公认的数据格式（还有一个是uBam格式）。

FASTA 是经过人为排列（头信息与序列信息）的序列文件。而FASTQ 则储存的是测序仪产生的原始数据。它由测序的图像数据转换而来，也是文本文件。

依照测序量或测序深度的不同，FASTQ文件大小也存在许多的差异，小到几M，大到百十G。一般文件的后缀为.fastq, .fq 或gz 压缩的.fq.gz。

第一行

第一行通常以@ 开头，是该段read 的名称，这段字符串是每一条read的唯一标识符，同一份FASTQ文件中不会重复出现，甚至不同的FASTQ文件里也不会有重复。

通常这段字符串由根据测序时的状态信息转换形成，中间不存在空格。

@DJB775P1:248:D0MDGACXX:7:1202:12362:49613

第二行

第二行为该段read 的测序序列，由A、T、C、G、N构成。N表示测序时无法被识别出来的碱基。

CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG

第三行

以+开头，一般什么信息都没有。

+

第四行

测序read 的质量值，和碱基信息一一匹配，表示每个测序碱基的可靠程度，用ASCII码表示，也非常重要。

IIIIIIIIIIIIIIIHHHHHHFFFFFFEECCCCBCECCCCCCCCCCCCCCCC

质量值

顾名思义，碱基质量值即用来衡量一个碱基好坏程度的一个数值。如果测得准，则其出错的概率就越低，质量值就越高；反之，其出错的概率就越高，质量值就越低。

质量值的计算公式如下：

Q = -10log(p_error)

质量值是测序错误率的对数（10为底数）乘以-10（并取整）。

若出错的概率为0.1，即p_error=0.1，则 Q = 10。

那为什么要用ASCII 存储质量值？

好看。

如果用数字表示，不仅由于数字的位数不同，不对齐，且需要通过分隔符号将每个质量值分开，也消耗内存。

但问题来了，小于33 的ASCII 码都是不可见的字符，如空格、占位符等。如何避开这些字符呢？

最简单的方法就是全部加上一个固定整数。你质量值是18， 加33就好了。

但起初科学家们都很随心所欲，你想从33 开始，我偏偏要从50开始。

但现在来说，一般广为应用的只有一种了：

即 质量数 1 对应在FASTQ 文件中显示为ASCII 码 33 对应的字符串!。

质量值体系检查

这里碱基矿工提供了一个shell 脚本以进行质量值体系的查看:

less $1 | head -n 1000 | awk '{if(NR%4==0) printf("%s",$0);}' \
| od -A n -t u1 -v \
| awk 'BEGIN{min=100;max=0;} \
  {for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END \
  {if(max<=126 && min<59) print "Phred33"; \
  else if(max>73 && min>=64) print "Phred64"; \
  else if(min>=59 && min<64 && max>73) print "Solexa64"; \
  else print "Unknown score encoding"; \
  print "( " min ", " max, ")";}'

可以通过将以上代码，储存在shell 文件中，用于FASTQ 类型文件的质量体系的检测。原理是截取前1000行数据，并抽取出质量值所在的行，接着分别计算，得出最大与最小的ASCII 值。

选定一个.fq 文件进行检验。

$ sh fq_qual_type.sh untreated.fq
Phred33
( 34, 67 )

ASCII 码转换

得到了目的FASTQ文件的质量体系，接着就可以进行数值的转换。当然如果你精通ASCII 码，直接将对应的文件内的质量值的字符串转换为其对应的ASCII码，接着减去使用的质量体系的下限值便可以得到该碱基的质量值。
如'J' 对应ASCII 码为74，而如果质量体系为Phred33，则通过减去33，则可知'J'对应的质量值为41。

如果嫌麻烦，也可以借助python 的ord()函数，将字符串内容输入进去，借助循环函数，转换后减33即可。