【问题标题】:Filtering blast results by evalue but only if unique按 evalue 过滤爆炸结果,但前提是唯一的
【发布时间】:2017-04-07 03:55:47
【问题描述】:

我正在开发一个管道,该管道在某些时候会生成数百个不同格式的文件(我在我不关心的字段中写 X):

 id1   X   X   X   X   X   X   X   X   X  evalue1   X
 id2   X   X   X   X   X   X   X   X   X  evalue2   X     
 ...

我必须过滤此文件,以便对于每个 ID,根据 evalue 获取最佳结果(越小越好),但如果最佳 evalue 使用相同的 ID 重复,则不计算该 ID。

例如,如果输入文件是:

 id1   X   X   X   X   X   X   X   X   X  3e-07   X
 id1   X   X   X   X   X   X   X   X   X  3e-04   X
 id2   X   X   X   X   X   X   X   X   X  3e-07   X     
 id3   X   X   X   X   X   X   X   X   X  3e-04   X     
 id3   X   X   X   X   X   X   X   X   X  3e-04   X     
 id3   X   X   X   X   X   X   X   X   X  1e-02   X     

预期的输出是:

 id1   X   X   X   X   X   X   X   X   X  3e-07   X
 id2   X   X   X   X   X   X   X   X   X  3e-07   X     

在 id1 的两次命中之间,最差的被删除,并且由于 id3 的最佳 evalue 不是唯一的,因此不存储 ID。

我尝试过调整 blast 命令行工具,但最接近的选项是将最大命中数设置为 1,但输出中仍保留 id3 之类的情况。所以我的解决方案是一个 python 脚本,但是文件的数量使这个过程非常耗时。

有没有办法使用 bash 工具(awk?)过滤这些文件以提高效率?

每个文件都有唯一的标识符,因此同一 ID 不能出现在多个文件中。

提前致谢

更新 1:

这里是示例文件:

 D00733:159:CA65UANXX:8:1104:7340:77245  gi|13507739|ref|NC_000912.1|    100.00  24      0       0       1       24      529212  529189  3e-07   44.6
 D00733:159:CA65UANXX:8:2303:18019:72377 gi|13507739|ref|NC_000912.1|    100.00  20      0       0       1       20      622755  622736  2e-05   37.4
 D00733:159:CA65UANXX:8:2103:11030:25200 gi|13507739|ref|NC_000912.1|    95.24   21      1       0       1       21      321813  321833  3e-04   33.7
 D00733:159:CA65UANXX:8:2103:11030:25200 gi|13507739|ref|NC_000912.1|    95.24   21      1       0       1       21      495963  495943  3e-04   33.7
 D00733:159:CA65UANXX:8:2103:11030:25200 gi|13507739|ref|NC_000912.1|    95.00   20      1       0       2       21      613871  613852  0.001   31.9

使用@karafka 建议的解决方案后,输出为:

 D00733:159:CA65UANXX:8:2303:18019:72377 gi|13507739|ref|NC_000912.1|   100.00  20  0   0   1   20  622755  622736  2e-05   37.4
 D00733:159:CA65UANXX:8:2103:11030:25200 gi|13507739|ref|NC_000912.1|   95.00   20  1   0   2   21  613871  613852  0.001   31.9
 D00733:159:CA65UANXX:8:1104:7340:77245  gi|13507739|ref|NC_000912.1|   100.00  24  0   0   1   24  529212  529189  3e-07   44.6

似乎最后一个 id 以 0.001 为最小值。

我正在使用 GNU Awk 3.1.5

更新 2

强制数值转换并不能解决 awk 3.1.5 中的问题,唯一的解决方案是:将 awk 更新到 >= 3.1.8

【问题讨论】:

  • 无法复制您的结果。我仍然得到两个记录。也许旧版本有一些特殊性。您可以通过添加 0 来强制进行数字转换。
  • 没问题 ;) 感谢您的宝贵时间
  • 在 GNU Awk 3.1.8 中按预期工作,但在 3.1.5 版本中输出错误...
  • 升级时间)

标签: bash awk blast


【解决方案1】:

awk 来救援!

awk '!($1 in min) || $11<min[$1] {min[$1]=$11; line[$1]=$0} 
     END {for(k in line) print line[k]}' file

 id1   X   X   X   X   X   X   X   X   X  3e-07   X
 id2   X   X   X   X   X   X   X   X   X  3e-07   X
 id3   X   X   X   X   X   X   X   X   X  3e-04   X

这不依赖于条目的顺序,也不能保证输出顺序。

sort 协助的另一种解决方案

sort -k1,1 -k11g file | awk '!a[$1]++'

 id1   X   X   X   X   X   X   X   X   X  3e-07   X
 id2   X   X   X   X   X   X   X   X   X  3e-07   X
 id3   X   X   X   X   X   X   X   X   X  3e-04   X

仅在最小值唯一时打印

awk '!($1 in min) || $11<=min[$1] {min[$1]=$11; line[$1]=$0; c[$1,$11]++}
    END {for(k in line) if(c[k,min[k]]==1) print line[k]}' file

 id1   X   X   X   X   X   X   X   X   X  3e-07   X
 id2   X   X   X   X   X   X   X   X   X  3e-07   X

要强制进行数字转换,您可以将 0 添加到值 ($11)。例如

... $11+0<=min[$1] {min[$1]=$11+0; line[$1]=$0; c[$1,$11+0]++}...

【讨论】:

  • 不错!这两个选项都满足取最佳值的第一个条件,但在 id3 的情况下,它不应该返回任何值,因为最佳值不是唯一的,有什么方法可以修改它?
  • 这就是我需要的!谢谢
  • 为什么如果其中一个数字没有'e'浮点表示法,它会将没有它的值作为最小值?例如,如果我们添加一个 evalue 为 0.01 的额外 id3,它会返回:id3 X X X X X X X X X 0.01 X
  • 没关系,我们是在比较数字,无论格式如何。
  • 这就是我的想法,但如果我的输入是:``` id1 X X X X X X X X X 3e-07 X id2 X X X X X X X X X 3e-07 X id3 X X X X X X X X 0.001 X id3 X X X X X X X X X 3e-04 X id3 X X X X X X X X X 3e-04 X ` `` 最后的解决方案返回:``` id1 X X X X X X X X X 3e-07 X id2 X X X X X X X X X 3e-07 X id3 X X X X X X X X X 0.001 X ```
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2019-07-24
  • 2019-04-24
  • 2021-09-10
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-10-05
相关资源
最近更新 更多