做生信的人几乎天天和 FASTA 打交道。无论是下载基因组、提取 CDS还是做引物设计第一步往往都是把 FASTA 文件里的序列读进内存然后做各种变换碱基计数、GC 含量、DNA 转录 RNA、RNA 翻译蛋白、反向互补。这些操作单个看都不难但拼在一起的时候很多人会在小地方翻车换行符没处理干净、把 DNA 密码子拿去查 RNA 密码子表、反向互补只互补忘了反向。这篇文章算是我把自己的常用工具箱整理了一遍。标题里的几个功能看起来基础但把它们写成规范、可复用的函数能省下后面大量重复劳动。适合刚入门的生信学习者也适合想把零散脚本整理成小模块的开发者。我会直接给出 Python 实现并把每个函数的设计理由、使用边界和容易踩的坑讲清楚。1. FASTA解析与序列提取把文件变成可计算的字符串1.1 FASTA 格式回顾为什么不能直接open().read()FASTA 格式本身不复杂规定也很简单每条序列以开头后面跟着序列标识符和可选的描述文字下一行起是序列本身序列可以换行可以全是小写也可以包含N等模糊碱基。比如chr1_gene001 partial sequence ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG GGCATTGTAATGGCTTAA问题在于“序列可以换行”。很多人第一次写读取脚本时直接用open(file).read()然后去掉所有之后的内容结果序列中间残留了换行符后面的计数、翻译全乱套。所以解析 FASTA 的第一步不是读字符串而是把每一条记录的注释行和序列行分开再把属于同一条记录的序列行拼接起来。1.2 一个通用的多序列读取函数我平时最常用的版本是返回一个字典键是序列 ID值是大写后的完整序列。这样方便后续按 ID 随便查。def read_fasta(filepath): seqs {} current_id None seq_lines [] with open(filepath, r, encodingutf-8) as f: for line in f: line line.strip() if not line: continue if line.startswith(): if current_id is not None: seqs[current_id] .join(seq_lines) current_id line[1:].split()[0] seq_lines [] else: seq_lines.append(line.upper()) if current_id is not None: seqs[current_id] .join(seq_lines) return seqs这个函数有四个细节值得注意第一line[1:].split()[0]只取第一个以空白字符分隔的 token。FASTA 的注释行经常写成id description如果直接存整行后面查序列 ID 时还得再做处理。只取第一个 token 是最省事的做法。第二line.strip()会去掉首尾空白包括换行符和回车符。这能兼容 Windows、Linux、Mac 三种换行风格也能处理行尾多余空格。如果不 strip拼接出来的序列可能夹杂\n和空格后面len(seq)统计长度都会出错。第三line.upper()统一转为大写。序列文件里经常出现小写碱基如果不统一大写后面做seq.count(G)时会把g漏掉。第四跳过空行。有些 FASTA 文件在记录之间有空行如果不跳过空行会触发else分支往seq_lines里拼进去一个空字符串虽然不会改变最终序列内容但属于无效操作。更危险的是文件末尾没有换行时循环结束后需要再检查一次current_id否则最后一条记录会被丢掉。如果文件里有多条序列且 ID 不重复返回字典非常方便。但要注意一点字典的键如果重复后面的记录会覆盖前面的记录。如果你确实拿到一个 ID 重复的 FASTA应该改用列表保存或者至少先打印一条警告。1.3 只要第一条序列的简易读取很多时候我们只需要验证一段序列能不能翻译并不关心文件里有多少条记录。这时候再走完整字典流程有点重。我习惯再写一个轻量函数遇到第二条就停下来只返回第一条序列。def read_first_seq(filepath): seq_lines [] started False with open(filepath, r, encodingutf-8) as f: for line in f: line line.strip() if not line: continue if line.startswith(): if started: break started True continue if started: seq_lines.append(line.upper()) return .join(seq_lines)这个函数没有 ID 返回只解决“我想赶紧拿到第一条序列跑下游分析”的需求。实际调试时可以少写很多临时文件操作。要注意的是如果 FASTA 文件里第一条序列本身跨了好些行这个函数会一直读到第二条注释行才开始 break所以拼接结果是完整的不会只拿第一行。2. 核苷酸计数与GC含量计算别小看这几个百分号2.1 用 Counter 一次完成四类碱基计数拿到序列字符串以后最简单的计数方式是连续调用str.count。比如seq.count(A)一次、seq.count(T)一次。但这样要扫描四遍序列写起来也啰嗦。Python 的collections.Counter可以一次性统计所有字符一行搞定。from collections import Counter def count_nucleotides(seq): seq seq.upper() counts Counter(seq) return counts用的时候直接c count_nucleotides(ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG) print(c) # Counter({G: 11, A: 9, T: 8, C: 8})但这里有个隐藏问题Counter会把所有字符都统计进去包括空格、N、-甚至换行符。所以传入之前必须保证序列是干净的最好在解析 FASTA 时就已经转成大写并且过滤掉非碱基字符。如果你要统计的序列可能含有模糊碱基建议还是写一个只统计ATCGU的版本def count_bases(seq): seq seq.upper() counts {A: 0, T: 0, U: 0, G: 0, C: 0} for base in seq: if base in counts: counts[base] 1 return counts这个版本的优点是输入 RNA 时U会单独统计输入 DNA 时U为 0。如果你在具体分析中不需要区分 DNA 和 RNA可以去掉U这一项。我通常把两个函数都放在工具模块里count_nucleotides做全量字符快速预览count_bases做严格的核酸计数。2.2 GC 含量计算分母到底用谁GC 含量指鸟嘌呤和胞嘧啶占所有碱基的比例单位是百分比。计算公式并不复杂GC% (G的数量 C的数量) / 总碱基数 * 100但“总碱基数”有分歧。如果序列里有N未知碱基、-缺口、甚至U分母到底要不要包含它们我在实际项目中通常提供两个参数def gc_content(seq, ignore_ambiguousTrue): seq seq.upper() if len(seq) 0: return 0.0 if ignore_ambiguous: # 只统计明确的 ATCG total seq.count(A) seq.count(T) seq.count(G) seq.count(C) gc seq.count(G) seq.count(C) else: total len(seq) gc seq.count(G) seq.count(C) if total 0: return 0.0 return gc / total * 100当ignore_ambiguousTrue时分母只包含 A、T、G、CN和-被排除。这样适合比对序列、组装序列因为未知碱基不该影响 GC 占比。当ignore_ambiguousFalse时分母就是序列总长度适合给完整基因组或者完整转录本算 GC 含量。默认值怎么设要看项目习惯。我个人的默认是True因为这样得到的结果在比较不同来源序列时不至于被模糊位点干扰。至于性能直接用四次seq.count()是没问题的Python 的字符串count是 C 层实现速度很快。如果非要一个循环内同时统计 G 和 C也可以写def gc_content_fast(seq): seq seq.upper() gc total 0 for base in seq: if base G or base C: gc 1 total 1 elif base A or base T: total 1 return gc / total * 100 if total else 0.0这段代码比连续count更省内存遍历但对大多数几千到几万碱基的序列差异基本感觉不到。我的体会是可读性优先除非你要处理一整条染色体级别的序列否则不用刻意优化。2.3 滑动窗口 GC 分布一条序列的局部特征整体 GC 含量掩盖了序列内部的不均匀性。比如一个基因组平均 GC 是 40%可能某些区段只有 25%某些区段高达 65%。这种局部特征在引物设计、CpG 岛识别、宏基因组 binning 里都很关键。配合上面的gc_content可以写一个滑动窗口def gc_window(seq, window100, stepNone): if step is None: step window seq seq.upper() result [] for i in range(0, len(seq) - window 1, step): result.append(gc_content(seq[i:iwindow])) return result默认stepwindow表示窗口之间不重叠适合快速扫描整个序列。想让窗口有重叠就传一个比窗口小的step比如窗口 100、步长 50。用的时候gc_values gc_window(ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG, window10, step5)这里要注意循环的边界range(0, len(seq) - window 1, step)保证最后一个窗口的结束位置不超过len(seq)不会切片出长度不足的尾巴。如果你希望统计到最后剩余的一小段可以在循环外面再补一个尾巴判断。不过对于 GC 分布可视化来说丢掉末尾不足一个窗口的序列通常影响不大。3. DNA转录、反向序列与反向互补正反都要会3.1 转录 DNA 到 RNA最容易被“优雅”写法坑到的地方DNA 转录成 RNA核心操作就是 T 替换成 U。Python 字符串的replace是最直接的实现def dna_to_rna(seq): return seq.upper().replace(T, U)这个函数很简单但要注意一个顺序陷阱如果输入序列已经含有Ureplace(T, U)不会碰U结果还是保留原来的U。这在严格区分 DNA / RNA 的场景下通常是合理的因为输入如果混有 U说明本身就不是标准 DNA。还有一个更高级的写法是用str.translate这样能一次性完成替换_DNA_TO_RNA_TRANS str.maketrans(ATGC, AUGC) def dna_to_rna_translate(seq): return seq.upper().translate(_DNA_TO_RNA_TRANS)用translate的好处是查表速度更快同时如果以后想同时处理多个字符映射不用写一堆replace。不过这个函数没有保留非 ATCG 字符比如N会原样保留这是符合预期的。我的建议是日常用的短序列用replace完全够可读性最好如果你在写一个高频调用的库函数用预编译的translate表更稳。3.2 反向序列和反向互补两个方向别混了标题里提到“反向序列获取”这在实际使用中有两个概念单纯的反向reverse_seq把序列倒过来。反向互补reverse_complement先取互补链再反向。做分子生物学实验时大家口中的“取反向”很多时候其实是反向互补。比如设计 PCR 引物时下游引物是模板链的反向互补序列BLAST 比对时常用反向互补来匹配另一条链。所以这两个函数必须同时存在并且命名要清晰。def reverse_seq(seq): return seq.upper()[::-1] def reverse_complement(seq): complement {A: T, T: A, G: C, C: G, N: N, U: A} seq seq.upper() return .join(complement.get(base, N) for base in reversed(seq))需要注意几个处理细节函数里用了reversed(seq)直接在迭代顺序上从末尾走到开头省了先 reverse 再遍历的一步。complement字典里把U的互补也定义为A这样输入 RNA 时也能得到合理的互补链。对于其他非标准字符如R、Y、-字典里没有会统一返回N这在严格处理时其实不够好因为R的互补是Y。如果你要支持 IUPAC 模糊码需要再补一套对应关系。另一个更 Python 的实现是使用str.translate和[::-1]_COMPLEMENTS str.maketrans(ATCGN, TAGCN) def reverse_complement_fast(seq): return seq.upper().translate(_COMPLEMENTS)[::-1]这个版本性能更高但模板里的字符映射写死遇到U会被原样保留不是A。所以如果你的数据是 RNA请用字典版本或者额外补一个maketrans(ATCGUN, TAGCNA)。我在实际项目里踩过一次坑拿一个 RNA 序列做反向互补用的却是 DNA 的maketrans结果U没有互补成A导致下游引物设计全错。从那以后我养成了习惯反向互补函数里必须显式写明“输入是 DNA 还是 RNA”至少要在注释里说清楚默认处理方式。3.3 编码链和模板链谁转录谁DNA 转录 RNA 听起来只是 T 到 U但有一个容易搞混的生物学前提真核细胞转录时用的是模板链RNA 序列与编码链序列一致只是 T 变 U与模板链反向互补。也就是说我们从 FASTA 里拿到的基因组序列通常默认是编码链的序列。因此用dna_to_rna直接把 T 变成 U 是合理的。但如果你手里拿的是某个生物信息软件输出的“模板链序列”直接 T 变 U 就反了。需要先取反向互补再转录。这一步很多新手栽过。所以我通常会在函数注释里写清楚def dna_to_rna(seq, is_templateFalse): seq seq.upper() if is_template: seq reverse_complement(seq) return seq.replace(T, U)这个参数不是生物意义上严格建模只是给实际调用者一个提醒转录前想清楚这条链是编码链还是模板链。4. RNA翻译为蛋白质遗传密码与阅读框4.1 标准遗传密码表手工建立还是直接引入库RNA 翻译成蛋白质需要一张密码子表。64 个密码子对应 20 种氨基酸加 3 个终止密码子。有没有必要自己把整张表写出来我的看法是如果你要练习底层的翻译逻辑就自己写如果你的目标是产出稳定结果直接用 BioPython 的CodonTable更省心。但考虑到标题说的是“相关函数”而且很多人并不想为了一个小工具引入重量级依赖我下面给出一份可直接用的标准表。RNA_CODON_TABLE { UUU: F, UUC: F, UUA: L, UUG: L, UCU: S, UCC: S, UCA: S, UCG: S, UAU: Y, UAC: Y, UAA: *, UAG: *, UGU: C, UGC: C, UGA: *, UGG: W, CUU: L, CUC: L, CUA: L, CUG: L, CCU: P, CCC: P, CCA: P, CCG: P, CAU: H, CAC: H, CAA: Q, CAG: Q, CGU: R, CGC: R, CGA: R, CGG: R, AUU: I, AUC: I, AUA: I, AUG: M, ACU: T, ACC: T, ACA: T, ACG: T, AAU: N, AAC: N, AAA: K, AAG: K, AGU: S, AGC: S, AGA: R, AGG: R, GUU: V, GUC: V, GUA: V, GUG: V, GCU: A, GCC: A, GCA: A, GCG: A, GAU: D, GAC: D, GAA: E, GAG: E, GGU: G, GGC: G, GGA: G, GGG: G, }这张表以 RNA 密码子为键氨基酸单字母缩写为值终止密码子用*表示。注意这里用的是 U 而不是 T所以翻译函数的入参必须是 RNA 序列。如果你手头是 DNA 序列第一步要先把 T 替换成 U。4.2 翻译函数按三个碱基一组切片标准化翻译逻辑其实就是一个固定间隔的循环从某个位置开始每三个碱基取一个密码子查表遇到终止密码子停止把氨基酸追加到列表里。def translate_rna(rna, reading_frame0): rna rna.upper().replace(T, U) amino_acids [] for i in range(reading_frame, len(rna) - 2, 3): codon rna[i:i3] aa RNA_CODON_TABLE.get(codon, X) if aa *: break amino_acids.append(aa) return .join(amino_acids)这里有几个实现细节值得展开说。首先循环边界用len(rna) - 2保证每次切片都是完整的三碱基密码子。比如序列长度是 10那么在i6时切片rna[6:9]是三碱基i9时rna[9:12]越界循环条件已经拦住。最后剩下的一个或两个碱基会直接被忽略。要不要警告末尾不完整的问题我一般不会在函数里加返回值而是在调用前检查len(rna) % 3如果是非零就打印提示。其次RNA_CODON_TABLE.get(codon, X)的意义是防止密码子在表中查不到。正常情况下 64 个密码子都在表里但如果输入序列里有N比如NNU这个切片就不在表中返回X可以明确标识未知氨基酸。我见过有人在这里直接加continue导致蛋白序列长度比预期短反而不容易发现异常。第三遇到终止密码子直接 break不会把*加入蛋白序列。这是最常规的做法。如果你需要保留终止符号可以改写成amino_acids.append(aa)而不是break函数只返回字符串由调用者决定是否要保留*。不同分析场景对终止密码子的处理不一样预测 CDS 时通常去掉做全基因组注释时有时候保留。4.3 起始密码子和开放阅读框从 AUG 开始更合理上面的translate_rna默认从序列的第 0 位开始翻译。但真实 CDS 翻译通常从第一个 AUG 起始密码子开始前面如果没有完整的 AUG翻译出来的蛋白会是一堆乱码。所以我还常用一个变体def translate_orf(rna): rna rna.upper().replace(T, U) start_idx rna.find(AUG) if start_idx -1: return return translate_rna(rna, start_idx)这个函数先找到第一个AUG然后从那里开始按三碱基翻译。它会自动忽略起始密码子前面的序列直到下一个终止密码子结束。这样拿到的序列基本对应一个开放阅读框。不过rna.find(AUG)有个局限如果序列里存在多个阅读框每个阅读框的 AUG 起始位置不同这个函数只会找到最靠前的那个。要完整预测所有可能的开放阅读框需要遍历三个阅读框依次寻找AUG并翻译然后比较长度或按阈值筛选。这种预测更接近真实基因预测但计算量也会上去。如果拿一段 DNA 序列直接调用翻译最好把转录和翻译合并成一个总流程def translate_dna(dna, start_from_atgTrue): rna dna_to_rna(dna) if start_from_atg: return translate_orf(rna) return translate_rna(rna)用一段示例序列跑一遍seq ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG rna dna_to_rna(seq) protein translate_rna(rna) print(seq) # ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG print(rna) # AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCG AUAG?我实际跑的结果是RNA AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCG AUAG翻译出来是MAIVMGR*。这里需要注意序列里有两个 AUG 起始密码子而translate_rna从第 0 位开始翻译遇到第一个终止密码子就停了。如果你想要从第二个 AUG 开始的 ORF就得用translate_orf或者指定reading_frame。这样的坑在真实基因序列中很常见。所以我的建议是如果目标是预测一个基因的蛋白优先写translate_orf而不是translate_rna默认从头翻译。5. 常见问题与排查技巧实录写这套函数的过程中我遇到过不少看起来莫名其妙的报错和结果异常。整理成下面的速查表能帮大家少走弯路。典型现象可能原因处理方法翻译返回空字符串RNA 序列里找不到 AUG检查序列方向尝试反向互补后再找GC 含量比预期低很多分母包含了 N 或 gap调用gc_content时设置ignore_ambiguousTrue反向互补结果和预期完全不一样先互补了但没有反向确认使用的是reverse_complement而不是complement密码子表 KeyError传入的是 DNA 而不是 RNA翻译前先replace(T,U)序列长度统计多了几个字符读取时没去掉换行符解析 FASTA 时对每行strip()重复的 ID 造成序列被覆盖FASTA 存在相同 ID改用列表结构或在读入时检测重复蛋白序列末尾比预期短很多提前遇到终止密码子确认是否选对了阅读框或序列方向5.1 最容易被忽略的大小写问题几乎所有函数我都统一转成大写但有一个隐患最早版本的dna_to_rna如果没做.upper()传入小写序列时replace(T,U)不会替换t会留下小写 t。后续统计和翻译时这个小写 t 既不在预期位置也查不到密码子最终结果出现奇怪的X或直接长度不对。现在我养成了习惯在函数的入口统一seq.upper()不要在调用方控制。如果哪天上游数据混入了小写碱基所有函数都会自动规避这个问题。5.2 内存不够大 FASTA 的流式处理如果你的输入是几百 MB 甚至几个 GB 的基因组 FASTA前面那个read_fasta返回完整字典会把内存吃满。此时应该用生成器一次只返回一条记录。def read_fasta_stream(filepath): seq_id None seq_parts [] with open(filepath, r, encodingutf-8) as f: for line in f: line line.strip() if line.startswith(): if seq_id is not None: yield seq_id, .join(seq_parts) seq_id line[1:].split()[0] seq_parts [] elif line: seq_parts.append(line.upper()) if seq_id is not None: yield seq_id, .join(seq_parts)调用方式for seq_id, seq in read_fasta_stream(large_genome.fasta): gc gc_content(seq) print(seq_id, gc)这里用yield而不是return所以一次只在内存中保留一条序列。如果你只对每条染色体做 GC 统计这种模式是完美的。要注意的是单条染色体序列也可能非常大拼接的时候仍然会占用这条染色体长度的内存但这是不可避免的因为我们需要完整序列做全局计算。5.3 方向搞反一个信息学模型之外的生物学问题很多次收到同学反馈“我的反向互补怎么跟参考基因组对不上”。我第一反应是问你手里的序列是正链还是负链FASTA 里的序列通常代表正链但一个基因如果位于基因组负链它的 CDS 序列在 FASTA 里未必以编码方向存储。此时你直接翻译得到的蛋白是不对的需要先对原始序列做反向互补再找起始密码子。这已经不是 Python 语法层面的问题而是生物学方向问题。我的建议是在写任何一条工具链前先明确输入数据的链方向并在函数命名附近写好注释。否则这类错误最难排查因为代码没有任何报错只是结果不对。另外还有一个小细节如果序列是 RNA 且带U进行反向互补我用的是U - A的映射。标准配对上RNA 互补链对应 DNA 模板链时A互补为U还是T需要根据下游需求决定。大部分引物设计工具会直接按 DNA 处理所以如果你传入 RNA 给我这个函数要仔细核对输出是否符合预期。收尾两个实用小技巧最后分享两个我在实际使用中沉淀下来的小技巧。第一个如果你想把所有函数打包成一个小模块建议别再单独编写count_nucleotides、count_bases、gc_content等一堆零散函数而是统一设计一个SequenceRecord类把序列和所有方法装在一起。这样调用起来更顺手class SequenceRecord: def __init__(self, seq_id, seq): self.id seq_id self.seq seq.upper() self.rna None self.protein None def gc(self, ignore_ambiguousTrue): return gc_content(self.seq, ignore_ambiguous) def transcribe(self): self.rna dna_to_rna(self.seq) return self.rna def translate(self): if self.rna is None: self.transcribe() self.protein translate_rna(self.rna) return self.protein类的好处是减少传递参数的次数尤其是当序列和它的转录产物、翻译产物需要同时存在时状态管理更清晰。如果只是写一次性脚本类可能显得有点重但长期维护的项目里非常值得。第二个技巧是自己写assert校验。在函数里加入几条验证逻辑能提前拦截大部分错误输入assert set(seq.upper()).issubset({A, T, G, C, N, U}), 序列中包含非核酸字符这条断言能帮你在统计 GC 之前就发现序列里混入了空格、数字或者其他字符。不过要注意真实测序数据偶尔会有小写字符和 IUPAC 模糊码直接断言会误伤。你可以把断言放在调试模式下或者只对已知的干净数据启用。这些函数单独看都不复杂但当它们组合成一套工具就能支撑起一个简单的基本分析流程。我是建议每个人在入门生信时都亲手实现一遍再换 BioPython这样你对后面大框架里“替你偷懒”的部分会有更清楚的认识。
