FASTA与FASTQ格式详解：从序列结构到质量值解析与实战

做了几年生信分析每天打交道最多的不是代码而是各种格式的文件。fastq和fasta这两个名字几乎所有和测序数据沾边的人都听过但说实话很多人对它们的理解停留在“一个是带质量值的一个是不带的”这种层面一旦遇到格式解析报错、质量值乱码、序列莫名多了一截这类问题还是容易一头雾水。这篇文章我想把fastq和fasta这两种最基础的序列文件格式彻底讲透。不光是格式规则还包括它们背后的设计逻辑、常见工具的解析方式、以及我实际项目中踩过的坑。如果你是刚入行的生信工程师、准备做NGS数据分析的学生或者只是偶尔需要处理测序数据的科研人员这篇文章应该能帮你省下不少排查问题的时间。1. 先搞清楚这两种格式到底解决了什么问题1.1 它们是什么能做什么fasta和fastq本质上是两种“纯文本约定的序列存储格式”。简单说就是用特定的符号规则把DNA、RNA或蛋白质序列以一种计算机容易读取、人类也能直接看懂的纯文本形式保存下来。fasta格式诞生得最早大概是1988年前后由Pearson和Lipman提出最初是为了配合FASTP和FASTA这些序列比对程序使用。它的核心思路极其简洁一行描述信息加若干行序列。描述行用大于号“”开头后面跟着序列的标识符和注释信息序列行就是纯粹的碱基或氨基酸字母没有特殊符号没有空格。fastq格式则要晚一些它是在高通量测序技术普及之后逐渐成为行业标准的。因为测序仪吐出来的原始数据里除了碱基序列本身还有一个非常重要但fasta无法表达的信息——每个碱基的测序质量。于是fastq在fasta的基础上额外加了两行用“”开头记录标识符用“”号作分隔最后一行用ASCII字符编码质量分数。我的理解是fasta是“序列的档案”fastq是“测序的原始凭证”。两者解决的问题不同不能简单地说谁比谁高级。做数据库存储、参考基因组比对这种对质量值不敏感的场景用fasta就够做变异检测、表达定量这类必须考虑测序可靠性的分析就必须保留fastq。1.2 这两类格式直接影响下游分析的正确性很多项目最终的结论出错根源往往不是算法选得不对而是输入的文件格式就没搞对。我见过有人把fastq文件直接改个扩展名当fasta用结果程序报错说“非法字符”也见过有人从NCBI下载数据时选错了SRA转fastq的参数导致质量值编码版本不一致后续质控报告里所有Q20/Q30比率的统计全部失真。所以理解fasta和fastq的格式细节不是“多此一举”的知识储备而是每一个要做序列分析的人都绕不开的基本功。接下来我分两部分详细拆解这两种格式的构成和内在逻辑。2. FASTA格式全拆解从一行“”开始2.1 格式规则与约定fasta格式的基本结构非常直观就两个部分描述行也叫标题行以“”开头后面紧跟该序列的唯一标识符ID空格之后的部分可用来写物种、功能注释等辅助信息。序列行由连续的碱基或氨基酸字母组成可以占一行也可以按固定长度折成多行解析时程序会自动忽略换行符把多行拼接成一个完整序列。用一套最简单的示例chr1 人类一号染色体片段 ACGTACGTACGT ACGTACGTACGT ACGTACGTACGT这里“chr1”是序列ID“人类一号染色体片段”是注释。从“”到下一个“”之间的所有字符都被视为同一条序列。这里有一个关键细节序列行中不允许出现数字、空格和除特定符号以外的任何字符。核酸序列标准字母是A、C、G、T还有代表不确定碱基的简并碱基符号如N、R、Y、S、W、K、M、B、D、H、V蛋白质序列则是20种氨基酸的单字母缩写加上B、Z、X等特殊符号。比对文件中还会出现“-”gap或“.”分别表示缺失或空位。2.2 不同来源的FASTA有什么差异平时用到的fasta文件虽然格式规则一致但来源不同细节差别还挺大的参考基因组fasta通常一条染色体对应一个序列条目ID可能类似“chr1”“NC_000001.11”或者“1”。注释行里常见物种名、组装版本信息。这种文件的序列行常见60或80个字符换行也有用100字符甚至更长的。转录本或蛋白fasta来自NCBI RefSeq、Ensembl、UniProt等数据库ID规则各成体系比如“NM_001301717.2”“ENST00000456328.2”“sp|P12345|ALBU_HUMAN”等。解析时最好只取第一个连续字段作为ID不要硬套某种固定规则。自定义小fasta很多脚本工具、引物设计软件的输出ID往往很简单甚至只有一行序列没有注释。在写解析代码时我一般建议按“第一个大于号到第二个大于号之间”来分条目绝不假设序列一定占一行也不假设注释行一定以某个固定模式结尾。2.3 容易被忽略的细节与陷阱换行问题。很多新手写Python解析fasta时直接用for line遍历文件然后判断line[0] 却忘了把每行的换行符strip掉结果序列末尾混入“\n”后续比对或者翻译全文报错。正确的做法是用rstrip而不是strip因为序列内部可能有合法的空字符虽然建议避免。ID重复问题。同一个fasta文件里出现两个相同ID的序列这在有些工具看来是合法但危险的。如果你的下游分析是按ID去重或索引重复ID会导致结果不可复现。我写代码时通常加一个查重步骤。空序列条目。早期GFF/GTF注释文件配套的fasta里偶尔会有“chrUn”这种没有任何序列的条目或者只包含“N”的序列。这种数据看着无害一旦涉及序列长度统计、k-mer计数就可能引入偏差。所以批量处理时建议统计每个条目的长度分布和GC含量提前发现异常。3. FASTQ格式全拆解四行记录一个碱基的完整信息3.1 四行结构的严格约定fastq每条read固定由四行组成EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG CTACGTACGTACGTACGTACGTACGTACGT IIIIIIIIIIIIIIIIIIIIIIIIIIIIII第一行以“”开头是read的标识符。这个ID和fasta的“”作用类似但fastq的标识符往往包含flowcell、lane、坐标、索引等测序仪信息字段多且长。需要注意第一行的ID和第三行的“”第三行通常可以是空ID也可以重复第一行ID视测序仪和工具而定。第二行碱基序列。长度就是read长度。第三行以“”开头作为分隔。有些文件这里重复一次ID有些直接就是一个“”解析时不能假设它一定有内容。第四行质量字符串。长度必须和第二行序列长度完全一致每个字符对应序列中对应位置碱基的质量编码。一个很常见的报错就是“sequence and quality lengths differ”说明某条read的第二行和第四行长度对不上。这种文件多半是拼接、裁剪或者传输过程中出了问题工具会直接拒绝读取。3.2 Phred质量分数与ASCII编码版本fastq的核心创新在于质量值。测序仪对每个碱基的判定都有一个错误概率我们要用一个数值来表示“这个碱基有多可靠”。业界默认用Phred质量分数公式是Q -10 × log10(P)其中P是碱基判定错误的概率。Q10表示错误率0.190%正确率Q20表示错误率0.0199%正确率Q30表示错误率0.00199.9%正确率。测序数据的质控报告里经常提到的Q20、Q30指的就是质量分数达到20或30以上的碱基占比。但fastq文件里第四行存的不是数字而是ASCII字符。因为如果直接存数字每个碱基至少占1到2个字节文件会膨胀用字母和符号来映射数字每个碱基质量值只占1个字符紧凑又便于快速读取。映射规则是把质量分数加上一个固定的偏移量再转成对应的ASCII字符。不同时期的测序平台用过的偏移量不一样这是很多历史遗留问题的根源。Phred33编码也叫Phred33ASCII码33-126现代Illumina设备1.8版本、Ion Torrent、以及绝大多数当前数据都用这个。质量值ASCII码值-33。比如字符“I”的ASCII码是73所以对应的Q 73 - 33 40。Phred64编码ASCII码64-126早期IlluminaGA II等1.3-1.7版本用过。质量值ASCII码值-64。这个编码下质量值的下限是0上限62左右。如果你拿Phred33的规则去解析Phred64的文件质量值会整体偏高约31读出来会比较“虚高”。Solexa编码更早期的Solexa平台用一套略有差异的映射公式是Q_solexa -10 × log10(P/(1-P))现在已经很少见基本只存在于一些考古级数据里。判断一个fastq是Phred33还是Phred64最土的办法是看质量字符串里有没有出现ASCII码小于64的符号。Phred33允许的字符包含大量标点符号如 !#$%()*,-./0123456789:;?ABCDEFGHI等等而Phred64编码里小于64的符号比如“!”到“?”都不合法。如果文件里出现了“!”、“#”、“$”这类符号基本可以断定不是Phred64。如果所有质量字符都是大写字母或某些符号且在到I之间大概率是Phred33。现代数据几乎统一用Phred33但处理旧数据时还是要多留个心眼。3.3 怎么快速查看和验证编码拿到一个陌生fastq文件我建议先不要直接丢给流程花几秒钟确认下面几个信息# 查看前几条read head -n 20 sample.fastq # 查看ASCII码值范围 head -n 400 sample.fastq | awk {if(NR%40){for(i1;ilength($0);i){print ord(substr($0,i,1))}}} | sort -n | uniqawk里没有内置ord函数可以改用FastQC或者一个小Python脚本import sys from collections import Counter with open(sys.argv[1], r) as f: cnt Counter() for i, line in enumerate(f): if i 400: break if i % 4 3: line line.rstrip(\n) cnt.update(ord(c) for c in line) print(sorted(cnt))对比输出结果的最小值最小值在33到73之间且包含低于64的符号通常是Phred33最小值在64以上且没有低于64的符号可能是Phred64。这是排查问题前“先看数据再动手”的好习惯能省去后面一整轮质控误判。4. FASTA与FASTQ的实战对比与转换操作4.1 核心差异对照做一个直接的对比方便抄作业维度FASTAFASTQ每序列行数不固定1行描述N行序列固定4行描述符号“”开头第一行“”第三行“”质量信息无第四行ASCII质量字符串主要来源参考基因组、数据库蛋白/核酸序列测序仪原始输出、比对前数据典型扩展名.fa、.fasta、.fna、.pep.fastq、.fq使用场景比对索引、序列搜索、存储组装结果质控、比对定量输入、变异检测输入文件大小相对较小通常大4倍以上多两行质量标识符是否允许折行允许序列换行每行固定内容不允许序列折行fastq是固定的四行结构所以相对而言更容易按行解析fasta因为序列可以折行解析反而要麻烦一点需要记录当前是描术行还是序列行的状态机。4.2 常用转换工具与命令在实际项目中fastq转fasta的需求非常常见比如做参考序列提取、做BLAST搜索输入或者某些只接受fasta的工具。推荐几个稳定可靠的工具。seqtk这是我用最多的工具轻量、快、不依赖库# fastq转fasta seqtk seq -A in.fastq out.fasta # 同时过滤低质量read质量值低于20的碱基占比超过10%的read不要 seqtk seq -A -q20 -n 0.1 in.fastq out.clean.fasta参数解释-A输出fasta格式-q20表示质量值阈值的Phred分数-n 0.1表示允许的N碱基最大比例或者叫允许低质量碱基的最大比例。BioPython适合需要在Python里进一步做操作的场景from Bio import SeqIO with open(input.fastq, r) as fq, open(output.fasta, w) as fa: for record in SeqIO.parse(fq, fastq): seq_str str(record.seq) fa.write(f{record.id}\n{seq_str}\n)这里要注意record.description可能包含多余空格和注释直接用record.id更干净。samtools适合BAM文件与fastq/fasta互转的场景不直接做fastq到fasta的转换但如果你手里是BAM格式想提取参考序列或者其他信息时它很有用。另外还有一个容易被忽略的点不要用Excel或文本编辑器直接打开大型fastq/fasta文件。这类文件动辄几个GBExcel根本打不开即便强行用某些支持打开大文件的编辑器打开也容易因为内存占用直接卡死电脑。处理这些文件永远用命令行工具。这和很多其他场景下“Excel打不开某些格式文件”的问题很像——不是文件损坏而是工具根本不适配但生物信息学数据一旦用Excel打开保存过文件就可能真的废了Excel会自动把某些看起来像数字的字符串改写成科学计数法破坏序列行内容。4.3 转换时需要避开的坑质量信息不可逆。fastq转fasta后质量值信息就永久丢失了。如果还想做下游质量过滤、去除接头后的二次质控必须保留原始fastq文件。所以转换前一定确认好真的不需要质量值了吗ID冲突风险。有些fastq里的注释行包含空格如果直接取整个第一行作为fasta的描述行后续某些工具可能因为ID含空格而解析失败。我建议在fastq转fasta时只保留后第一个连续字段作为ID其他的注释信息可要可不要。重复ID去重。尤其在做转录组数据时一个基因可能对应多条转录本有些工具生成的fastq里ID会有重复。如果后续要用CD-HIT或者其他聚类工具重复ID会影响聚类结果。可以用seqkit rmdup这类工具做去重或者直接在转换时加上随机后缀。压缩格式。常见的.fastq.gz、.fasta.gz是用gzip压缩的。命令行工具大多支持直接读取gz文件但如果自己写脚本处理gz一定要用gzip.open而不是普通open否则读出来全是乱码。5. 我踩过的坑和排查技巧实录5.1 高频问题速查与解决方法日常和这两种格式打交道遇到的无非就是下面这几类问题我整理成了一张速查表现象可能原因快速排查与解决程序报错“sequence and quality lengths differ”fastq第二行与第四行长度不一致文件传输中断或工具处理出错用awk检查是否有行数不是4倍数的记录重新下载或重新生成数据质量值普遍偏高或偏低质控曲线怪异Phred编码识别错误Phred64和Phred33混用看质量字符串ASCII最小值小于64视为Phred33大于等于64视为Phred64用FastQC自动判断序列中出现非法字符如数字、*文件损坏或者本来就不是标准fasta/fastqgrep -v file.fasta解析fasta时序列缺失或乱序序列折行导致状态机判断错误或者文件里有非ASCII隐藏字符检查是否有\rWindows换行用dos2unix或sed清理保证状态机正确处理“”ID重复导致下游比对失败同一ID出现多次与参考数据库不匹配用awk统计重复ID确认后加后缀或去重fastq文件第一行开头序列行也出现某些平台read ID后跟的注释列里包含符号按四行结构严格解析不要靠行首字符判断压缩包gz损坏下载不完全或磁盘空间不够用gzip -t测试完整性重新下载文件名或扩展名看起来对但工具不认文件内容与扩展名不匹配比如把Excel输出改名成fasta用file命令查看真实类型重新导出为纯文本这些问题里最关键的一个习惯是拿到文件先做“体检”不要着急跑主流程。5.2 一个真实排查案例质量值异常虚高有一次我在处理一批公开数据时发现质控报告里Q30比例高得离谱几乎每个碱基都是Q40。经验告诉我这不正常正常的二代测序数据虽然好但不可能高成这样。我临时写了个小脚本统计质量字符串的ASCII码范围果然最小值都在64以上几乎全是“h”到“~”这种字符。判断这是Phred64编码的旧数据而我的下游工具按Phred33去解析了所以把质量值整体抬高了31左右。解决办法很简单用seqtk或FastQC的编码修正参数重新解析# 显式指定Phred64 seqtk seq -Q64 -A in.fastq out.fasta或者直接在FastQC里指定编码类型后再看质控结果。恢复正确编码后Q30比例回到正常范围后续变异检测的假阳性率也明显改善。这个案例想说明的是格式问题经常不会让程序直接崩溃而是悄无声息地污染后续结果。如果对数据来源和编码版本心里没底一定要在分析开始前做好验证。5.3 文件处理习惯与工作流建议作为还算有经验的从业者我总结了一套处理fasta/fastq文件的固定工作流分享给大家参考第一步确认编码与完整性。对fastq执行head查看若干条read用小脚本看ASCII范围用gzip -t验证压缩文件完整性。第二步确认ID规则。查看描述行的字段分隔符确定ID是取前空格还是前竖线这一步决定了后续所有按ID合并、去重、过滤的操作是否正确。第三步制定格式规范。下游分析前统一把fastq转成Phred33编码把fasta统一成单行序列或固定80字符折行避免因为“不同工具对折行容忍度不同”而出问题。第四步数据备份与保留。fastq原始数据是分析的唯一凭证无论中间产出多少fasta、bam、vcf都建议把原始fastq或它的md5校验和保存好任何下游结果出问题都能回溯。第五步自动化检查写入流程。用snakemake或nextflow这类流程管理器的话建议在规则里加上对输入文件格式的检查步骤比如用seqkit stats输出序列数、总长度、GC含量用断言保证reads数不小于某个阈值这样流程能更早感知到输入异常。6. 工具选型分享这几个命令足够覆盖日常80%的需求说到处理fasta和fastq没有必要依赖很重的可视化软件命令行工具足够高效可靠。我个人平时固定使用下面几个组合起来能覆盖绝大多数日常操作seqtk格式转换、随机抽样、序列清洗、长度过滤。轻量高效几乎每个生信环境都会装。seqkit处理fasta/fastq的瑞士军刀。除了转换它还能统计每条序列的长度分布、GC含量、去重、取子集、滑窗过滤。输出统计信息非常方便。samtools处理SAM/BAM/CRAM也支持从reference fasta中提取特定区域序列。BioPython做更复杂的自定义解析比如同时读取fasta和注释表做联动处理或者批量修改序列ID。Python的pysam当你需要以编程方式处理比对后数据并输出fastq/fasta子集时pysam比BioPython更底层更高效。命令示例# 用seqkit统计基本信息 seqkit stats input.fastq # 按长度过滤取序列长度大于等于100的read seqkit seq -m 100 input.fastq input.filtered.fastq # 随机富集1万条read适合做快速测试 seqkit sample -n 10000 input.fastq subset.fastq工具的使用技巧很多但万变不离其宗先弄清输入文件的准确格式信息再选合适的工具。不要一上来就套一个万能脚本否则很容易把格式错误也一起“处理”掉。7. 最后记录一些个人体会接触生信这几年我最大的感受是格式解析这种看似基础的工作处理不好是整个分析链路上最耗时的坑。有些项目卡了整整两天最后发现只是因为fastq文件里夹了几行Windows换行符\r工具把这些\r当成了序列的一部分导致所有reads长度增加1个碱基下游比对全部报错。所以我建议所有做测序数据处理的朋友都把这个习惯刻进脑子里拿到任何fasta/fastq文件永远先检查三件事——行数是否符合预期格式、序列长度是否一致、质量值ASCII范围是否符合编码版本。这三个检查用不了两分钟但能帮你避免后续无数的“灵异事件”。如果你手头正好有一批数据不知道怎么判断编码或者写解析脚本时总在处理边界情况欢迎把具体报错发出来一起讨论。格式解析本质上是对“数据基础”的理解基础越扎实后面的分析越不容易翻车。