正则表达式：FASTA解析、日志提取、awk/sed

2700 字
14 分钟
正则表达式：FASTA解析、日志提取、awk/sed
2025-04-15
技术
文本处理
/
生信
/
Linux

生信数据大多是纯文本——FASTA、FASTQ、SAM、GFF、BED、VCF。正则表达式是处理这些文本格式的最短路径。一行 sed 或 grep 即可替代十几步骤的手工操作。本文覆盖生信中的正则实战：FASTA 序列头解析、FASTQ 质量提取、GFF 注释挖掘、日志错误提取、awk/sed 替换。

实测环境：Debian 12，GNU grep 3.11，GNU sed 4.9，GNU awk 5.3。

1. 正则基础——只用这 10 个元字符#

不用学全。生信中 90% 的场景只需要这 10 个：

元字符含义生信示例
.任意单个字符chr. 匹配 chr1, chr2, chrX…
*前一个字符重复0次或多次A* 匹配空、A、AA、AAA…
+前一个字符重复1次或多次[AT]+ 匹配序列
?前一个字符可选chr?1 匹配 chr1 或 ch1
^行首^> 匹配 FASTA 头
$行尾fastq.gz$ 文件名尾
[]字符集[ATCG] 匹配任一碱基
[^]排除字符集[^>] 匹配非>字符
()捕获组>(\w+) 提取序列名
|或fastq|fq 两种后缀
\s空白字符匹配空格、tab
\w单词字符[a-zA-Z0-9_]

1.1 BRE vs ERE vs PCRE——别混用#

Linux 文本工具对正则的支持不一样：

工具默认模式用 + 需要用 () 需要
grepBRE（基本正则）\+\(\)
grep -EERE（扩展正则）+()
grep -PPCRE（Perl兼容）+()
sedBRE\+\(\)
sed -EERE+()
awkERE+()

踩坑最多的是 grep 和 sed 默认是 BRE。 我现在的习惯是：grep 一律加 -P（PCRE），sed 一律加 -E。因为生信正则里 +、?、| 太常用，不想每次转义。

2. FASTA/FASTQ 格式——序列文件的正则操作#

2.1 解析 FASTA 序列头#

NCBI RefSeq 的 FASTA 头和 ENSEMBL 的完全不同：

# NCBI RefSeq
>NM_001126114.2 Homo sapiens tumor protein p53 (TP53), mRNA
>NC_000001.11 Homo sapiens chromosome 1, GRCh38.p14


# ENSEMBL
>ENST00000269305.9 cdna chromosome:GRCh38:17:7661779:7687538:1
>ENSG00000141510 gene_biotype:protein_coding


# UCSC
>chr1_1000_2000

用正则提取各种 FASTA 头的关键字段：

Terminal window
# ===== 从不同类型FASTA头中提取信息 =====


# NCBI：提取 accession、物种、基因名
grep '^>' refseq.fasta | sed -E 's/^>([^ ]+) ([^ ]+) (.+)/ACCN=\1 SPECIES=\2 DESC=\3/'


# 输出：ACCN=NM_001126114.2 SPECIES=Homo DESC=sapiens tumor protein p53 (TP53), mRNA


# ENSEMBL：提取转录本ID、染色体、坐标
grep '^>' ensembl.fasta | grep -oP '>(ENST\d+)'                    # 转录本ID
grep '^>' ensembl.fasta | grep -oP 'chromosome:[^:]+:([0-9]+)'     # 坐标


# UCSC：提取染色体和区间
grep '^>' ucsc.fasta | grep -oP '>(chr\w+)_(\d+)_(\d+)'
# 用 sed 分别提取
sed -n 's/^>\(chr\w\+\)_\([0-9]\+\)_\([0-9]\+\)/CHROM=\1 START=\2 END=\3/p' ucsc.fasta

2.2 序列碱基统计——不用 seqkit 也能快速检查#

Terminal window
# 统计 ATGC 组成（去掉序列头行）
grep -v '^>' sequence.fasta | fold -w1 | sort | uniq -c


# 输出示例：
#  15234 A
#   9876 C
#  10021 G
#  14892 T
#      3 N       ← 有3个不确定碱基，注意！


# 计算 GC 含量
total=$(grep -v '^>' sequence.fasta | tr -d '\n' | wc -c)
gc=$(grep -v '^>' sequence.fasta | tr -d '\n' | grep -oP '[GCgc]' | wc -l)
echo "GC content: $(echo "scale=2; $gc * 100 / $total" | bc)%"

GCGC\\% = \\frac{G + C}{A + T + G + C} \\times 100\\%

2.3 FASTQ 质量线提取#

FASTQ 每 4 行一组：@header / sequence / + / quality。提取质量信息：

Terminal window
# 提取所有质量行（每4行的第4行）
awk 'NR % 4 == 0' sample.fastq | head -5


# 统计 Q>30 的碱基比例（质量字符 ASCII-33=Q值）
awk 'NR % 4 == 0 {
    len = length($0)
    for (i = 1; i <= len; i++) {
        q = ord(substr($0, i, 1)) - 33
        if (q >= 30) q30++
        total++
    }
}
END {
    printf("Q30: %.1f%% (%d/%d)\n", q30*100/total, q30, total)
}
function ord(c) { return index("!\"#$%&\047()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~", c) + 32 }
' sample.fastq

当然，生产环境用 seqkit stats 或 fastp 更快，但理解这个正则逻辑能帮你在没有专用工具时自救。

3. GFF/GTF 注释文件的挖掘#

GFF（General Feature Format）是基因注释的标准格式。9 列，tab 分隔：

chr1  HAVANA  gene  11869  14409  .  +  .  gene_id=ENSG00000223972;gene_name=DDX11L1
chr1  HAVANA  exon  11869  12227  .  +  .  gene_id=ENSG00000223972;transcript_id=ENST00000456328

3.1 提取特定类型——如所有 exon#

Terminal window
grep '\texon\t' annotation.gtf | head -5

\t 是你最好的朋友——GFF/GTF 全是 tab 分隔，用 \t 精确匹配列，避免 exon 出现在第 9 列注释里被误匹配。

3.2 从第9列提取关键字段#

GFF 第 9 列用 ; 分隔键值对。提取所有基因的 gene_id 和 gene_name：

Terminal window
# 从 GTF 中提取基因的 ID 和名称
grep '\tgene\t' annotation.gtf | \
    grep -oP 'gene_id "([^"]+)"|gene_name "([^"]+)"' | \
    paste - - | \
    sed 's/gene_id "//g;s/"//g;s/gene_name "//g;s/"\t/ /g'

分解这个管道：

  1. grep '\tgene\t' → 只取 gene 行
  2. grep -oP 'gene_id "([^"]+)"|gene_name "([^"]+)"' → 只输出匹配部分（-o = only matching）
  3. paste - - → 把两行（gene_id、gene_name）合并到一行
  4. sed → 清理引号

输出：

ENSG00000223972 DDX11L1
ENSG00000227232 WASH7P
ENSG00000278267 MIR6859-1

3.3 awk 解析 GFF 的复杂统计#

Terminal window
# 统计每个染色体上各类 feature 的数量
awk -F'\t' '$3 == "gene" {
    split($9, attrs, ";")
    for (a in attrs) {
        if (attrs[a] ~ /gene_biotype/) {
            split(attrs[a], bt, " ")
            gsub(/"/, "", bt[2])
            count[$1][bt[2]]++
        }
    }
}
END {
    for (chr in count) {
        for (bt in count[chr]) {
            print chr, bt, count[chr][bt]
        }
    }
}' annotation.gtf | sort -k1,1 -k3,3rn | head -20

这里用到的正则 ~ /gene_biotype/ 是 awk 的语法糖，等价于 match。

4. 日志解析——从生信流程日志中提取关键信息#

4.1 提取比对率#

STAR、BWA、HISAT2 的输出日志都有比对率信息，但格式各不同：

Terminal window
# STAR 日志
grep -oP 'Uniquely mapped reads % \| \K[0-9.]+' star_Log.final.out


# BWA 日志（从 stderr）
bwa mem ... 2>&1 | grep -oP '[0-9.]+%'


# HISAT2 日志
grep -oP 'overall alignment rate: \K[0-9.]+%' hisat2.log

\K 是 PCRE 的特性：丢弃 \K 左边匹配的内容，只保留右边。非常方便。

4.2 批量提取所有样本的比对率——一行搞定#

Terminal window
# 假设每个样本的 STAR 日志在 logs/SAMPLE_Log.final.out
for log in logs/*_Log.final.out; do
    sample=$(basename "$log" _Log.final.out)
    rate=$(grep -oP 'Uniquely mapped reads % \| \K[0-9.]+' "$log")
    echo "$sample $rate"
done

4.3 从 fastp JSON 报告提取过滤率（用 jq + 正则）#

Terminal window
# fastp JSON 里过滤统计
jq -r '.summary.before_filtering.total_reads, .summary.after_filtering.total_reads' fastp_report.json


# 计算过滤率
before=$(jq '.summary.before_filtering.total_reads' fastp_report.json)
after=$(jq '.summary.after_filtering.total_reads' fastp_report.json)
python3 -c "print(f'Filter rate: {(1 - $after/$before)*100:.1f}%')"

4.4 从命令行历史中提取你常用的长命令#

Terminal window
# 找出过去一周用得最多的10个生信命令（要求：命令长度>30字符）
history | sed -E 's/^[ ]*[0-9]+[ ]+//' | \
    grep -P '^.{30,}' | \
    sort | uniq -c | sort -rn | head -10

5. 实用 Python 正则#

Shell 正则够用但长文本处理有时力不从心。Python 的 re 模块配合命名捕获组（named groups）更清晰：

#!/usr/bin/env python3
"""regex bioinfo examples - Python re module"""


import re


# ===== 1. 解析 FASTA 头（各种格式） =====
fasta_headers = [
    ">NM_001126114.2 Homo sapiens tumor protein p53 (TP53), mRNA",
    ">ENST00000269305.9 cdna chromosome:GRCh38:17:7661779:7687538:1",
    ">chr1_1000_2000",
    ">ENSG00000141510|TP53|protein_coding",
]


# 模式1：NCBI RefSeq
ncbi_pattern = re.compile(
    r'>(?P<accession>[A-Z]+_\d+(?:\.\d+)?) '
    r'(?P<species>\w+ \w+) '
    r'(?P<description>.+)'
)


# 模式2：ENSEMBL 转录本
ensembl_transcript = re.compile(
    r'>(?P<transcript_id>ENST\d+(?:\.\d+)?) '
    r'.*chromosome:(?P<build>\w+):(?P<chr>\w+):(?P<start>\d+):(?P<end>\d+)'
)


# 模式3：UCSC 坐标
ucsc_pattern = re.compile(
    r'>(?P<chrom>chr\w+)_(?P<start>\d+)_(?P<end>\d+)'
)


for header in fasta_headers:
    for pattern in [ncbi_pattern, ensembl_transcript, ucsc_pattern]:
        m = pattern.match(header)
        if m:
            print(f"Matched: {m.groupdict()}")
            break
    else:
        print(f"No match: {header}")

输出：

Matched: {'accession': 'NM_001126114', 'species': 'Homo sapiens',
          'description': 'tumor protein p53 (TP53), mRNA'}
Matched: {'transcript_id': 'ENST00000269305', 'build': 'GRCh38',
          'chr': '17', 'start': '7661779', 'end': '7687538'}
Matched: {'chrom': 'chr1', 'start': '1000', 'end': '2000'}
No match: >ENSG00000141510|TP53|protein_coding

5.1 从 VCF 提取变异信息#

# ===== 2. 解析 VCF 行 =====
vcf_line = "chr1\t123456\trs123\tA\tG\t100\tPASS\tDP=50;AF=0.35\tGT:AD:DP\t0/1:20,10:30"


# VCF 格式（8列固定 + 第9列INFO + 样本列）
vcf_pattern = re.compile(r'''
    ^(?P<chrom>\S+)\t
    (?P<pos>\d+)\t
    (?P<id>\S+)\t
    (?P<ref>[ATCG]+)\t
    (?P<alt>[ATCG,]+)\t
    (?P<qual>\S+)\t
    (?P<filter>\S+)\t
    (?P<info>.+?)(?:\t|$)
''', re.VERBOSE)


m = vcf_pattern.match(vcf_line)
if m:
    info = dict(kv.split('=') for kv in m.group('info').split(';') if '=' in kv)
    print(f"Variant: {m.group('chrom')}:{m.group('pos')} "
          f"{m.group('ref')}>{m.group('alt')}")
    print(f"  DP={info.get('DP', 'N/A')}, AF={info.get('AF', 'N/A')}")

5.2 批量 GTF 第9列解析#

# ===== 3. 解析 GTF 第9列属性 =====
gtf_attr = 'gene_id "ENSG00000223972"; gene_version "5"; '\
           'gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";'


# 把键值对提取为字典
attr_dict = dict(re.findall(r'(\w+) "([^"]+)"', gtf_attr))
print(attr_dict)
# {'gene_id': 'ENSG00000223972', 'gene_version': '5',
#  'gene_name': 'DDX11L1', 'gene_source': 'havana',
#  'gene_biotype': 'transcribed_unprocessed_pseudogene'}


# 获取特定属性
gene_id = attr_dict.get("gene_id", "N/A")
gene_name = attr_dict.get("gene_name", "N/A")
biotype = attr_dict.get("gene_biotype", "N/A")

正则 (\w+) "([^"]+)" 的解读：

  • (\w+) → 捕获属性名（如 gene_id）
  • "([^"]+)" → 捕获双引号内的值

6. sed 实战——生信中的数据清洗#

6.1 批量修改文件名（在文本列表里）#

Terminal window
# 把样本列表从 "/path/to/S01_R1.fastq.gz" 改成 "S01"
sed -E 's|.*/(S[0-9]+)_R1\.fastq\.gz|\1|' file_list.txt

6.2 去 FASTA 文件的换行（不规范的FASTA）#

Terminal window
# 有的FASTA每行碱基数不同，需要标准化为每行60个
# 1. 先变成单行序列
awk '/^>/ {if (seq) print seq; print; seq=""; next}
     {seq=seq $0}
     END {if (seq) print seq}' input.fasta > oneline.fasta


# 2. 再把非头行按60字符折行
fold -w 60 oneline.fasta > formatted.fasta

6.3 批量清理日志中的 ANSI 颜色码#

Terminal window
# 从程序输出中去掉颜色码（\033[...m）
sed -E 's/\x1b\[[0-9;]*m//g' program.log > clean.log

7. 踩坑记录#

坑1：贪婪匹配 vs 非贪婪匹配#

这是正则里最常见的坑。生信中特别容易在解析 FASTA 头时碰到：

Terminal window
# 贪婪匹配（默认）——会吃太多
echo ">geneA description; geneB description" | grep -oP '>.*;'
# 输出：>geneA description; geneB description;  ← 吃了整行！


# 非贪婪匹配——加 ?
echo ">geneA description; geneB description" | grep -oP '>.*?;'
# 输出：>geneA description;  ← 正确停在第一个分号

quantifier_greedy⇒maximal match，quantifier_lazy⇒minimal matchquantifier\_greedy \Rightarrow maximal\ match， quantifier\_lazy \Rightarrow minimal\ match

.* = 尽可能多匹配，.*? = 尽可能少匹配。解析分隔符分隔的文本时，一律用非贪婪 .*?（或 [^分隔符]+）。

坑2：grep 不带 -P 时 \d 不生效#

Terminal window
# 用 grep（默认 BRE）找数字
grep '\d+' file.txt       # 不匹配！BRE 不认识 \d


# 正确做法
grep -P '\d+' file.txt    # PCRE 模式
grep -E '[0-9]+' file.txt # ERE 模式用传统字符集

注意： \d、\w、\s 这些简写只在 PCRE 里有效。grep 和 sed 默认不认。

坑3：sed 中 -E 和 \1 的配合#

Terminal window
# 错误：-E 模式下引用捕获组仍用 \1
sed -E 's/>(ENSG[0-9]+).*/Gene: \1/'   # 正确


# 不要写成 $1（那是 Perl 的语法）
sed -E 's/>(ENSG[0-9]+).*/Gene: $1/'   # 错误！sed 用 \1

坑4：跨行匹配——grep 天生只在一行内匹配#

Terminal window
# 你用 grep 找跨两行的模式——永远找不到
grep 'START.*END' file.txt  # 只匹配 START 和 END 在同一行


# 跨行匹配用 pcregrep 或 awk
pcregrep -M 'START\n.*END' file.txt


# 或者用 awk
awk '/START/,/END/' file.txt

这个坑在处理带换行的 FASTA 序列时最常见。比如你想找包含特定 motif 的序列，但这个 motif 恰好在换行边界。建议先把 FASTA 转成单行再匹配。

坑5：过度使用 .* 导致灾难性回溯（catastrophic backtracking）#

# 这个正则看起来无害但能在特定输入上卡死几分钟
import re
pattern = re.compile(r'(A+)+B')
text = "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC"  # 没有 B
pattern.match(text)  # 可能卡死！回溯次数指数级增长

生信里如何避免： 当重复模式（如 [ATCG]+）嵌套使用时，用 (?> ... ) 原子组（atomic group）或改用 finditer() 而不是 match()。

# 安全写法
pattern = re.compile(r'A+B')  # 别嵌套重复

8. 生信正则速查表#

场景正则工具
匹配 FASTA 头^>grep '^>'
非 FASTA 头行^[^>]grep -v '^>'
提取 NCBI accession>(\w+_\d+(\.\d+)?)grep -oP
提取 ENSEMBL IDENS[GT]\d{11}grep -oP
GTF tab 分隔列\tawk -F'\t'
GTF 第9列属性(\w+) "([^"]+)"grep -oP
匹配质量字符（Q≥30）[?@A-I]（以 ! = 0 为准）grep -oP
染色体名chr(\w+)grep -oP
基因组坐标(\d+)[-:](\d+)sed -E
日志百分比([0-9]+(?:\.[0-9]+)?)%grep -oP
非贪婪匹配.*?PCRE
提取 tab 文件第n列cut -f n（比正则简单）cut

9. 总结#

正则不是学一次就会的。我用了五年，每次碰到新格式还是会查手册。但有几个路径能帮你少走弯路：

  1. 先用简单命令验证——grep -c 看匹配到了多少行，确认没多匹配也没漏
  2. 逐步构建正则——先写最简单的 ^>，确认能匹配 FASTA 头；再加 ^>\w+，再加捕获
  3. grep -oP 是最好的调试器——-o 只输出匹配部分，让你直观看到正则到底吃掉了什么
  4. 生信专用工具优先——能用 seqkit 就别手写 awk 正则，能用 bedtools 就别手动解析 BED

本文于 2025-10-10 在 Debian 12 上实测完成。GNU grep 3.11, GNU sed 4.9, GNU awk 5.3, Python 3.10。所有正则均在生信数据上实测验证。

文章分享

如果这篇文章对你有帮助，欢迎分享给更多人！

正则表达式：FASTA解析、日志提取、awk/sed
https://fg.ink/posts/regex-bioinformatics/
作者
风观
发布于
2025-04-15
许可协议
CC BY-NC-SA 4.0
FastQC质检与fastp清洗：测序数据预处理
jq与csvkit：JSON/CSV数据处理
相关文章 智能推荐
1
Vim配置与技巧：语法高亮、折叠、宏、多文件
技术 Vim在生信脚本编辑中的实用配置与技巧，覆盖语法高亮、代码折叠、宏录制和多文件操作。
2
文件编码与换行符：UTF-8/CRLF/BOM诊断与转换
技术 UTF-8、CRLF、BOM等文件编码与换行符问题的诊断与转换方法，覆盖iconv、dos2unix在生信场景下的实战。
3
生信常见报错排查：20条错误与解决方案
技术 生信分析中20个常见报错的排查手册，按文件权限、内存溢出、软件依赖和流程逻辑分类，每条附错误原文与解决方案。
4
Markdown与Pandoc：文档写作与格式转换
技术 Markdown写作与Pandoc格式转换的完整工作流，覆盖实验记录模板、文献引用和学术级PDF输出。
5
命令行效率工具：bat/fd/ripgrep/jq
技术 bat、fd、ripgrep和jq四个命令行效率工具的生信实战，覆盖文件浏览、内容搜索和JSON数据解析。
随机文章 随机推荐
Profile Image of the Author
风观
风有来路，观有所思
分类
标签
站点统计
文章
50
分类
1
标签
29
总字数
61,837
运行时长
0 天
最后活动
0 天前
风观 © 2018-2026   风有来路，观有所思
© 2026 风观. All Rights Reserved. / RSS / Sitemap
Powered by Astro & Firefly

文章目录