在MacOS、Unix系统中用awk从Stockholm格式文件提取两个独立字符串

假设你需要提取的是序列ID（比如WP_002855993.1/5-168）和对应的纯氨基酸序列（去除所有非字母字符），或者是提取#=GS行中的功能描述（DE字段）和序列ID，下面分两种常见场景给出解决方案：

场景1：提取序列ID与纯氨基酸序列 #

这个awk命令会自动识别非注释行（行首不是#的行），提取第一个字段作为序列ID，再把该行里的短横线、点、空格等非字母字符全部清除，得到干净的氨基酸序列：

/^[^#]/ {
    id = $1
    seq = substr($0, length(id)+1)
    gsub(/[^A-Za-z]/, "", seq)
    print "ID: " id
    print "Sequence: " seq
    print "---"
}

运行方式 #

把上述代码保存为extract_stockholm_seq.awk，然后在终端执行：

awk -f extract_stockholm_seq.awk your_stockholm_file.sto

针对你提供的示例内容，输出结果会是：

ID: WP_002855993.1/5-168
Sequence: LEHNGKKYSDKDLIDAFYQLGIKRGDILCVHTELmkfgKALLTKNNDFLKTLLECFFKVLGKEGTLLMPTFTYSFCKNEVYDKVHSKGKVGVLNEFFRTSGgGVRRTSDPIFSFAVKGAKADIFLKENSSCFGKDSVYEILTREGGKFMLLG
---

场景2：提取#=GS行的功能描述与对应序列ID #

如果需要提取注释行里的功能说明和关联的序列ID，用这个awk命令就能实现：

/^#=GS/ {
    id = $2
    # 从第4个字段开始拼接完整的功能描述
    desc = ""
    for(i=4; i<=NF; i++) {
        desc = desc " " $i
    }
    # 去掉描述开头多余的空格
    sub(/^ /, "", desc)
    print "ID: " id
    print "Description: " desc
    print "---"
}

针对你提供的示例内容，输出结果会是：

ID: WP_002855993.1/5-168
Description: [subseq from] MULTISPECIES: AAC(3) family N-acetyltransferase [Campylobacter]
---
ID: WP_002856586.1/5-166
Description: [subseq from] MULTISPECIES: aminoglycoside N(3)-acetyltransferase [Campylobacter]
---

如果你的需求是提取其他类型的字符串，可以补充说明具体要提取的内容，我再调整命令~

内容的提问来源于stack exchange，提问作者D.Parker