计算多个范围中唯一整数的数量以计算爆炸结果的水平覆盖范围[关闭]

Question

当我尝试将 DNA 序列读数映射/对齐到蛋白质数据库时，我找不到许多提供水平参考/主题序列覆盖的映射工具，但是blastx/diamond确实在每个查询读数上对齐了参考基因的开始和结束位置。因此，我需要计算每个参考基因的多个范围 [subjectStart..subjectEnd] 中唯一整数的数量，以获得比对序列的长度。

blastx/diamond 输出被修剪为 3 列，用于参考基因 ID、开始和结束，使用

awk -F"\t" '{if ($2 <= $3) {print $0} else {print $1"\t"$3"\t"$2}}'

修复反向序列。每个行 ID 可能有也可能没有彼此重叠或重复的范围，并且这些重复的数字只需要计算一次。

输入

ID  Start   End
A   1   50
A   2   45
A   25  150
A   50  150
A   155 200
A   205 300
B   5   50
B   61  70
B   81  100
C   1   500`

所需的输出。

ID  count
A   292
B   76
C   500

我在这个问题的第一版中请求帮助以获得所需的输出，我可以立即从许多专家那里得到完美的解决方案！ 以下内容非常适合我，我学到了很多东西。我衷心感谢大家！！！

Ruby by @dawg

ruby -lane 'BEGIN{h=Hash.new { |hash, key| hash[key] = Set.new() }}
    h[$F[0]].merge(($F[1].to_i..$F[2].to_i)) if $.>1
END{
    puts "ID\tCount"
    h.each{|k,v| puts "#{k}\t#{v.length}"}
}
' Input.file > Output.file

Perl 作者：@zdim

 perl -MData::Dumper -MList::Util=uniq -wnE'
    ($id, $beg, $end) = split; 
    next if not $beg or $beg =~ /[^0-9]/ or not $end or $end =~ /[^0-9]/;
    push @{$res{$id}}, $beg..$end; 
    }{ 
    for (keys %res) { $res{$_} = uniq sort { $a <=> $b } @{$res{$_}} };
    say Dumper \%res
' Input.file > Output.file

@pmf 的 Awk

 awk 'NR>1 {s[$1] += $3 - ($2 <= b[$1] ? ($3 <= b[$1] ? $3 : b[$1]) + 1 : $2) + 1; b[$1] = b[$1] <= $3 ? $3 : b[$1]} END {OFS="\t"; print "ID", "count"; for (i in s) {print i, s[i]}}' Input.file > Output.file

@dawg 的 awk

 awk '
FNR>1{for(i=$2;i<=$3;i++) ss[$1 "|" i]}
END{
    print "ID", "Count"
    for (e in ss) {
        split(e,idx,"\|")
        cnt[idx[1]]++
    }
    for (e in cnt) print e, cnt[e]
}
' OFS="\t" Input.file > Output.file

Answer 1

鉴于：

cat file
ID  Start   End
A   1   50
A   2   45
A   25  150
A   50  150
A   155 200
A   205 300
B   5   50
B   61  70
B   81  100
C   1   500

这里有一个 Ruby 可以做到这一点：

ruby -lane 'BEGIN{h=Hash.new { |hash, key| hash[key] = Set.new() }}
    h[$F[0]].merge(($F[1].to_i..$F[2].to_i)) if $.>1
END{
    puts "ID\tCount"
    h.each{|k,v| puts "#{k}\t#{v.length}"}
}
' file

或者这个 awk：

awk '
FNR>1{for(i=$2;i<=$3;i++) ss[$1 "|" i]}
END{
    print "ID", "Count"
    for (e in ss) {
        split(e,idx,"\|")
        cnt[idx[1]]++
    }
    for (e in cnt) print e, cnt[e]
}
' OFS="\t" file 

任一打印（但 awk 可能是无序的......）：

ID  Count
A   292
B   76
C   500

---

根据评论，B 50 5 而不是 B 5 50，我们可以让它工作吗？

ruby -lane 'BEGIN{h=Hash.new { |hash, key| hash[key] = Set.new() }}
    v1,v2=[$F[1],$F[2]].map(&:to_i)
    h[$F[0]].merge(v1..v2) if $.>1 && v2>v1
END{
    puts "ID\tCount"
    h.each{|k,v| puts "#{k}\t#{v.length}"}
}
' file

但这应该不是必要的。

..

运算符创建一个添加到唯一值集合中的范围。范围

50..5

是空值，因此不会添加任何内容。

这是在 IRB（Ruby 交互式 shell）中进行的测试：

irb(main):043:0> s=Set.new()
=> #<Set: {}>
irb(main):044:0> s.merge(50..5)
=> #<Set: {}>
irb(main):045:0> s.merge(1..5)
=> #<Set: {1, 2, 3, 4, 5}>

所以你可以看到没有添加任何东西。

Answer 2

单向

perl -MData::Dumper -MList::Util=uniq -wnE'
    ($id, $beg, $end) = split; 
    next if not $beg or $beg =~ /[^0-9]/ or not $end or $end =~ /[^0-9]/;
    push @{$res{$id}}, $beg..$end; 
    }{ 
    for (keys %res) { @{$res{$_}} = uniq sort { $a <=> $b } @{$res{$_}} };
    say Dumper \%res
' data.txt

此命令行程序（“one-liner”）是一个演示程序，请在正常程序中重写。

语法

}{

标记

END

块的开始——一旦文件中的所有行都已处理完毕，代码就会运行。可以在这里写

END { ... }

。

这会为每个 ID 生成一个唯一且已排序的整数列表。如果您特别需要计数，只需将列表分配给标量即可

for (keys %res) { $res{$_} = uniq sort { $a <=> $b } @{$res{$_}} };
   

注释

• 检查字段是否为数字是使用正则表达式完成的，因为它们被假定为整数。否则，最好使用

  Scalar::Util

  中的 looks_like_number

• 在较新的 Perls 中，可以使用 postfix dereference 语法

  来完成取消引用

  push $res{$_}->@*, $beg..$end;
  

• 使用 Key::Sort 排序更容易（特别是在更复杂的情况下）

  use Key::Sort qw(usort);  # use "isort" if negative ints are possible
  
  foreach my $id (sort keys %res) {
      $res{$id}->@* = uniq usort $res{$id}->@*;
  }
  

  或者，获取每个 ID 的计数而不是数字

      $res{$id} = uniq usort $res{$id}->@*;
  

  还要注意库的

  _inplace

  方法，如果没有

  uniq

  ，这里会很好

• 我首先

  sort

  然后

  uniq

  假设

  uniq

  本身在内部使用排序（因为

  uniq

  “保留唯一元素的顺序”这会更快）；但如果

  uniq

  不 排序，那么首先过滤掉重复项然后排序会更快。