提取包含两个模式的行答案

【问题标题】：Extract lines containing two patterns提取包含两个模式的行
【发布时间】：2018-12-06 22:28:06
【问题描述】：

我有一个文件，其中包含如下几行：

>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>
>header3
<pattern_1>ATGGCCACCAACAACCAGAGCTCCC
>header4
GACCGGCACGTACAACCTCCAGGAAATCGTGCCCGGCAGCGTGTGGATGGAGAGGGACGTG
>header5
TGCCCCCACGACCGGCACGTACAAC<pattern_2>

我想提取包含标题行的所有行。

我尝试过使用 grep，但它只提取序列行而不提取标题行。

grep <pattern_1> | grep <pattern_2> input.fasta > output.fasta

如何在 Linux 中提取同时包含模式和标题的行？图案可以出现在线条中的任何位置。不限于行的开头或结尾。

预期输出：

>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>

【问题讨论】：

我想提取所有包含标题行的行。对我来说毫无意义。请为该输入添加预期输出（或更具体）。
也许是另一种选择？ grep -E "pattern|header" file
注意fasta文件可以有多行序列。

标签： linux awk grep fasta

【解决方案1】：

$ grep -A 1 header[12] file
>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>

man grep:

   -A NUM, --after-context=NUM
          Print  NUM  lines  of  trailing  context  after  matching lines.
          Places  a  line  containing  a  group  separator  (--)   between
          contiguous  groups  of  matches.  With the -o or --only-matching
          option, this has no effect and a warning is given.

   -B NUM, --before-context=NUM
          Print NUM  lines  of  leading  context  before  matching  lines.
          Places   a  line  containing  a  group  separator  (--)  between
          contiguous groups of matches.  With the  -o  or  --only-matching
          option, this has no effect and a warning is given.

grep -B 1 pattern_[12]也可以工作，但是样本数据中有多个pattern_1s，所以...这次不行。

【讨论】：

【解决方案2】：

您可以像这样使用 awk 轻松做到这一点：

awk '/^>/{h=$0;next}
     /<pattern_1>/&&/<pattern_2>/{print h;print}' input.fasta > output.fasta

这是一个 sed 解决方案，它也可以产生所需的输出：

sed -n '/^>/{N;/<pattern_1>/{/<pattern_2>/p}}' input.fasta > output.fasta

如果可能存在多行记录，您可以使用：

awk -v pat1='<pattern_1>' -v pat2='<pattern_2>' '
/^>/ {r=$0;p=0;next}
!p {r=r ORS $0;if(chk()){print r;p=1};next}
p

function chk(   tmp){
    tmp=gensub(/\n/,"","g",r)
    return (tmp~pat1&&tmp~pat2)
}' input.fasta > output.fasta

【讨论】：

多记录解决方案无效。 fasta 文件中的多行代表一个序列，只是分成多行以便于阅读（80 个字符宽），我相信在这种情况下和应该在完整的多行序列中的某个位置。
解决方案仍然无效，因为pattern1 或pattern2 可以拆分为多行。
@kvantour 好吧，我放弃了
@kvantour 我喝醉了。不过我稍后再试
@kvantour 再次更新。不能再进一步了。

【解决方案3】：

您可能对BioAwk 感兴趣，它是 awk 的改编版本，用于处理 fasta 文件

bioawk -c fastx -v seq1="pattern1" -v seq2="pattern2" \
       '($seq ~ seq1) && ($seq ~ seq2) { print ">"$name; print $seq }' file.fasta

如果要seq1开头，seq2结尾，可以改成：

bioawk -c fastx -v seq1="pattern1" -v seq2="pattern2" \
       '($seq ~ "^"seq1) && ($seq ~ seq2"$") { print ">"$name; print $seq }' file.fasta

这对于处理 fasta 文件非常实用，因为序列通常分布在多行中。上面的代码很容易处理这个问题，因为变量 $seq 包含完整的序列。

如果您不想安装 BioAwk，可以使用以下方法处理您的 FASTA 文件。它将允许多行序列并执行以下操作：

一次读取一条记录（假设标题中没有>，第一个字符除外）
从记录中提取标题并将其存储在name中（不是真的需要）
将完整序列合并到单个字符串中，删除所有换行符和空格。这样可以确保在模式拆分为多行时搜索 pattern1 或 pattern2 不会失败。
如果找到匹配项，则打印记录。

以下 awk 执行请求：

awk -v seq1="pattern1" -v seq2="pattern2" \
    'BEGIN{RS=">"; ORS=""; FS="\n"}
     { seq="";for(i=2;i<=NF;++i) seq=seq""$i; gsub(/[^a-zA-Z0-9]/,"",seq) }
     (seq ~ seq1 && seq ~ seq2){print ">" $0}' file.fasta

如果记录头包含其他不在行首的> 字符，您必须采取稍微不同的方法（除非您使用 GNU awk）

awk -v seq1="pattern1" -v seq2="pattern2" \
    '/^>/ && (seq ~ seq1 && seq ~ seq2) {
         print name
         for(i=0;i<n;i++) print aseq[i]
     }
     /^>/ { seq=""; delete aseq; n=0; name=$0; next }
     { aseq[n++] = $0; seq=seq""$0; sub(/[^a-zA-Z0-9]*$/,"",seq) }
     END { if (seq ~ seq1 && seq ~ seq2) {
              print name
              for(i=0;i<n;i++) print aseq[i]
            }
     }' file.fasta

注意：我们在这里使用sub，以防在fasta文件中引入意外字符（例如空格/制表符或CR (\r)）

注意： BioAwk 基于Brian Kernighan's awk，记录在"The AWK Programming Language", by Al Aho, Brian Kernighan, and Peter Weinberger (Addison-Wesley, 1988, ISBN 0-201-07981-X) 中。我不确定这个版本是否兼容POSIX。

【讨论】：

克万图尔。 Tks 分享 Awk 书。
@OXXO 请注意，这是一本非常古老的书，仅供参考。更好的书籍可以在stackoverflow.com/tags/awk/info 上找到
感谢您的建议。克万图尔

【解决方案4】：

如果您的输入文件与您的帖子中描述的完全一样，那么您可以使用：

grep -B1 '^<pattern_1>.*<pattern_2>$' input 
>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>

-B1 将在匹配行的顶部显示它之前的行。使用的正则表达式基于您的 2 个模式在行首和行尾以确切顺序定位的假设。如果不是这种情况：使用'.*<pattern_1>.*<pattern_2>.*'。最后但同样重要的是，如果不总是遵守 2 种模式的顺序，那么您可以使用：'^.*<pattern_1>.*<pattern_2>.*$\|^.*<pattern_2>.*<pattern_1>.*$'

在以下输入文件上：

cat input
>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>
>header2b
<pattern_2>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_1>
>header3
<pattern_1>ATGGCCACCAACAACCAGAGCTCCC
>header4
GACCGGCACGTACAACCTCCAGGAAATCGTGCCCGGCAGCGTGTGGATGGAGAGGGACGTG
>header5
TGCCCCCACGACCGGCACGTACAAC<pattern_2>

输出：

grep -B1 '^.*<pattern_1>.*<pattern_2>.*$\|^.*<pattern_2>.*<pattern_1>.*$' input 
>header1
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCGGGCCTCTTTTCCTGACGGCCGCCCCCACTGCCCCCACGACCGGCCCGTACAAC<pattern_2>
>header2
<pattern_1>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_2>
>header2b
<pattern_2>CGGCGGGCAGATGGCCACCAACAACCAGAGCTCCCTGGCCTGCAATCACTACTCGTGTTTTGCCACCACTGCCCCCACGACCGGCACGTACAAC<pattern_1>

【讨论】：

【解决方案5】：

如果您希望 grep 打印匹配周围的行，请使用 -B 标志表示之前的行，使用 -A 表示之后的行，使用 -C 表示匹配之前和之后的行。

在您的情况下， grep -B 1 似乎可以完成这项工作。

【讨论】：