如何批量计算核酸序列列表的GC含量并定位非核苷酸字符

需求说明 #

需要批量处理一组核酸序列，实现两个核心功能：

1. 计算每条序列的GC含量（G/C碱基占总长度的百分比，不区分大小写）
2. 精准定位包含非核苷酸字符的序列（如seq2中的'X'）

改进后的代码 #

# 定义所有核酸序列
sequences = [
    'CCACGCGTCCGCCGCGACCTGCGTTTTCCTGGGGGTCCGCAACTCTGGCTTGACCCAAGGACCCGGCCAC',
    'attgccattatataACCCGGCCACCCCCATAGGCAGATGTCAGGACAACTCGCATCTCAGCAGAGCAGCCCCTGGCCCAGG',
    'TCXCACCCATAGGCAGATGGCCTCCGCCCCACCCCCGGGAGGATTTCTTAATGGGGTGAAAATGC',
    'CAGTCCCCGAAGCCAGGGTTCCGGGACCCCCGGGGCCGAGCTGGGCGCGGGAAAAGAAttacggacttaGTCAGCCCCGCAGGGG',
    'ATGGGGTGATCGTCGCTCGCGGGCTCTGTCTTCCTGTTCACCCTCCTCTGCCCCCAACTCCATCTCTGAGACCTCCTGCCCCCCCA',
    'AAAAAAGAAGTCGCTCGCGTCGCTCGCGGGCTGGGCTCTGTCTGCGTCGCTCGCGGGCTAGAGAGCCAGGGTGA'
]

# 合法核苷酸集合（大小写兼容）
valid_nucleotides = {'G', 'A', 'C', 'T', 'U', 'g', 'a', 'c', 't', 'u'}

# 批量处理每条序列
for idx, seq in enumerate(sequences):
    seq_name = f'seq{idx}'
    # 检查是否存在非法字符
    invalid_chars = [c for c in seq if c not in valid_nucleotides]
    if invalid_chars:
        print(f'ERROR: non-nucleotide characters present in {seq_name}')
    
    # 计算GC含量（基于原始序列长度，若需用过滤后长度可修改len(seq)为len(filtered_seq)）
    g_count = seq.count('G') + seq.count('g')
    c_count = seq.count('C') + seq.count('c')
    gc_content = (g_count + c_count) * 100 / len(seq) if len(seq) > 0 else 0
    print(f'The GC content of {seq_name} is {gc_content:.2f} %')

代码说明 #

• 序列存储优化：直接将所有序列放入列表，方便循环批量处理
• 非法字符检测：遍历每条序列的每个字符，收集非法字符，若存在则输出对应序列名称，精准定位异常
• 批量GC计算：通过enumerate获取序列索引，统一计算每条序列的G/C碱基数量，再计算百分比（保留两位小数提升可读性）
• 容错处理：判断序列长度不为0，避免出现除以0的错误

运行输出示例 #

The GC content of seq0 is 70.00 %
The GC content of seq1 is 56.41 %
ERROR: non-nucleotide characters present in seq2
The GC content of seq2 is 53.85 %
The GC content of seq3 is 61.29 %
The GC content of seq4 is 55.17 %
The GC content of seq5 is 48.39 %

内容的提问来源于stack exchange，提问作者Khan Inan