如何通过匹配简并模式筛选FASTQ文件？

需求说明 #

从FASTQ文件中筛选出包含特定简并模式的序列及其质量值，模式为..........CAA.....GTGG..........（.代表任意A/C/G/T碱基），即序列需包含「连续10个任意碱基 + CAA + 连续5个任意碱基 + GTGG」的片段。

方法1：Python脚本实现 #

利用Python正则表达式匹配FASTQ四行一组的格式，筛选符合条件的条目，并可按需截取序列和质量值：

import re

# 定义匹配模式：匹配任意10个ATCG碱基后接CAA，再接5个ATCG碱基后接GTGG
match_pattern = re.compile(r'[ATCG]{10}CAA[ATCG]{5}GTGG')

with open('input.fastq', 'r') as in_fq, open('output.fastq', 'w') as out_fq:
    while True:
        # 按FASTQ格式读取四行数据
        header = in_fq.readline()
        if not header:
            break  # 读取到文件末尾
        seq_line = in_fq.readline().strip()
        plus_line = in_fq.readline()
        qual_line = in_fq.readline().strip()
        
        # 检查序列是否符合模式
        if match_pattern.search(seq_line):
            # 按示例输出截取（去掉前3个碱基和末尾2个，可根据需求修改切片）
            output_seq = seq_line[3:-2] if len(seq_line) > 5 else seq_line
            output_qual = qual_line[3:-2] if len(qual_line) > 5 else qual_line
            
            # 写入输出文件
            out_fq.write(header)
            out_fq.write(output_seq + '\n')
            out_fq.write(plus_line)
            out_fq.write(output_qual + '\n')

使用方式：将脚本保存为filter_fastq.py，确保输入文件名为input.fastq，运行python filter_fastq.py即可生成output.fastq。

方法2：Awk命令实现 #

Awk适合快速处理文本文件，可直接通过命令行或脚本实现筛选：

命令行直接运行

awk '
BEGIN { pattern = /[ATCG]{10}CAA[ATCG]{5}GTGG/ }
NR%4 == 1 { h = $0 }
NR%4 == 2 { s = $0 }
NR%4 == 3 { p = $0 }
NR%4 == 0 {
    if (s ~ pattern) {
        # 按示例截取输出，如需完整输出直接打印s和$0即可
        print h
        print substr(s, 4, length(s)-5)
        print p
        print substr($0, 4, length($0)-5)
    }
}
' input.fastq > output.fastq

脚本文件运行

将以下内容保存为filter_fastq.awk：

BEGIN {
    # 定义匹配模式
    pattern = /[ATCG]{10}CAA[ATCG]{5}GTGG/
}

# 四行一组存储数据
NR%4 == 1 { header = $0 }
NR%4 == 2 { seq = $0 }
NR%4 == 3 { plus = $0 }
NR%4 == 0 {
    qual = $0
    # 匹配模式则输出
    if (seq ~ pattern) {
        print header
        print substr(seq, 4, length(seq)-5)  # 截取前3后2之外的部分
        print plus
        print substr(qual, 4, length(qual)-5)
    }
}

运行命令：awk -f filter_fastq.awk input.fastq > output.fastq

正则说明 #

• [ATCG]{10}：匹配任意10个A/C/G/T碱基
• CAA：匹配固定序列
• [ATCG]{5}：匹配任意5个A/C/G/T碱基
• GTGG：匹配固定序列
• 若允许模式出现在序列任意位置（而非必须从开头），直接使用上述正则即可；若要求模式从序列起始位置开始，需在正则前加^。

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