☰

将文件 a 的第一列与文件 b 的段落相匹配

tag-icon • tag-icon text-processing
将文件 a 的第一列与文件 b 的段落相匹配

我有2个文件。第一个，fileA看起来像

TCONS_00000066  XLOC_000030     -       u       q1:XLOC_000030|TCONS_00000066|0|0.000000|0.000000|0.000000|0.000000|-
TCONS_00000130  XLOC_000057     -       u       q1:XLOC_000057|TCONS_00000130|0|0.000000|0.000000|0.000000|0.000000|-
TCONS_00000395  XLOC_000206     -       u       q1:XLOC_000204|TCONS_00000393|0|0.000000|0.000000|0.000000|0.000000|-

FileB好像：

>TCONS_00000001 gene=XLOC_000001
AGATGAGCTGGTGGGGATGCTCTAAGAGAACGAGAGAAGCACAGAGCAGATAAACCACACCCACAGGCAC
CACCGTCCTTGTTGGTAATGAAGAAGACGAGACGACGACTTCCCCACTAGGAAACACGACGGAGGCGGAG
ATGATCGACGGCGGAGAGAGCTACAGAAACATCGATGCCTCCTGTCCAATCCCCCCATCCCATTCGGTAG
TTGGATTGAAGACTACCGAATAAGAGAAGCAGGCAGGCAGACAAACCCTTGAACCAAGGAGTCCTCGCTG
AGGAAGCTTTGGATCCACGACGCAGCTATGGCCTCCCCGCCCACCAGGCCGCCAGCCACAACCAGCTGAC
TAGGTCGCATGCATCATCAGATTTCAATCTCCCTTCGTTCCCTGTCCCTAATCCAATACCAATAGGGAGC
AATCAGCTGCTCCTCGACGGCGAGGGAGATGTCGTCGGCCGCGGGCCAAGACAACGGAGATACCGCTGGG
GACTACATCAAGTGGATGTGCGGCGCCGGTGGCCGTGCGGGCGGCGCCATGGCCAACCTCCAGCGCGGCG
TTGGCTCCCTCGTCCGTGACATTGGCGACCCCTGCCTCAACCCATCCCCCGTTAAGGGGAGCAAAATGCT
CAAACCGGAAAAATGGCACACATGTTTTGATAATGATGGAAAGGTCATAGGTTTCCGTAAAGCCCTAAAA
TTCATTGTCTTAGGGGGTGTGGATCCCACTATTCGAGCTGAAGTTTGGGAATTTCTTCTTGGCTGCTATG
CCTTGAGTAGTACCTCAGAGTATAGGAGGAAACTAAGAGCTGTTAGAAGGGAAAAATATCAAATTTTAGT
TAGACAGTGCCAGAGCATGCACCCAAGCATTGGTACAGGTGAGCTTGCTTACGCTGTTGGATCAAAGCTA

现在，fileA在第一列中包含选定的转录本编号，并fileB包含所有转录本的序列。我想扫描fileB第一列fileA并打印匹配转录本的尾随序列以及转录本编号。

答案1

sed -n 's|^\(TCONS_[0-9]*\) .*|\
         /^>\1 /,/\\n>/P|p' fileA |
sed -e '$!N' -f - -e D fileB >outfile

...这应该可以实现你所追求的。对于 fileA 中以字符串TCONS_和任何数字开头且后跟空格字符的每一行，第一sed行将打印如下行：

/^>TCONS_000001 /,/\n>/P

...其中000001是生成线的数字序列。

第二个sed给出了三个脚本 - 所有这些脚本都应用于其命名的输入文件 - fileB。

  1. 第一个命令$!N指示它将Next 输入行附加到模式空间中!除最后一行之外的每一行$。

  2. 接下来是 stdin - -f -- 第一个sed为其构造的，如前所述。

    • 第一个在第二个中打印的每一行sed都是一个 2 地址范围，//,//指示第二个在模式空间中为落在两个地址之间的每一行sed打印P到第一个ewline。\n

  3. 最后一个脚本只是D指示sed删除模式空间中D第一个出现的\newline 并再次尝试这三个脚本。

结果是，sed当块标题位于模式空间的头部时，第二个脚本的第一个脚本的任何范围的匹配都会开始- 将其拉入 w/ ext 并删除上一行^之后的迭代- 并在下一个行时结束块标题首先出现，并且仍然尾随由ewline 字符分隔的模式空间。因为第二个永远不会比该分隔符打印得更远，所以会以单行先行方式滑过 fileB - 打印以 fileA 第 1 列中找到的字符串开头的所有块，而不打印其他内容。ND\nsedPsed

其他（更复杂/更高效）方法

抛光：

also_fasta()( IFS='
';  get_match_lines(){
        tr   -s  '[:space:]'    \\n  |
        grep "$1" | grep -nFf - "$3"
    }
    get_script(){
        set \\ '\1,$p;n\' '1,\1b\1\' \
               '/^>/!b\1|p' '$c\' q
        sed -n "s|\([^:]*\).*|:\1$*"
    }
    get_match_lines "$@" < "$2"      |
    get_script | sed -nf - "$3"
)

如果您在 shell 中定义该函数并像这样调用它：

also_fasta 'TCONS_[0-9]*' fileA fileB

...它应该可以解决问题。

我在做了一些测试后想到了这一点。虽然第一个肯定有效，但它也绝对是慢的对于大输入。

将示例中的所有[ACGT]*行复制到剪贴板后，我创建了示例数据，例如：

seq -w -s " gene=XLOC_dummy
$(        xsel -bo | sed 's/ *$//')
>TCONS_"  0 512 99999999 >/tmp/temp

上面用于seq编写如下块：

>TCONS_99993600 gene=XLOC_dummy
AGATGAGCTGGTGGGGATGCTCTAAGAGAACGAGAGAAGCACAGAGCAGATAAACCACACCCACAGGCAC
CACCGTCCTTGTTGGTAATGAAGAAGACGAGACGACGACTTCCCCACTAGGAAACACGACGGAGGCGGAG
ATGATCGACGGCGGAGAGAGCTACAGAAACATCGATGCCTCCTGTCCAATCCCCCCATCCCATTCGGTAG
TTGGATTGAAGACTACCGAATAAGAGAAGCAGGCAGGCAGACAAACCCTTGAACCAAGGAGTCCTCGCTG
AGGAAGCTTTGGATCCACGACGCAGCTATGGCCTCCCCGCCCACCAGGCCGCCAGCCACAACCAGCTGAC
TAGGTCGCATGCATCATCAGATTTCAATCTCCCTTCGTTCCCTGTCCCTAATCCAATACCAATAGGGAGC
AATCAGCTGCTCCTCGACGGCGAGGGAGATGTCGTCGGCCGCGGGCCAAGACAACGGAGATACCGCTGGG
GACTACATCAAGTGGATGTGCGGCGCCGGTGGCCGTGCGGGCGGCGCCATGGCCAACCTCCAGCGCGGCG
TTGGCTCCCTCGTCCGTGACATTGGCGACCCCTGCCTCAACCCATCCCCCGTTAAGGGGAGCAAAATGCT
CAAACCGGAAAAATGGCACACATGTTTTGATAATGATGGAAAGGTCATAGGTTTCCGTAAAGCCCTAAAA
TTCATTGTCTTAGGGGGTGTGGATCCCACTATTCGAGCTGAAGTTTGGGAATTTCTTCTTGGCTGCTATG
CCTTGAGTAGTACCTCAGAGTATAGGAGGAAACTAAGAGCTGTTAGAAGGGAAAAATATCAAATTTTAGT
TAGACAGTGCCAGAGCATGCACCCAAGCATTGGTACAGGTGAGCTTGCTTACGCTGTTGGATCAAAGCTA

...到文件/tmp/temp，其中字符串的数字部分以每块 512 的间隔>TCONS_[0-9]*递增。整个文件大小约为 185mbs。0000051299999999

然后我做了：

grep -F 00\  </tmp/temp | cut -d\> -f2 >/tmp/tempA

...对于模式文件（这只是将选择范围缩小到任何TCONS_[0-9]*00<space>匹配项）。

两个文件的行数为：

wc -l /tmp/temp*
  2734369 /tmp/temp
     7813 /tmp/tempA
  2742182 total

虽然我没有sed针对甚至那个大小的输入运行两个 s - 我尝试的最大建议sed | sed是一个数据文件，一个该大小的什一税 - 以及一个选择的模式文件2\- 不过，它确实接近于此。

不管怎样，我思考了发生了什么，我意识到，当模式文件接近 8000 个正则表达式时，每次一行都会D它都必须通过一遍又一遍地检查之前出现的所有正则表达式来工作。这大概是不是最佳地完成。

至少在我生成的输入中，我确实有一件事情适合我——它或多或少是连续的。在该线程之后，我意识到如果我可以按行号而不是正则表达式工作，那么它甚至不必是连续的 - 所以我转向grep.

我运行的命令（以及我上面的功能的基础）是：

sed  -n 's|^\(TCONS_[0-9]*\) .*|>\1 |p
'    < /tmp/tempA              |
grep -nFf - /tmp/temp          |
sed  -n 's|\([^:]*\):.*|:\1\
        \1,$p;n;1,\1b\1\
        /^>/!b\1|p;$c\' -e q   |
sed  -nf - /tmp/temp

如果你使用这个，你应该替换fileBfor/tmp/temp和fileAfor /tmp/tempA。

grep -F使用固定字符串而不是正则表达式 - 并且速度要快得多（特别是考虑到我们正在与一线比赛一起工作）- 其余的基本上忽略了所有的grep结果，除了它在每个结果的开头返回的行号。sed中间的内容以这样的方式处理输出grep，即sed最终处理数据文件的程序将永远不会回溯其脚本，即使一次 - 它将在处理其输入时处理其脚本文件。

中间的每行打印的sed最后都会写下类似以下内容：sedgrep

:2732801                     #define branch label :LINENO
2732801,$p                   #if LINENO thru last line print
n                            #overwrite current line w/ next line
1,2732801b2732801            #if first line thru LINENO branch to :LINENO
/^>/!b2732801                #if !not /^>/ branch to :LINENO

因此，对于 的每一个grep匹配项sed都实现了两个微小的读取循环：第一个是一个空操作，其中用sed下一行覆盖当前行，直到将当前行号增加到 的grep下一个匹配项。第二个是打印循环，其中sed不断打印当前行，然后用下一行覆盖它，直到有一行匹配/^>/。通过这种方式，sed与其内文件同步处理其脚本 - 永远不会在脚本中前进到下一个匹配行号允许的范围之外。

这比其他脚本好几个数量级。它处理 185mbs 就像......

( also_fasta 'TCONS_[0-9]*' /tmp/tempA /tmp/temp; ) \
   3.93s user 0.39s system 100% cpu 4.281 total

...另一个处理输入大小的 10% ...

( sed -n 's|^\(TCONS_[0-9]*\) .*|/>\1 /,/\\n>/P|p' /tmp/tempA  |sed -e  -f... ) \
108.58s user 0.04s system 100% cpu 1:48.56 total

更具体地说，其他脚本输入的行数为：

wc -l /tmp/temp*
  273435 /tmp/temp
     782 /tmp/tempA
  274217 total

答案2

如果您要经常做这种事情，我建议您使用我拥有的几个工具，它们对于此类事情非常有用：

  1. FastaToTbl:

    #! /bin/sh
gawk '{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else 
            printf "%s\t", substr($0,2,length($0)-1)
        else 
                printf "%s", $0
}END{printf "\n"}'  "$@"

    将该文件另存为FastaToTblin~/bin并使其可执行chmod 755 ~/bin/FastaToTbl。

  2. TblToFasta

    #! /bin/sh
gawk '{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1

  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld
      if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
   while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}' "$@"

    将该文件另存为TblToFastain~/bin并使其可执行chmod 755 ~/bin/TblToFasta。

现在您已经设置了这两个脚本，您可以使用FastaToTbl将 FASTA 序列更改为 Tbl 格式（序列名称、TAB、序列），这样更容易搜索：

FastaToTbl fileB | grep -f <(awk '{print $1}' fileA)  | TblToFasta

它使用FastaToTbl脚本将序列和 ID 放在同一个文件上。打印andawk '{print $1}' fileA的第一个字段fileA，因为我们使用<()它可以作为“文件”传递给 grep 的-f选项，该选项要求要搜索的模式文件。最后TblToFasta返回FASTA格式。

最后，你可能想看看我的retrieveseqs.pl脚本可以执行此操作以及许多其他操作。例如，在这种情况下，您可以这样做：

retrieveseqs.pl -f fileB fileA

1感谢最初编写这两个脚本的 JF Abril。

相关内容