使用 awk 提取两个单独的字符串

Question

MacOS、Unix

所以我有一个以下斯德哥尔摩格式的文件：

# STOCKHOLM 1.0

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

WP_002855993.1/5-168         ------LEHNGKKYSDKDLIDAFYQLGIKRGDILCVHTELmkfgKALLT.K...NDFLKTLLECFFKVLGKEGTLLMP-TF---TYSF------CKNE------VYDKVHSKG--KVGVLNEFFRTSGgGVRRTSDPIFSFAVKGAKADIFLKEN--SSCFGKDSVYEILTREGGKFMLLGLNYG-HALTHYAEE-----
#=GR WP_002855993.1/5-168 PP ......6788899999***********************9333344455.6...8999********************.33...3544......4555......799999975..68********98626999****************999865..689*********************9875.456799996.....
WP_002856586.1/5-166         ------LEFENKKYSTYDFIETFYKLGLQKGDTLCVHTEL....FNFGFpLlsrNEFLQTILDCFFEVIGKEGTLIMP-TF---TYSF------CKNE------VYDKINSKT--KMGALNEYFRKQT.GVKRTNDPIFSFAIKGAKEELFLKDT--TSCFGENCVYEVLTKENGKYMTFGGQG--HTLTHYAEE-----
#=GR WP_002856586.1/5-166 PP ......5566677788889999******************....**9953422246679*******************.33...3544......4455......799998876..589**********.******************99999886..689******************999765..5666***96.....
#=GC PP_cons                 ......6677788899999999*****************9....77675.5...68889*******************.33...3544......4455......799999976..689*******998.8999**************99999876..689******************9998765.466699996.....
#=GC RF                      xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx....xxxxx.x...xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

WP_002855993.1/5-168         -----------------------------------------------------------------------------------------------------
#=GR WP_002855993.1/5-168 PP .....................................................................................................
WP_002856586.1/5-166         -----------------------------------------------------------------------------------------------------
#=GR WP_002856586.1/5-166 PP .....................................................................................................
#=GC PP_cons                 .....................................................................................................
#=GC RF                      xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
//

我创建了一个脚本来提取我想要的 ID，在本例中为 WP_002855993.1 和 WP_002856586.1，然后搜索另一个文件以提取 DNA 序列适当的 ID。脚本如下：

#!/bin/bash

for fileName in *.sto;
do
protID=$(grep -o "WP_.\{0,11\}" $fileName | sort | uniq)
echo $protID
file=$(echo $fileName | cut -d '_' -f 1,2,3)
file=$(echo $file'_protein.faa')
echo $file 
if [ -n "$protID" ]; then
gawk "/^>/{N=0}/^.*$protID/{N=1} {if(N)print}" $file >> 
sequence_protein.file
fi
done

这是我正在查看的文件类型的示例：

>WP_002855993.1 MULTISPECIES: AAC(3) family N-acetyltransferase [Campylobacter]
MKYFLEHNGKKYSDKDLIDAFYQLGIKRGDILCVHTELMKFGKALLTKNDFLKTLLECFFKVLGKEGTLLMPTFT
>WP_002856586.1 MULTISPECIES: aminoglycoside N(3)-acetyltransferase [Campylobacter]
MKYLLEFENKKYSTYDFIETFYKLGLQKGDTLCVHTELFNFGFPLLSRNEFLQTILDCFFEVIGKEGTLIMPTFT
YSFCKNEVYDKINSKTKMGALNEYFRKQTGVKRTNDPIFSFAIKGAKEELFLKDTTSCFGENCVYEVLTKENGKY
>WP_002856595.1 MULTISPECIES: acetyl-CoA carboxylase biotin carboxylase subunit [Campylobacter]
MNQIHKILIANRAEIAVRVIRACRDLHIKSVAVFTEPDRECLHVKIADEAYRIGTDAIRGYLDVARIVEIAKACG

如果我有一个 ID，这个脚本就可以工作，但在某些情况下我有两个 ID，我会收到一个错误，因为我认为它正在寻找像 "WP_002855993.1 WP_002856586.1" 这样的 ID。有没有办法修改此脚本，使其查找两个单独的事件？我想这与 gawk 命令有关，但我不确定到底是什么。提前致谢！

Answer 1

考虑到您的输出文件是测试文件。

使用下面的命令只给你文件名：

>>cat text | awk '{print }' | grep -R 'WP*' | cut -d":" -f2

给我输出：

WP_002855993.1/5-168
WP_002856586.1/5-166
WP_002855993.1/5-168
WP_002856586.1/5-166

你想要这样的输出吗？

Answer 2

对原脚本的更新：

#!/usr/bin/env bash

for file_sto in *.sto; do
   file_faa=$(echo $file_sto | cut -d '_' -f 1,2,3)
   file_faa=${file_faa}"_protein.faa"

   awk '(NR==FNR) { match([=10=],/WP_.\{0,11\}/);
                    if (RSTART > 0)  a[substr([=10=],RSTART,RLENGTH)]++ 
                    next; }
        ( in a){ print RS [=10=] }' $file_sto RS=">" $file_faa >> sequence_protein.file
done

awk 部分甚至可以简化为：

awk '(NR==FNR) { if ([=11=] ~ /^WP_/) a[]++; next }
     ( in a) { print RS [=11=] }' FS='/' $file_sto FS=" " RS=">" $file_faa

此 awk 脚本执行以下操作：

将字段分隔符 FS 设置为 / 并读取文件 $file_sto.
读取$file_sto时记录号NR与文件记录号FNR相同。
(NR==FNR) { if ([=21=] ~ /^WP_/) a[]++; next }: 由于前面的条件，这一行只工作了一个$file_sto。它检查该行是否以 WP_ 开头。如果是，它将第一个字段 </code>（由 <code>FS 分隔，即 /）存储在数组 a 中；然后跳到文件中的下一条记录 (next)。
如果我们完成读取文件 $file_sto，我们将字段分隔符设置回单个 space FS=" "（参见 section Regular expression）和记录分隔符 RS 到 > 并开始读取文件 $file_faa 后者意味着 [=34=] 将包含 > 和第一个字段 </code> 之间的所有行是 <code>protID.
读取$file_faa，文件记录号FNR从1重新开始，而NR不复位。因此，第一个 awk 行被跳过。
( in a){ print RS [=42=] }如果第一个字段在数组a中，打印记录，记录分隔符在它前面。

修复原始脚本：

如果您想保留原始脚本，可以将 protID 存储在列表中，然后循环列表：

#!/bin/bash

for fileName in *.sto; do
    protID_list=( $(grep -o "WP_.\{0,11\}" $fileName | sort | uniq) )
    echo ${protID_list[@]}
    file=$(echo $fileName | cut -d '_' -f 1,2,3)
    file=$(echo $file'_protein.faa')
    echo $file 
    for protID in ${protID_list[@]}; do
       if [ -n "$protID" ]; then
          gawk "/^>/{N=0}/^.*$protID/{N=1} {if(N)print}" $file >> 
          sequence_protein.file
       fi
    done
done

使用 awk 提取两个单独的字符串

Using awk to extract two separate strings

unix

bash

awk

grep

bioinformatics